source: trunk/MagicSoft/Mars/mmontecarlo/MCollArea.cc@ 745

#include "MCollArea.h" 

#include <MLog.h>
#include <TH2.h> 

ClassImp(MCollArea)



MCollArea::MCollArea(const char *name, const char *title) 
{ 
  //
  //   default constructor
  //

  //   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  : "MCollArea";
  *fTitle = title ? title : "Data to Calculate Coll-Area";
  
  
  fHistAll = new TH2D("collAreaAll", "all showers - Radius vs log(E) distribution", 
                      50, 0., 5., 
                      50, 0., 500. ) ;
  
  fHistSel = new TH2D("collAreaSel", "selected showers - Radius vs log(E) distribution", 
                      50, 0., 5., 
                      50, 0., 500. ) ; 

  fHistCol = new TH1D("collArea", "Collection Area",
                      50, 0., 5.) ;
  
} 

MCollArea::~MCollArea()
{ 
  delete fHistAll ; 
  delete fHistSel ; 
  delete fHistCol ;
} 

void MCollArea::FillAll(Float_t log10E, Float_t radius)
{ 
  fHistAll->Fill(log10E, radius ) ; 
} 

void MCollArea::FillSel(Float_t log10E, Float_t radius)
{ 
  fHistSel->Fill(log10E, radius ) ; 
} 

void MCollArea::DrawAll() 
{ 
  fHistAll->Draw() ; 
} 

void MCollArea::DrawSel() 
{ 
  fHistSel->Draw() ; 
} 

void MCollArea::Draw(Option_t* option) 
{ 
  fHistCol->Draw(option) ;
} 

void MCollArea::CalcEfficiency()
{ 
    // Description!

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

    for (Int_t ix=1; ix<=iBiny; ix++ )
    {
        for (Int_t iy=1; iy<=iBiny; iy++ )
        {
            const Float_t N    = fHistSel->GetCellContent(ix, iy);
            const Float_t Nall = fHistAll->GetCellContent(ix, iy);

            if ( Nall <= 0 ) {
                // cout << ix << " " << iy << endl ;
                continue;
            }

            const Double_t eff = N / Nall ;
            const Double_t err = sqrt(Nall + Nall*N - N*N - N) / (Nall*Nall);
            /*
             cout << ix << " " << iy
             << " N " << N
             << " Nall " << Nall
             << " effi  " << eff
             << " error " << err
             << endl ;
             */
            fHistSel->SetCellContent(ix, iy, eff);
            fHistSel->SetCellError(ix, iy, err);
        }
    }

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

        for (Int_t iy=1; iy<=iBiny; iy++ )
        {
            const Double_t r1  = ((TAxis *)fHistSel->GetYaxis())->GetBinLowEdge(iy);
            const Double_t r2  = ((TAxis *)fHistSel->GetYaxis())->GetBinLowEdge(iy+1);

            const Double_t A   = kPI * (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));
    }
}