Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc	(revision 6465)
+++ trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.cc	(revision 6486)
@@ -47,7 +47,11 @@
 using namespace std;
 
-const Int_t fNTbins = 2;
-const Double_t fTbin[fNTbins+1] = {0,10,20};
-
+const Int_t ntbins=1;                      // default number of theta bins
+const Double_t tbin[ntbins+1] = {0,90};    // default theta bins bounds
+const Int_t nebins = 10;                   // default number of energy bins
+const Float_t emin = 10.;                  // default minimum energy value
+const Float_t emax = 10000.;               // default maximum energy value
+const Int_t nsubbins = 20;                 // default number of subbins per bin
+const Char_t* deff = "4.e9*pow(x,-2.6+1)"; // default spectrum function
 
 // -------------------------------------------------------------------------
@@ -56,12 +60,11 @@
 //
 MEffAreaAndCoeffCalc::MEffAreaAndCoeffCalc()
-  : fSpec(NULL), fHorig(NULL), fEmin(10.), fEmax(10000.), fEbins(10), fEsubbins(20),
-    fCoeff(NULL), fEffA(NULL), fFile(NULL)
+  : fSpec(NULL),  fEmin(emin), fEmax(emax), fEbins(nebins), fEsubbins(nsubbins), 
+    fWeight(NULL), fCoeff(NULL), fEffA(NULL)//, fFile(NULL)
 {
   // set the default function
-  SetFunction("4.e9*pow(x,-2.6+1)");  
+  SetFunction(deff);  
 
   // create the TChains 
-  // OJO!! make sure they read the appropriate branch
   fCini = new TChain("OriginalMC");
   fCcut = new TChain("Events");
@@ -76,7 +79,10 @@
   fCcut->SetBranchAddress("MHillas.",&fHillas);
   fCcut->SetBranchAddress("MMcEvt.",&fMcEvt);
+
+  // initial value of the zenith angle binning
+  SetThetaBinning(ntbins,tbin);
   
   // borra
-  fFile = new TFile("test.root","RECREATE");
+  //  fFile = new TFile("coeftest.root","RECREATE");
 }
 
@@ -90,4 +96,7 @@
     delete fSpec;
 
+  if(fTbin)
+    delete [] fTbin;
+
   if(fWeight)
     delete [] fWeight;
@@ -99,9 +108,20 @@
     delete fEffA;
 
-  if(fHorig) 
-    delete fHorig;
 
   delete fCini;
   delete fCcut;
+}
+
+// -------------------------------------------------------------------------
+//
+// Set the binning in zenith angle
+//
+void MEffAreaAndCoeffCalc::SetThetaBinning(Int_t n, const Double_t* binlist)
+{
+  fNTbins=n;
+  if(fTbin) delete [] fTbin;
+  fTbin = new Double_t[n+1];
+  for(Int_t i=0;i<n+1;i++)
+    fTbin[i] = binlist[i];
 }
 
@@ -144,13 +164,4 @@
 void MEffAreaAndCoeffCalc::FillOriginalSpectrum()
 {  
-  const Double_t logemin = TMath::Log10(fEmin);
-  const Double_t logemax = TMath::Log10(fEmax);
-  const Int_t esbins = fEbins*fEsubbins; // total number of subbins
-  
-  if(fHorig) delete fHorig;
-  fHorig   = new TH1F("fhorig","Original energy spectrum",esbins,logemin,logemax);
-  fCini->Draw("logenergy>>fhorig","","goff");
-  // borra
-  fHorig->Write();
 }
 
@@ -161,37 +172,54 @@
 void MEffAreaAndCoeffCalc::ComputeWeights()
 {  
-  const Double_t logemin = TMath::Log10(fEmin);
-  const Double_t logemax = TMath::Log10(fEmax);
-  const Double_t de = (logemax-logemin)/fEbins; // bin size (in log)
-  const Double_t desub = de/fEsubbins; // subbin size (in log)
-  const Int_t esbins = fEbins*fEsubbins; // total number of subbins
-
+  // OJO!! maybe this should be hard-coded somewhere else
+  const Float_t abslogmin=0.;
+  const Float_t abslogmax=5.;
+
+  const Float_t logemin = TMath::Log10(fEmin);
+  const Float_t logemax = TMath::Log10(fEmax);
+  const Float_t de    = (logemax-logemin)/fEbins; // bin size (in log)
+  const Float_t desub = de/fEsubbins; // subbin size (in log)
+  const Int_t esbins   = fEbins*fEsubbins; // total number of subbins
+
+  // compute the binning for weights histogram
+  const Int_t nmindist = (logemin>abslogmin)? Int_t((logemin-abslogmin)/desub) : 0;
+  const Int_t nmaxdist = (logemax<abslogmax)? Int_t((abslogmax-logemax)/desub) : 0;
+
+  fLogEWmin = logemin-desub*nmindist;
+  fLogEWmax = logemax+desub*nmaxdist;
+  fEWbins   = nmindist+esbins+nmaxdist;
+  
   // reset the weights array
   if(fWeight)
     delete [] fWeight;
-  fWeight = new Double_t[esbins];
+  fWeight = new Double_t[fEWbins];
     
-  if(!fHorig)
-    FillOriginalSpectrum();
-  
-  // borra
-  TH1F hw("hw","hw",esbins,fEmin,fEmax);
-  for(Int_t i=0;i<esbins;i++)
-    {
-      const Double_t denom = fHorig->GetBinContent(i+1);              // number of events in initial spectrum
-      const Double_t ew    = TMath::Power(10,logemin+(i+0.5)*desub); // real energy
-      const Double_t numer = fSpec->Eval(ew);                        // number of events for the required spectrum
-      if(denom)
+
+  TH1F* horigs = new TH1F("horigs","Original energy spectrum",fEWbins,fLogEWmin,fLogEWmax);
+  fCini->Draw("logenergy>>horigs","","goff");
+  // borra
+  //  horigs->Write();
+
+  // borra
+  //  TH1F hw("hw","hw",fEWbins,fLogEWmin,fLogEWmax);
+  for(Int_t i=0;i<fEWbins;i++)
+    {
+      const Float_t denom = horigs->GetBinContent(i+1);               // number of events in initial spectrum
+      const Float_t ew    = TMath::Power(10,fLogEWmin+(i+0.5)*desub); // real energy
+      const Float_t numer = fSpec->Eval(ew);                          // number of events for the required spectrum
+      if(denom>10)
 	fWeight[i]=numer/denom;
       else
 	{
-	  cout << "MEffAreaAndCoeffCalc::ComputeWeights Warning: no statistic to compute weight for energy " << ew << ", setting it to -1 "  << endl;
+	  //	  cout << "MEffAreaAndCoeffCalc::ComputeWeights Warning: no statistic to compute weight for energy " << ew << ", setting it to -1 "  << endl;
 	  fWeight[i]=-1;
 	}
       // borra
-      hw.SetBinContent(i+1,fWeight[i]);
-    }
-  // borra
-  hw.Write();
+      //      hw.SetBinContent(i+1,fWeight[i]);
+    }
+  // borra
+  //  hw.Write();
+
+  delete horigs;
 }
 
@@ -208,9 +236,8 @@
     }
 
-  const Double_t logemin = TMath::Log10(fEmin);
-  const Double_t logemax = TMath::Log10(fEmax);
-  const Double_t de = (logemax-logemin)/fEbins; // bin size (in log)
-  const Double_t desub = de/fEsubbins; // subbin size (in log)
-  const Int_t esbins = fEbins*fEsubbins; // total number of subbins
+  const Float_t logemin = TMath::Log10(fEmin);
+  const Float_t logemax = TMath::Log10(fEmax);
+  const Float_t de = (logemax-logemin)/fEbins; // bin size (in log)
+  const Float_t desub = de/fEsubbins; // subbin size (in log)
   const Int_t nentries = Int_t(fCcut->GetEntries());
  
@@ -220,7 +247,7 @@
 
   // 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);
+  //  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
@@ -238,23 +265,33 @@
       Float_t emc   = fMcEvt->GetEnergy();
       Float_t estim = fHillas->GetSize()/0.18/15.;
+
+      if((emc<fEmin && estim<fEmin) || (emc>fEmax && estim>fEmax) ||
+	 (emc<fEmin && estim>fEmax) || (emc>fEmax && estim<fEmin))
+	continue;
+
       Float_t theta = fMcEvt->GetTheta()*180./3.14159; // zenith angle (in degrees)
 
       // compute the bin where the weight is taken from
-      Int_t wbin = Int_t((TMath::Log10(emc)-logemin)/desub);
+      Int_t wbin = Int_t((TMath::Log10(emc)-fLogEWmin)/desub);
 
       // fill the histograms
-      if(wbin<esbins)
+      if(wbin<fEWbins)
 	{
-	  hest->Fill(TMath::Log10(estim),theta,fWeight[wbin]); 
-	  hmc->Fill(TMath::Log10(emc),theta,fWeight[wbin]);
-	  // borra
-	  if(theta<fTbin[1])
+	  if(fWeight[wbin]>0)
 	    {
-	      hest1->Fill(TMath::Log10(estim),fWeight[wbin]); 
-	      hmc1->Fill(TMath::Log10(emc),fWeight[wbin]);
+	      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 no computed weight (lack of MC statistics), skipping" << endl;
 	}
       else
-	cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has energy out of bounds, skipping" << endl;
+      	cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has energy out of bounds, skipping" << endl;
     }
 
@@ -266,12 +303,12 @@
 	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);
+// 	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);
 	    //borra
-	    if(j==0 && den1)
-	      coef1->SetBinContent(i+1,num1/den1);
+// 	    if(j==0 && den1)
+// 	      coef1->SetBinContent(i+1,num1/den1);
 	  }
 	else
@@ -284,7 +321,7 @@
 
   //borra
-  hmc1->Write();
-  hest1->Write();
-  coef1->Write();
+//   hmc1->Write();
+//   hest1->Write();
+//   coef1->Write();
 
   delete hmc;
@@ -305,12 +342,11 @@
 
   //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 = 1.; //total area (in cm^2)
-  const Double_t logemin = TMath::Log10(fEmin);
-  const Double_t logemax = TMath::Log10(fEmax);
+  const Float_t radius = 30000.; // generation radius (in cm)
+  const Float_t Atot = 3.14159*radius*radius; //total area (in cm^2)
+  const Float_t logemin = TMath::Log10(fEmin);
+  const Float_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)
+  const Float_t de = (logemax-logemin)/fEbins; // bin size (in log)
+  const Float_t desub = de/fEsubbins; // subbin size (in log)
 
   // reset the effective areas
@@ -319,5 +355,5 @@
   fEffA = new TH2F("feffa","Effective area",fEbins,logemin,logemax,fNTbins,fTbin);
   //borra
-  TH1F eff("eff","Effective area",fEbins,logemin,logemax);
+//   TH1F eff("eff","Effective area",fEbins,logemin,logemax);
   
   // fill the spectrum of the survivors
@@ -343,6 +379,6 @@
 	    else
 	      {
-		const Double_t ew = TMath::Power(10,logemin+(i*fEsubbins+k+0.5)*desub);
-		const Double_t ww = fSpec->Eval(ew);
+		const Float_t ew = TMath::Power(10,logemin+(i*fEsubbins+k+0.5)*desub);
+		const Float_t ww = fSpec->Eval(ew);
 		effarea+=Atot*numerator/denominator*ww;
 		wtot   += ww;
@@ -357,14 +393,15 @@
 	  {
 	    fEffA->SetBinContent(i+1,j+1,effarea/wtot);	
-	    if(j==0)
-	      {
-		cout << "*****" << i << ": " << effarea/wtot << endl;
-		eff.SetBinContent(i+1,effarea/wtot);
-	      }
+	    // borra
+// 	    if(j==0)
+// 	      {
+// 		//		cout << "*****" << i << ": " << effarea/wtot << endl;
+// 		eff.SetBinContent(i+1,effarea/wtot);
+// 	      }
 	  }
       }
 
   // borra
-  eff.Write();
+//   eff.Write();
   
   delete hpass;
Index: trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.h
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.h	(revision 6465)
+++ trunk/MagicSoft/Mars/mtemp/mifae/library/MEffAreaAndCoeffCalc.h	(revision 6486)
@@ -19,10 +19,16 @@
 
   TF1* fSpec;        // function used to parametrize the spectrum
-  TH1F* fHorig;      // histogram with the original sample energy spectrum
 
-  Double_t fEmin;    // Minimum energy in GeV
-  Double_t fEmax;    // Maximum energy in GeV
+  Float_t fEmin;    // Minimum energy in GeV
+  Float_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...)
+
+  Float_t fLogEWmin; // Log Minimum energy for weights (in GeV)
+  Float_t fLogEWmax; // Log Maximum energy for weights (in GeV)
+  Int_t    fEWbins;   // Number of bins for weights
+
+  Int_t fNTbins;      // Number of bins in zenith angle
+  Double_t* fTbin;     // array containing bin boundaries (size must be fNTbins+)
 
   Double_t* fWeight; // array containing weights
@@ -58,4 +64,6 @@
   void SetEmax(Float_t x)   {fEmax=x;}
 
+  void SetThetaBinning(Int_t n, const Double_t* binlist);
+
   void AddFile(const Char_t* name) {fCini->Add(name); fCcut->Add(name);}
 
