Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc	(revision 6368)
+++ trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc	(revision 6465)
@@ -38,4 +38,5 @@
 #include "MMcEvt.hxx"
 #include "TH2F.h"
+#include "TFile.h"
 
 #include "MLog.h"
@@ -55,29 +56,27 @@
 //
 MEffAreaAndCoeffCalc::MEffAreaAndCoeffCalc()
-  : fSpec(NULL), fHorig(NULL), fEbins(10), fEsubbins(20), fEmin(10.), fEmax(10000.),
-    fCoeff(NULL), fEffA(NULL)
+  : fSpec(NULL), fHorig(NULL), fEmin(10.), fEmax(10000.), fEbins(10), fEsubbins(20),
+    fCoeff(NULL), fEffA(NULL), fFile(NULL)
 {
   // set the default function
-  // OJO!! this should be changed
   SetFunction("4.e9*pow(x,-2.6+1)");  
 
   // create the TChains 
   // OJO!! make sure they read the appropriate branch
-  fCini = new TChain("Events");
+  fCini = new TChain("OriginalMC");
   fCcut = new TChain("Events");
 
-  // define some usefull aliases
-  fCini->SetAlias("logenergy","log10(MMcEvt.fEnergy)");
-  fCini->SetAlias("theta","MMcEvt.fTheta*180./3.14159");
-
-  // OJO!! the estimated energy must be changed for a real one
+  // define some useful aliases
+  fCini->SetAlias("logenergy","log10(MMcEvtBasic.fEnergy)");
+  fCini->SetAlias("theta","MMcEvtBasic.fTelescopeTheta*180./3.14159");
+
   fCcut->SetAlias("logenergy","log10(MMcEvt.fEnergy)");
-  fCcut->SetAlias("theta","MMcEvt.fTheta*180./3.14159");
-
-  fCcut->SetBranchStatus("*",0);
-  fCcut->SetBranchStatus("MHillas.*",1);
-  fCcut->SetBranchStatus("MMcEvt.*",1);
+  fCcut->SetAlias("theta","MMcEvt.fTelescopeTheta*180./3.14159");
+
   fCcut->SetBranchAddress("MHillas.",&fHillas);
   fCcut->SetBranchAddress("MMcEvt.",&fMcEvt);
+  
+  // borra
+  fFile = new TFile("test.root","RECREATE");
 }
 
@@ -119,4 +118,9 @@
       emin = fEmin;
       emax = fEmax;
+    }
+  else
+    {
+      fEmin = emin;
+      fEmax = emax;
     }
   fSpec = new TF1("fspec",chfunc,emin,emax);
@@ -147,4 +151,6 @@
   fHorig   = new TH1F("fhorig","Original energy spectrum",esbins,logemin,logemax);
   fCini->Draw("logenergy>>fhorig","","goff");
+  // borra
+  fHorig->Write();
 }
 
@@ -169,4 +175,6 @@
     FillOriginalSpectrum();
   
+  // borra
+  TH1F hw("hw","hw",esbins,fEmin,fEmax);
   for(Int_t i=0;i<esbins;i++)
     {
@@ -181,5 +189,9 @@
 	  fWeight[i]=-1;
 	}
-    }
+      // borra
+      hw.SetBinContent(i+1,fWeight[i]);
+    }
+  // borra
+  hw.Write();
 }
 
@@ -207,4 +219,9 @@
   TH2F* hmc  = new TH2F("hmc","MC energy",fEbins,logemin,logemax,fNTbins,fTbin);
 
+  // borra
+  TH1F* hest1  = new TH1F("hest1","Estimated energy",fEbins,logemin,logemax);
+  TH1F* hmc1   = new TH1F("hmc1","MC energy",fEbins,logemin,logemax);
+  TH1F* coef1  = new TH1F("coef1","coefficients",fEbins,logemin,logemax);
+
   // reset the coefficients
   if(fCoeff)
@@ -220,5 +237,5 @@
       // OJO!! Estimated energy will be taken directly from some input container
       Float_t emc   = fMcEvt->GetEnergy();
-      Float_t estim = fHillas->GetSize()/15.;
+      Float_t estim = fHillas->GetSize()/0.18/15.;
       Float_t theta = fMcEvt->GetTheta()*180./3.14159; // zenith angle (in degrees)
 
@@ -231,7 +248,13 @@
 	  hest->Fill(TMath::Log10(estim),theta,fWeight[wbin]); 
 	  hmc->Fill(TMath::Log10(emc),theta,fWeight[wbin]);
+	  // borra
+	  if(theta<fTbin[1])
+	    {
+	      hest1->Fill(TMath::Log10(estim),fWeight[wbin]); 
+	      hmc1->Fill(TMath::Log10(emc),fWeight[wbin]);
+	    }
 	}
       else
-	cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has energy out of bounds, skpping" << endl;
+	cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has energy out of bounds, skipping" << endl;
     }
 
@@ -242,6 +265,14 @@
 	const Float_t num = hmc->GetBinContent(i+1,j+1);
 	const Float_t den = hest->GetBinContent(i+1,j+1);
+	//borra
+	const Float_t num1 = hmc1->GetBinContent(i+1);
+	const Float_t den1 = hest1->GetBinContent(i+1);
 	if(den)
-	  fCoeff->SetBinContent(i+1,j+1,num/den);
+	  {
+	    fCoeff->SetBinContent(i+1,j+1,num/den);
+	    //borra
+	    if(j==0 && den1)
+	      coef1->SetBinContent(i+1,num1/den1);
+	  }
 	else
 	  {
@@ -252,4 +283,9 @@
 
 
+  //borra
+  hmc1->Write();
+  hest1->Write();
+  coef1->Write();
+
   delete hmc;
   delete hest;
@@ -270,14 +306,19 @@
   //OJO!! radius should be read from somewhere!
   const Double_t radius = 30000.; // generation radius (in cm)
-  const Double_t Atot = 3.14159*radius*radius; //total area (in cm^2)
+  //  const Double_t Atot = 3.14159*radius*radius; //total area (in cm^2)
+  const Double_t Atot = 1.; //total area (in cm^2)
   const Double_t logemin = TMath::Log10(fEmin);
   const Double_t logemax = TMath::Log10(fEmax);
   const Int_t esbins = fEbins*fEsubbins; // total number of subbins
+  const Double_t de = (logemax-logemin)/fEbins; // bin size (in log)
+  const Double_t desub = de/fEsubbins; // subbin size (in log)
 
   // reset the effective areas
   if(fEffA)
     delete fEffA;
-  fEffA = new TH2F("feffa","Effective area",esbins,logemin,logemax,fNTbins,fTbin);
-
+  fEffA = new TH2F("feffa","Effective area",fEbins,logemin,logemax,fNTbins,fTbin);
+  //borra
+  TH1F eff("eff","Effective area",fEbins,logemin,logemax);
+  
   // fill the spectrum of the survivors
   TH2F* hpass= new TH2F("hpass","Survivors",esbins,logemin,logemax,fNTbins,fTbin);
@@ -289,16 +330,41 @@
   // compute the effective areas
   for(Int_t j=0;j<fNTbins;j++)
-    for(Int_t i=0;i<esbins;i++)
+    for(Int_t i=0;i<fEbins;i++)
       {
-	if(horig->GetBinContent(i+1,j+1)<=0)
+	Float_t effarea =0;
+	Float_t wtot = 0;
+	for(Int_t k=0;k<fEsubbins;k++)
+	  {
+	    Float_t numerator = hpass->GetBinContent(i*fEsubbins+k+1,j+1);
+	    Float_t denominator = horig->GetBinContent(i*fEsubbins+k+1,j+1);
+	    
+	    if(denominator<=0)
+	      cout << "MEffAreaAndCoeffCalc::ComputeEffectiveAreas Warning: energy subbin " << i*fEsubbins+k <<", theta bin " << j << " contains no events" << endl;
+	    else
+	      {
+		const Double_t ew = TMath::Power(10,logemin+(i*fEsubbins+k+0.5)*desub);
+		const Double_t ww = fSpec->Eval(ew);
+		effarea+=Atot*numerator/denominator*ww;
+		wtot   += ww;
+	      }
+	  }
+	if(!wtot)
 	  {
 	    cout << "MEffAreaAndCoeffCalc::ComputeEffectiveAreas Warning: energy bin " << i <<", theta bin " << j << " contains no events setting effective area to 0" << endl;
 	    fEffA->SetBinContent(i+1,j+1,0);
+	  }	    
+	else
+	  {
+	    fEffA->SetBinContent(i+1,j+1,effarea/wtot);	
+	    if(j==0)
+	      {
+		cout << "*****" << i << ": " << effarea/wtot << endl;
+		eff.SetBinContent(i+1,effarea/wtot);
+	      }
 	  }
-	else
-	  fEffA->SetBinContent(i+1,j+1,Atot*hpass->GetBinContent(i+1,j+1)/horig->GetBinContent(i+1,j+1));
       }
-  
-
+
+  // borra
+  eff.Write();
   
   delete hpass;
Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.h
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.h	(revision 6368)
+++ trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.h	(revision 6465)
@@ -18,11 +18,11 @@
  private:
 
-  TF1* fSpec;        // function used to parametrized the spectrum
-  TH1F* fHorig;      // histogram with the original energy spectrum
+  TF1* fSpec;        // function used to parametrize the spectrum
+  TH1F* fHorig;      // histogram with the original sample energy spectrum
 
-  Int_t fEbins;      // number of bins to build spectrum
-  Int_t fEsubbins;   // number of subbins per big bin
   Double_t fEmin;    // Minimum energy in GeV
   Double_t fEmax;    // Maximum energy in GeV
+  Int_t fEbins;      // number of bins to build spectrum
+  Int_t fEsubbins;   // number of subbins per big bin (to compute weights, eff areas...)
 
   Double_t* fWeight; // array containing weights
@@ -30,9 +30,11 @@
   TH2F* fEffA;       // histogram containing effective areas
 
-  TChain* fCini;     // chain for initial MC events (even those not triggering)
+  TChain* fCini;     // chain for initial MC files (before trigger)
   TChain* fCcut;     // chain for surviving MC events (after cuts)
 
   MHillas* fHillas;  // pointer to the MHillas Branch
   MMcEvt*  fMcEvt;   // pointer to the MMcEvt Branch
+
+  TFile* fFile; // output file (for debugging only)
 
  protected:
@@ -44,4 +46,5 @@
 
  public:
+
   MEffAreaAndCoeffCalc();
 
@@ -55,6 +58,5 @@
   void SetEmax(Float_t x)   {fEmax=x;}
 
-  void AddFileToInitialMC(const Char_t* name) {fCini->Add(name);}
-  void AddFileToFinalMC(const Char_t* name) {fCcut->Add(name);}
+  void AddFile(const Char_t* name) {fCini->Add(name); fCcut->Add(name);}
 
   TH2F* GetEffectiveAreaHisto() {return fEffA;}
