source: trunk/MagicSoft/Mars/macros/unfoldCoeff.C@ 7005

/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Javier Rico, 12/2005 <mailto:jrico@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2005
!
!
\* ======================================================================== */

// 
// Example macro on how to calculate coefficients for unfolding and
// effective areas.
//.The input is a file of the MC events surviving a certain kind 
// of analysis. It must contain the branch MEnergyEst with energy 
// estimator.
// It must also contain a tree called "OriginalMC" with a branch
// "MMcEvtBasic" and one entry per original simulated shower used to produce
// the events in that input file. 
//
// The macro performs 2 loops.
// In the first one the histograms of the original MC events are filled,
// needed by the coefficient (MMcUnfoldCoeffCalc) and effective area
// (MMcCollectionAreaCalc) computation
// In the second one we deal with the survivors and compute effective areas
// and coefficients

void unfoldCoeff(TString filename="/data/19990101_10002_Q_MCGamTestLZA_E_10_5.root", TString outname="coeff.root")
{

  ///////////////////////////////////////////////////////////////////
  // First loop: over all events processed by corsika
  ///////////////////////////////////////////////////////////////////
  MParList  parlist;
  MTaskList tasklist;
  parlist.AddToList(&tasklist);
  
  // theta binning
  Int_t zbins = 2;
  TArrayD edges(zbins+1);
  edges[0] = 0;  
  edges[1] = 10;
  edges[2] = 20;

  MBinning binstheta("binsTheta");
  binstheta.SetEdges(edges);               // Theta [deg]
  
  // energy binning
  const Int_t ebins=10;
  const Float_t emin=10;
  const Float_t emax=1000;
  
  MBinning binsenergy("binsEnergy");  
  binsenergy.SetEdgesLog(ebins,emin,emax); // Energy [GeV]
  
  parlist.AddToList(&binsenergy);
  parlist.AddToList(&binstheta);
  
  // Tentative spectrum
  TF1 spectrum("func","pow(x,-1.6)",2.,20000.);
  
  // Setup out tasks:
  MReadMarsFile reader("OriginalMC", filename);
  reader.DisableAutoScheme();
  
  // Create unfolding coefficients object and add it to the parameter list
  MHMcCollectionArea aeff;
  MHMcUnfoldCoeff coeff;
  parlist.AddToList(&coeff);
  parlist.AddToList(&aeff);
  
  // Task which fills the necessary histograms in MHMcCollectionArea
  MMcUnfoldCoeffCalc coeffcalc;
  MMcCollectionAreaCalc aeffcalc;
  
  coeffcalc.SetSpectrum(&spectrum);
  aeffcalc.SetSpectrum(&spectrum);
  
  tasklist.AddToList(&reader);
  tasklist.AddToList(&coeffcalc);
  tasklist.AddToList(&aeffcalc);
  
  // set up the loop for the processing
  MEvtLoop magic;
  magic.SetParList(&parlist);
  
  if (!magic.Eventloop())
    return;
  
  tasklist.PrintStatistics();
  
  ///////////////////////////////////////////////////////////////////
  // Second loop: now fill the histograms of the survivors
  ///////////////////////////////////////////////////////////////////
    
  MParList parlist2;
  MTaskList tasklist2;
  parlist2.AddToList(&tasklist2);
  parlist2.AddToList(&coeff);
  parlist2.AddToList(&aeff);
  
  MReadMarsFile reader2("Events", filename);
  reader2.DisableAutoScheme();
  
  tasklist2.AddToList(&reader2);
  tasklist2.AddToList(&aeffcalc);
  tasklist2.AddToList(&coeffcalc);
  
  //
  // set up the loop for the processing
  //
  MEvtLoop magic2;
  magic2.SetParList(&parlist2);
  
  if (!magic2.Eventloop())
    return;
  
  tasklist2.PrintStatistics();
  
  
  
  // Dummy cross-check, see if we recover the original spectrum
  const UInt_t ncutfiles=1;
  Char_t* cutName[ncutfiles]={"/data/19990101_10002_Q_MCGamTestLZA_E_10_5.root"};
  
  TChain* ccut = new TChain("Events");
  for(Int_t i = 0; i < ncutfiles; i++)
    ccut->Add(cutName[i]);
  //  ccut->SetAlias("logestenergy","log10(MHillas.fSize/0.18/15.)");
  ccut->SetAlias("logestenergy","log10(MEnergyEst.fEnergy)");
  ccut->SetAlias("theta","MMcEvt.fTelescopeTheta*180./3.14159");
  
  const Double_t logemin = TMath::Log10(emin);
  const Double_t logemax = TMath::Log10(emax);
  TH1D* hspec = new TH1D("hspec","Spectrum",ebins,logemin,logemax);
  hspec->Sumw2();
  ccut->Draw("logestenergy>>hspec","theta<10","goff");
  
  for(Int_t i=0;i<ebins;i++)
    {
      const Float_t uncorrval = hspec->GetBinContent(i+1);
      const Float_t effa   = ((TH2D*) aeff.GetHistCoarse())->GetBinContent(i+1,1);
      const Float_t unfold = ((TH2D*)coeff.GetHistCoeff())->GetBinContent(i+1,1);

      const Float_t euncorrval = hspec->GetBinError(i+1);
      const Float_t eeffa   = ((TH2D*) aeff.GetHistCoarse())->GetBinError(i+1,1);
      const Float_t eunfold = ((TH2D*)coeff.GetHistCoeff())->GetBinError(i+1,1);

      Float_t corrval,ecorrval;
      if(effa>0)
        {         
          corrval  = uncorrval*unfold/effa;
          if(uncorrval>0 && effa>0 & unfold>0)      
            ecorrval = corrval*TMath::Sqrt(euncorrval/uncorrval*euncorrval/uncorrval+
                                           eeffa/effa*eeffa/effa+
                                           eunfold/unfold*eunfold/unfold);       
          else
            ecorrval = 0;
        }
      else
        {
          corrval = 0;
          ecorrval = 0;
        }
      hspec->SetBinContent(i+1,corrval);
      hspec->SetBinError(i+1,ecorrval);
    }
  
  // Write histogram to a file in case an output
  // filename has been supplie
  if (outname.IsNull())
    return;
  
  TFile f(outname, "recreate");
  if (!f)
    return;
  
  coeff.GetHistAll()->Write();
  coeff.GetHistWeight()->Write();
  coeff.GetHistMcE()->Write();
  coeff.GetHistEstE()->Write();
  coeff.GetHistCoeff()->Write();
  aeff.GetHistCoarse()->Write();
  aeff.Write();
  hspec->Write();
}