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

    } 
} 
