source: trunk/MagicSoft/Mars/mhist/MHMcCollectionArea.cc@ 1014

/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de)
!              Thomas Bretz  12/2000 (tbretz@uni-sw.gwdg.de)
!
!   Copyright: MAGIC Software Development, 2000-2001
!
!
\* ======================================================================== */

#include "MHMcCollectionArea.h" 

#include <TH2.h>
#include <TCanvas.h>

#include "MH.h"

ClassImp(MHMcCollectionArea);

// --------------------------------------------------------------------------
//
//  Creates the three necessary histograms:
//   - selected showers (input)
//   - all showers (input)
//   - collection area (result)
//
MHMcCollectionArea::MHMcCollectionArea(const char *name, const char *title)
{ 
    //   initialize the histogram for the distribution r vs E
    //
    //   we set the energy range from 1 Gev to 10000 GeV (in log 5 orders
    //   of magnitude) and for each order we take 10 subdivision --> 50 xbins
    //
    //   we set the radius range from 0 m to 500 m with 10 m bin --> 50 ybins

  
    fName  = name  ? name  : "MHMcCollectionArea";
    fTitle = title ? title : "Data to Calculate Collection Area";

    const Int_t  nbins = 50;
    const Float_t maxx = 5;
    const Float_t maxy = 500;

    Float_t *binsx = new Float_t[nbins+1];
    for (int i=0; i<nbins+1; i++)
        binsx[i] = pow(10, maxx*i/nbins);

    Float_t *binsy = new Float_t[nbins+1];
    for (int i=0; i<nbins+1; i++)
        binsy[i] = maxy*i/nbins;

    fHistAll = new TH2D("AllEvents", "All showers - Radius vs Energy distribution",
                        nbins, binsx, nbins, binsy);
    fHistSel = new TH2D("SelectedEvents", "Selected showers - Radius vs Energy distribution",
                        nbins, binsx, nbins, binsy);
    fHistCol = new TH1D("CollectionArea", "Collection Area vs. Energy",
                        nbins, binsx);

    delete binsx;
    delete binsy;

    fHistAll->SetDirectory(NULL);
    fHistSel->SetDirectory(NULL);
    fHistCol->SetDirectory(NULL);

    fHistAll->SetXTitle("E [GeV]");
    fHistAll->SetYTitle("r [m]");
    fHistAll->SetZTitle("N");

    fHistSel->SetXTitle("E [GeV]");
    fHistSel->SetYTitle("r [m]");
    fHistSel->SetYTitle("N");

    fHistCol->SetXTitle("E [GeV]");
    fHistCol->SetYTitle("A [m^{2}]");
}

// --------------------------------------------------------------------------
//
// Delete the three histograms
//
MHMcCollectionArea::~MHMcCollectionArea()
{
    delete fHistAll;
    delete fHistSel;
    delete fHistCol;
}

// --------------------------------------------------------------------------
//
// Fill data into the histogram which contains all showers
//
void MHMcCollectionArea::FillAll(Float_t energy, Float_t radius)
{
    fHistAll->Fill(energy, radius);
}

// --------------------------------------------------------------------------
//
// Fill data into the histogram which contains the selected showers
//
void MHMcCollectionArea::FillSel(Float_t energy, Float_t radius)
{
    fHistSel->Fill(energy, radius);
}

// --------------------------------------------------------------------------
//
// Draw the histogram with all showers
//
void MHMcCollectionArea::DrawAll(Option_t* option)
{
    if (!gPad)
        MH::MakeDefCanvas(fHistAll);

    fHistAll->Draw(option);

    gPad->SetLogx();

    gPad->Modified();
    gPad->Update();
}

// --------------------------------------------------------------------------
//
// Draw the histogram with the selected showers only.
//
void MHMcCollectionArea::DrawSel(Option_t* option)
{
    if (!gPad)
        MH::MakeDefCanvas(fHistSel);

    fHistSel->Draw(option);

    gPad->SetLogx();

    gPad->Modified();
    gPad->Update();
}

// --------------------------------------------------------------------------
//
// Creates a new canvas and draws the histogram into it.
// Be careful: The histogram belongs to this object and won't get deleted
// together with the canvas.
//
TObject *MHMcCollectionArea::DrawClone(Option_t* option) const
{
    TCanvas *c = MH::MakeDefCanvas(fHistCol);

    //
    // This is necessary to get the expected bahviour of DrawClone
    //
    gROOT->SetSelectedPad(NULL);

    fHistCol->DrawCopy(option);

    gPad->SetLogx();

    c->Modified();
    c->Update();

    return c;
}

void MHMcCollectionArea::Draw(Option_t* option)
{
    if (!gPad)
        MH::MakeDefCanvas(fHistCol);

    fHistCol->Draw(option);

    gPad->SetLogx();

    gPad->Modified();
    gPad->Update();
}

// --------------------------------------------------------------------------
//
//  Calculate the Efficiency (collection area)
//
void MHMcCollectionArea::CalcEfficiency()
{
    // Description!

    //
    //  first of all calculate the efficency
    //  do it here by hand to get the right error of efficency
    //
    const Int_t nbinx = fHistSel->GetXaxis()->GetNbins();
    const Int_t nbiny = fHistSel->GetYaxis()->GetNbins();

    for (Int_t ix=1; ix<=nbinx; ix++)
    {
        for (Int_t iy=1; iy<=nbiny; iy++)
        {
            const Float_t Na = fHistAll->GetCellContent(ix, iy);

            if (Na <= 0)
                continue;

            const Float_t Ns = fHistSel->GetCellContent(ix, iy);

            const Double_t eff = Ns/Na;
            const Double_t err = sqrt((1.-eff)*Ns)/Na;

            // old calculation from Harald:
            //  const Double_t eff = Ns/Na;
            //  const Double_t err = sqrt(Na + Na*Ns - Ns*Ns - Ns) / (Na*Na);

            fHistSel->SetCellContent(ix, iy, eff);
            fHistSel->SetCellError(ix, iy, err);
        }
    }

    //
    //  now calculate the Collection area for different
    //  energies
    //
    for (Int_t ix=1; ix<=nbinx; ix++)
    {
        Double_t errA = 0;
        Double_t colA = 0;

        for (Int_t iy=1; iy<=nbiny; iy++)
        {
            TAxis *yaxis = fHistSel->GetYaxis();

            const Double_t r1  = yaxis->GetBinLowEdge(iy);
            const Double_t r2  = yaxis->GetBinLowEdge(iy+1);

            const Double_t A   = TMath::Pi() * (r2*r2 - r1*r1);

            const Double_t eff = fHistSel->GetCellContent(ix, iy);
            const Double_t err = fHistSel->GetCellError(ix, iy);

            colA += eff*A;
            errA += A*A * err*err;
        }

        fHistCol->SetBinContent(ix, colA);
        fHistCol->SetBinError(ix, sqrt(errA));
    }
}