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

#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. ) ; 

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

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

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) 
{ 
  fHistColl->Draw(option) ; 
} 

void MCollArea::CalculateEffi()
{ 
  //  first of all calculate the efficency
  //  do it here by hand to get the right error of efficency
 
  Int_t iBinx = ( (TAxis *) fHistSel->GetXaxis())->GetNbins() ; 
  Int_t iBiny = ( (TAxis *) fHistSel->GetYaxis())->GetNbins() ; 
 
  Float_t N, Nall ; 
  Double_t effi, error ; 
  for (Int_t ix=1; ix<=iBiny; ix++ ) 
    {
      for (Int_t iy=1; iy<=iBiny; iy++ ) 
        { 
          effi = error = 0. ; 

          N    = fHistSel->GetCellContent(ix, iy) ; 
          Nall = fHistAll->GetCellContent(ix, iy) ; 
          
          if ( Nall > 0 ) { 
            effi   = N / Nall ; 
            error = sqrt ( Nall + Nall * N - N * N - N ) / (Nall * Nall ) ; 

            cout << ix << " " << iy 
                 << " N " << N 
                 << " Nall " << Nall
                 << " effi  " << effi 
                 << " error " << error 
                 << endl ; 

            fHistSel->SetCellContent(ix, iy, effi) ; 
            fHistSel->SetCellError(ix, iy, error) ; 
          } 
          else 
            cout << ix << " " << iy << endl ; 
            
          

          
        }
    } 
  
  // 
  //  now calculate the Collection area for different 
  //  energies
  //   
  
  
  Double_t  r1, r2, eff, errEff, A, errA, collA ; 

  for (Int_t ix=1; ix<=iBiny; ix++ ) 
    {
      A = errA = collA = errEff = 0. ; 
      
      for (Int_t iy=1; iy<=iBiny; iy++ ) 
        { 
          r1 = ( (TAxis *) fHistSel->GetYaxis())->GetBinLowEdge(iy) ; 
          r2 = ( (TAxis *) fHistSel->GetYaxis())->GetBinLowEdge(iy+1) ; 
          A  = 3.141592654 * ( r2*r2 - r1*r1 ) ; 
          eff    = fHistSel->GetCellContent(ix, iy) ;
          errEff = fHistSel->GetCellError(ix, iy) ; 
          collA += eff * A ; 
          
          errA  += ( A * A ) * errEff * errEff ; 
        }

      errA = sqrt( errA ) ; 

      cout << ix << " --> " << A
           << endl ; 
      
      fHistColl->SetBinContent(ix, collA ) ; 
      fHistColl->SetBinError(ix, errA ) ; 
      

    } 
}