source: trunk/MagicSoft/Mars/mtemp/mifae/macros/computeCoeff.C@ 6359

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;
}