Index: trunk/MagicSoft/Mars/mtemp/mifae/macros/computeCoeff.C
===================================================================
--- trunk/MagicSoft/Mars/mtemp/mifae/macros/computeCoeff.C	(revision 6359)
+++ trunk/MagicSoft/Mars/mtemp/mifae/macros/computeCoeff.C	(revision 6359)
@@ -0,0 +1,93 @@
+
+
+void computeCoeff()
+{    
+  MEffAreaAndCoeffCalc calc;
+
+  // initial MC files
+  calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_0_7_1000to1009_w0.root");
+  calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_0_7_1130to1139_w0.root");
+  calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_0_7_1520to1529_w0.root");
+  calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_90_7_1260to1269_w0.root");
+  calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_90_7_1390to1399_w0.root");
+  calc.AddFileToInitialMC("/data/camera_NCL/Gamma_zbin0_90_7_1650to1659_w0.root");
+  
+
+  // files with remaining events after cuts (or where to apply cuts) see warning above
+  calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_0_7_1000to1009_w0_cleaned4035.root");
+  calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_0_7_1130to1139_w0_cleaned4035.root");
+  calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_0_7_1520to1529_w0_cleaned4035.root");
+  calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_90_7_1260to1269_w0_cleaned4035.root");
+  calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_90_7_1390to1399_w0_cleaned4035.root");
+  calc.AddFileToFinalMC("/data/camera_CL/Gamma_zbin0_90_7_1650to1659_w0_cleaned4035.root");
+
+
+  // define the funtion of the desired spectrum
+  calc.SetFunction("4.e9*pow(x,-2.6+1)",10.,10000.);
+  calc.ComputeAllFactors();
+  
+  /************************************************/
+  /* Build spectrum                               */
+  /* now it's just a cross check with the same MC */
+  /************************************************/
+  const Int_t ebins = 10;              // number of bins to build spectrum
+  const Int_t esubbins = 20;           // number of subbins per big bin
+  const Double_t emin=10;
+  const Double_t emax=10000;
+  const Double_t logemin = TMath::Log10(emin);
+  const Double_t logemax = TMath::Log10(emax);
+  const Double_t de = (logemax-logemin)/ebins; // bin size (in log)
+  const Double_t desub = de/esubbins; // subbin size (in log)
+  const Int_t esbins = ebins*esubbins; // total number of subbins
+
+  // remaining events after cuts
+  const UInt_t ncutfiles=6;
+  Char_t* cutName[ncutfiles]={
+    "/data/camera_CL/Gamma_zbin0_0_7_1000to1009_w0_cleaned4035.root",
+    "/data/camera_CL/Gamma_zbin0_0_7_1130to1139_w0_cleaned4035.root",
+    "/data/camera_CL/Gamma_zbin0_0_7_1520to1529_w0_cleaned4035.root",
+    "/data/camera_CL/Gamma_zbin0_90_7_1260to1269_w0_cleaned4035.root",
+    "/data/camera_CL/Gamma_zbin0_90_7_1390to1399_w0_cleaned4035.root",
+    "/data/camera_CL/Gamma_zbin0_90_7_1650to1659_w0_cleaned4035.root"
+  };
+
+  TChain* ccut = new TChain("Events");
+  for(Int_t i = 0; i < ncutfiles; i++)
+    ccut->Add(cutName[i]);
+  ccut->SetAlias("logestenergy","log10(MHillas.fSize/15.)");
+  const Int_t nentries = Int_t(ccut->GetEntries());
+
+  MHillas*  hillas;
+  MMcEvt*   mcevt;
+  ccut->SetBranchStatus("*",0);
+  ccut->SetBranchStatus("MHillas.*",1);
+  ccut->SetBranchStatus("MMcEvt.*",1);
+  ccut->SetBranchAddress("MHillas.",&hillas);
+  ccut->SetBranchAddress("MMcEvt.",&mcevt);
+
+  TH1D* hspec = new TH1D("hspec","Spectrum",ebins,logemin,logemax);
+  hspec->Sumw2();
+
+  for(Int_t i=0;i<nentries;i++)
+    {
+      ccut->GetEntry(i);  
+      
+      // OJO!! Estimated energy will be taken directly from some input container
+      Float_t estim  = hillas->GetSize()/15.;
+
+      UInt_t effabin = UInt_t((TMath::Log10(estim)-logemin)/desub);
+      UInt_t coefbin = UInt_t((TMath::Log10(estim)-logemin)/de);
+      
+      Float_t effa  =  calc.GetEffectiveArea(effabin);
+      Float_t unfold = calc.GetCoefficient(coefbin);
+      
+      if(effa)
+	hspec->Fill(TMath::Log10(estim),unfold/effa*1e9);
+    }      
+
+  // SAVE RESULTS
+  TFile file("prueba.root","RECREATE");
+  hspec->Write();
+
+  return;
+}