[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / analysis_macro.C

#include "AliAODTrack.h"
#include <iostream>
#include "TFile.h"
#include "TH1I.h"
#include "TH2.h"
#include "AliSpectraAODHistoManager.h"
#include "AliSpectraAODEventCuts.h"
#include "AliSpectraAODTrackCuts.h"
#include "TCanvas.h"
#include "TLegend.h"

using namespace std;

void analysis_macro(const char * dataFile = "Pt.AOD.1._data_ptcut.root", const char * mcFile ="Pt.AOD.1._MC.root") {
   
  // gSystem->Load("libTree.so");
  // gSystem->Load("libGeom.so");
  // gSystem->Load("libVMC.so");
  // gSystem->Load("libPhysics.so");
  // gSystem->Load("libSTEERBase.so");
  // gSystem->Load("libESD.so");
  // gSystem->Load("libAOD.so");
  // gSystem->Load("libANALYSIS.so");
  // gSystem->Load("libANALYSISalice.so");
  // gSystem->Load("libANALYSIS");
  // gSystem->Load("libANALYSISalice");
  // gROOT->LoadMacro("$ALICE_ROOT/PWG0/CreateESDChain.C");
  // gSystem->AddIncludePath("-I$ALICE_ROOT/include");
  // gStyle->SetPalette(1);
  // gStyle->SetFillColor(kWhite);

  // gROOT->LoadMacro("AliSpectraAODTrackCuts.cxx+g");
  // gROOT->LoadMacro("AliSpectraAODEventCuts.cxx+g");
  // gROOT->LoadMacro("AliSpectraAODHistoManager.cxx+g");
  // gROOT->LoadMacro("AliAnalysisTaskSpectraAOD.cxx+g");
 
  const char * histoname[] =
    {
      // names of histos we want to draw. please observe order of histos, elsewise name will not make sense
      "#frac{MC_{#sigma, rec}}{Data}, P^{+}",
      "#frac{MC_{#sigma, rec}}{Data}, K^{+}",
      "#frac{MC_{#sigma, rec}}{Data}, #pi^{+}",
      "#frac{MC_{#sigma, rec}}{Data}, P^{-}",
      "#frac{MC_{#sigma, rec}}{Data}, K^{-}",
      "#frac{MC_{#sigma, rec}}{Data}, #pi^{-}",      
          
      "#frac{MC_{#sigma, rec, prim}}{Data}, P^{+}",
      "#frac{MC_{#sigma, rec, prim}}{Data}, K^{+}",
      "#frac{MC_{#sigma, rec, prim}}{Data}, #pi^{+}",
      "#frac{MC_{#sigma, rec, prim}}{Data}, P^{-}",
      "#frac{MC_{#sigma, rec, prim}}{Data}, K^{-}",
      "#frac{MC_{#sigma, rec, prim}}{Data}, #pi^{-}",     
         
      "#frac{MC_{#sigma, rec, sec}}{Data}, P^{+}",
      "#frac{MC_{#sigma, rec, sec}}{Data}, K^{+}",
      "#frac{MC_{#sigma, rec, sec}}{Data}, #pi^{+}",
      "#frac{MC_{#sigma, rec, sec}}{Data}, P^{-}",
      "#frac{MC_{#sigma, rec, sec}}{Data}, K^{-}",
      "#frac{MC_{#sigma, rec, sec}}{Data}, #pi^{-}",     
         
      "#frac{MC_{true, gen}}{MC_{#sigma, rec}}, P^{+}",
      "#frac{MC_{true, gen}}{MC_{#sigma, rec}}, K^{+}",
      "#frac{MC_{true, gen}}{MC_{#sigma, rec}}, #pi^{+}",
      "#frac{MC_{true, gen}}{MC_{#sigma, rec}}, P^{-}",
      "#frac{MC_{true, gen}}{MC_{#sigma, rec}}, K^{-}",
      "#frac{MC_{true, gen}}{MC_{#sigma, rec}}, #pi^{-}",     
              
      "#frac{MC_{true, gen}}{MC_{#sigma, rec}} x Data, P^{+}",
      "#frac{MC_{true, gen}}{MC_{#sigma, rec}} x Data, K^{+}",
      "#frac{MC_{true, gen}}{MC_{#sigma, rec}} x Data, #pi^{+}",
      "#frac{MC_{true, gen}}{MC_{#sigma, rec}} x Data, P^{-}",
      "#frac{MC_{true, gen}}{MC_{#sigma, rec}} x Data, K^{-}",
      "#frac{MC_{true, gen}}{MC_{#sigma, rec}} x Data, #pi^{-}",     

      "Data, P^{+}",
      "Data, K^{+}",
      "Data, #pi^{+}",
      "Data, P^{-}",
      "Data, K^{-}",
      "Data, #pi^{-}",      

      "#frac{MC_{rec #sigma}}{MC_{true, rec}}, P^{+}",
      "#frac{MC_{rec #sigma}}{MC_{true, rec}}, K^{+}",
      "#frac{MC_{rec #sigma}}{MC_{true, rec}}, #pi^{+}",
      "#frac{MC_{rec #sigma}}{MC_{true, rec}}, P^{-}",
      "#frac{MC_{rec #sigma}}{MC_{true, rec}}, K^{-}",
      "#frac{MC_{rec #sigma}}{MC_{true, rec}}, #pi^{-}",     

      "#frac{MC_{rec #sigma sec}}{MC_{#sigma, rec}}, P^{+}",
      "#frac{MC_{rec #sigma sec}}{MC_{#sigma, rec}}, K^{+}",
      "#frac{MC_{rec #sigma sec}}{MC_{#sigma, rec}}, #pi^{+}",
      "#frac{MC_{rec #sigma sec}}{MC_{#sigma, rec}}, P^{-}",
      "#frac{MC_{rec #sigma sec}}{MC_{#sigma, rec}}, K^{-}",
      "#frac{MC_{rec #sigma sec}}{MC_{#sigma, rec}}, #pi^{-}",     
    };

  const char * histotitle[] =
    {
      // titles of histos we want to draw. please observe order of histos, elsewise name will not make sense
      "MC #sigma, rec / Data, P^{+}",
      "MC #sigma, rec / Data, K^{+}",
      "MC #sigma, rec / Data, #pi^{+}",
      "MC #sigma, rec / Data, P^{-}",
      "MC #sigma, rec / Data, K^{-}",
      "MC #sigma, rec / Data, #pi^{-}",     
   
      "MC #sigma, rec, prim / Data, P^{+}",
      "MC #sigma, rec, prim  / Data, K^{+}",
      "MC #sigma, rec, prim  / Data, #pi^{+}",
      "MC #sigma, rec, prim  / Data, P^{-}",
      "MC #sigma, rec, prim  / Data, K^{-}",
      "MC #sigma, rec, prim  / Data, #pi^{-}",      
    
      "MC #sigma, rec, sec / Data, P^{+}",
      "MC #sigma, rec, sec / Data, K^{+}",
      "MC #sigma, rec, sec / Data, #pi^{+}",
      "MC #sigma, rec, sec / Data, P^{-}",
      "MC #sigma, rec, sec / Data, K^{-}",
      "MC #sigma, rec, sec / Data, #pi^{-}",      
          
      "Correction factor, P^{+}",
      "Correction factor, K^{+}",
      "Correction factor, #pi^{+}",
      "Correction factor, P^{-}",
      "Correction factor, K^{-}",
      "Correction factor, #pi^{-}",      

      "True Data, P^{+}",
      "True Data, K^{+}",
      "True Data, #pi^{+}",
      "True Data, P^{-}",
      "True Data, K^{-}",
      "True Data, #pi^{-}",      

      "Raw Data, P^{+}",
      "Raw Data, K^{+}",
      "Raw Data, #pi^{+}",
      "Raw Data, P^{-}",
      "Raw Data, K^{-}",
      "Raw Data, #pi^{-}",      

      "Identificaiton quality, P^{+}",
      "Identificaiton quality, K^{+}",
      "Identificaiton quality, #pi^{+}",
      "Identificaiton quality, P^{-}",
      "Identificaiton quality, K^{-}",
      "Identificaiton quality, #pi^{-}",     

      "Contribution of secondaries, P^{+}",
      "Contribution of secondaries, K^{+}",
      "Contribution of secondaries, #pi^{+}",
      "Contribution of secondaries, P^{-}",
      "Contribution of secondaries, K^{-}",
      "Contribution of secondaries, #pi^{-}",     
    };

  // Open root MC file and get classes
  cout << "Analysis Macro" << endl;
  cout << "  > Reading MC data" << endl;
  TFile *_mc = TFile::Open(mcFile);
  AliSpectraAODHistoManager* hman_mc = (AliSpectraAODHistoManager*) _mc->Get("SpectraHistos");
  AliSpectraAODEventCuts* ecuts_mc = (AliSpectraAODEventCuts*) _mc->Get("Event Cuts");
  AliSpectraAODTrackCuts* tcuts_mc = (AliSpectraAODTrackCuts*) _mc->Get("Track Cuts");
  // print info about mc track and Event cuts
  cout << " -- Info about MC -- "<< endl;
  ecuts_mc->PrintCuts();
  tcuts_mc->PrintCuts();
  // get the mc histo's necessary to write the contamination histos from the rootfile
  TH1F* spectrahistos_mc[30];
  cout << " -- Reading and normalizing histograms -- " << endl;
  for( Int_t i = 0 ; i <= AliSpectraNameSpace::kNPtSpecies ; i++ ) { 
    spectrahistos_mc[i] = dynamic_cast<TH1F*>(hman_mc->GetHistogram((AliSpectraNameSpace::AODPtHist_t)i)); 
    Double_t events_mc =  1. / ecuts_mc->NumberOfEvents();
    if(!spectrahistos_mc[i]) {
      continue;
    }    
    spectrahistos_mc[i]->Scale(events_mc, "width");
  }
  cout << "  > Reading real data " << endl;
  // proceed likewise for data
  TFile *_data = TFile::Open(dataFile);
  AliSpectraAODHistoManager* hman_data = (AliSpectraAODHistoManager*) _data->Get("SpectraHistos");
  AliSpectraAODEventCuts* ecuts_data = (AliSpectraAODEventCuts*) _data->Get("Event Cuts");
  AliSpectraAODTrackCuts* tcuts_data = (AliSpectraAODTrackCuts*) _data->Get("Track Cuts");
  // print info about track and Event cuts
  cout << " -- Info about data -- " << endl;
  ecuts_data->PrintCuts();
  tcuts_data->PrintCuts();
  // get the histo's necessary to write the contamination histos from the rootfile
  // MF 22/02/2012
  // Get the raw data distribution:
  // 1. Get DCA vs pt histo for each particle
  // 2. Loop over all pt bins, projecting on DCA
  // 3. Fit the DCA data distribution with MC templates for primaries, sec. w., sec. m. (TFractionFitter)
  // 4. Extract from the fit raw primaries vs pt in data
  // 5. Correct for the primaries efficiency estimated in MC
  TH1F* spectrahistos_data[30];
  cout << " -- Reading and normalizing histograms -- " << endl;
  for( Int_t i = 0 ; i <= AliSpectraNameSpace::kNPtSpecies ; i++ ) { 
    spectrahistos_data[i] = dynamic_cast<TH1F*>(hman_data->GetHistogram((AliSpectraNameSpace::AODPtHist_t)i));
    if(!spectrahistos_data[i]) {
      continue;
    }    
    Double_t events_data =  1. / ecuts_data->NumberOfEvents();
    //normalize the histos    
    spectrahistos_data[i]->Scale(events_data, "width");
  }
  //create output file and select it as the active file
  TFile * output = new TFile("analysis_output.root", "RECREATE");
  output->cd();
  // Write the scaled data and MC histos in the output file for convenience
  output->mkdir("Data");
  output->cd("Data");
  for(Int_t ihist = 0; ihist <= AliSpectraNameSpace::kNPtSpecies; ihist++){
      spectrahistos_data[ihist]->Write();
  }
  output->mkdir("MC");
  output->cd("MC");
  for(Int_t ihist = 0; ihist <= AliSpectraNameSpace::kNPtSpecies; ihist++){
      spectrahistos_mc[ihist]->Write();
  }
  output->mkdir("ANALYSIS");
  output->cd("ANALYSIS");
  // perform some analysis!
  //  next loop: MC_Sigma_Rec over Rec_Data
  TCanvas* c1 = new TCanvas("MC reconstructed over Data","MC reconstructed over Data");
  c1->Divide(3,2);
  TH1F* kRatioPtMCSigmaRecSigmaRec[6];
  for (Int_t i = 0 ; i < 6 ; i ++ )
    {
      kRatioPtMCSigmaRecSigmaRec[i] = (TH1F*)(spectrahistos_mc[i])->Clone(histotitle[i]);
      kRatioPtMCSigmaRecSigmaRec[i]->Divide(spectrahistos_data[i]);
      c1->cd(1+i);
      kRatioPtMCSigmaRecSigmaRec[i]->SetTitle(histoname[i]);
      kRatioPtMCSigmaRecSigmaRec[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
      kRatioPtMCSigmaRecSigmaRec[i]->Draw();
      kRatioPtMCSigmaRecSigmaRec[i]->Write();
    }
  // next loop: MC_sigma_primaries over Rec_data
  TCanvas* c2 = new TCanvas("MC reconstructed primaries over Data", "MC reconstructed primaries over Data");
  c2->Divide(3,2);
  TH1F* kRatioPtMCSigmaRecPrimarySigmaRec[6];
  for (Int_t i = 0 ; i < 6 ; i ++ )
    {
      kRatioPtMCSigmaRecPrimarySigmaRec[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaPrimaryProtonPlus])->Clone(histotitle[i+6]);
      kRatioPtMCSigmaRecPrimarySigmaRec[i]->Divide(spectrahistos_data[i]);
      c2->cd(1+i);
      kRatioPtMCSigmaRecPrimarySigmaRec[i]->SetTitle(histoname[i+6]);
      kRatioPtMCSigmaRecPrimarySigmaRec[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
      kRatioPtMCSigmaRecPrimarySigmaRec[i]->Draw();
      kRatioPtMCSigmaRecPrimarySigmaRec[i]->Write();
    }
  // next loop: MC_Sigma_secondaries over Rec_data
  TCanvas* c3 = new TCanvas("MC reconstructed secondaries over Data","MC reconstructed secondaries over Data");
  c3->Divide(3,2);
  TH1F* kRatioPtMCSigmaRecSigmaRecSecondary[6];
  for (Int_t i = 0 ; i < 6 ; i ++ )
    {
      kRatioPtMCSigmaRecSigmaRecSecondary[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaSecondaryProtonPlus])->Clone(histotitle[i+12]);
      kRatioPtMCSigmaRecSigmaRecSecondary[i]->Divide(spectrahistos_data[i]);
      c3->cd(1+i);
      kRatioPtMCSigmaRecSigmaRecSecondary[i]->SetTitle(histoname[12+i]);
      kRatioPtMCSigmaRecSigmaRecSecondary[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
      kRatioPtMCSigmaRecSigmaRecSecondary[i]->Draw();
      kRatioPtMCSigmaRecSigmaRecSecondary[i]->Write();
    }
  //
  TCanvas* c4 = new TCanvas("MC true generated over MC true reconstructed","MC true generated over MC true reconstructed");
  c4->Divide(3,2);
  TH1F* kRatioPtMCTrueMCRec[6];
  for (Int_t i = 0 ; i < 6 ; i ++ )
    {
      kRatioPtMCTrueMCRec[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus])->Clone(histotitle[i+18]);
      kRatioPtMCTrueMCRec[i]->Divide(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaPrimaryProtonPlus]);
      c4->cd(1+i);
      kRatioPtMCTrueMCRec[i]->SetTitle(histoname[18+i]);
      kRatioPtMCTrueMCRec[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
      kRatioPtMCTrueMCRec[i]->Draw();
      kRatioPtMCTrueMCRec[i]->Write();
    }
  
  TCanvas* c5 = new TCanvas("True Data","True Data");
  c5->Divide(3,2);
  TH1F* kRatioPtTrueData[6];
  for (Int_t i = 0 ; i < 6 ; i ++ )
    {
      kRatioPtTrueData[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus])->Clone(histotitle[i+24]);
      kRatioPtTrueData[i]->Divide(spectrahistos_mc[i]); // Michele: Since we are correcting all contributions at once, this should be exaclty the same as in the data!
      kRatioPtTrueData[i]->Multiply(spectrahistos_data[i]);
      c5->cd(1+i);
      kRatioPtTrueData[i]->SetTitle(histoname[24+i]);
      kRatioPtTrueData[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
      kRatioPtTrueData[i]->Draw();
      kRatioPtTrueData[i]->Write();
    }
 
  TCanvas* c6 = new TCanvas("Raw data", "Raw data");
  c6->Divide(3,2);
  TH1F* kRawData[6];
  for (Int_t i = 0 ; i < 6 ; i ++ )
    {
      kRawData[i] = (TH1F*)(spectrahistos_data[i+AliSpectraNameSpace::kHistPtRecSigmaProtonPlus])->Clone(histotitle[i+30]);
      c6->cd(1+i);
      kRawData[i]->SetTitle(histoname[30+i]);
      kRawData[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
      kRawData[i]->Draw();
      kRawData[i]->Write();
    }
  
  TCanvas* c7 = new TCanvas("contamination", "contamination");
  c7->Divide(3,2);
  TH1F* kContamination[6];
  for (Int_t i = 0 ; i < 6 ; i ++ )
    {
      kContamination[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaProtonPlus])->Clone(histotitle[i+36]);
      kContamination[i]->Divide(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecTrueProtonPlus]);
      c7->cd(1+i);
      kContamination[i]->SetTitle(histoname[36+i]);
      kContamination[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
      kContamination[i]->Draw();
      kContamination[i]->Write();
    }
  
  TCanvas* c8 = new TCanvas("contamination secondary", "contamination secondary");
  c8->Divide(3,2);
  TH1F* kSecContamination[6];
  for (Int_t i = 0 ; i < 6 ; i ++ )
    {
      kSecContamination[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaSecondaryProtonPlus])->Clone(histotitle[i+42]);
      kSecContamination[i]->Divide(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaProtonPlus]);
      c8->cd(1+i);
      kSecContamination[i]->SetTitle(histoname[42+i]);
      kSecContamination[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
      kSecContamination[i]->Draw();
      kSecContamination[i]->Write();
    }
   
  TCanvas* c9 = new TCanvas("true data and generated", "true data and generated");
  c9->Divide(3,2);
  TH1F* ktrueandgenerated[6];
  for (Int_t i = 0 ; i < 6 ; i ++ )
    {
      // kSecContamination[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaSecondaryProtonPlus])->Clone(histotitle[i+42]);
      // kSecContamination[i]->Divide(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaProtonPlus]);
      c9->cd(1+i);
      // kSecContamination[i]->SetTitle(histoname[42+i]);
      // kSecContamination[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
      kRatioPtTrueData[i]->Draw();
      spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus]->SetMarkerColor(kRed);
      spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus]->SetMarkerStyle(kOpenCircle);
      spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus]->SetLineColor(kRed);
      spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus]->Draw("same");

      // legend
      TLegend * leg = new TLegend(0.6,0.64,0.89,0.79);  //coordinates are fractions
      leg->AddEntry(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus],"Generated","p"); 
      leg->AddEntry(kRatioPtTrueData[i],"Corrected","p"); 
      leg->Draw();
    }
   
  cout << "Analysis complete" << endl;
  cout << " - Please exit root to close and save output file " << endl;
}
Commit	Line	Data
c88234ad	1	#include "AliAODTrack.h"
	2	#include <iostream>
	3	#include "TFile.h"
	4	#include "TH1I.h"
	5	#include "TH2.h"
	6	#include "AliSpectraAODHistoManager.h"
	7	#include "AliSpectraAODEventCuts.h"
	8	#include "AliSpectraAODTrackCuts.h"
	9	#include "TCanvas.h"
	10	#include "TLegend.h"
	11
	12	using namespace std;
	13
	14	void analysis_macro(const char * dataFile = "Pt.AOD.1._data_ptcut.root", const char * mcFile ="Pt.AOD.1._MC.root") {
	15
	16	// gSystem->Load("libTree.so");
	17	// gSystem->Load("libGeom.so");
	18	// gSystem->Load("libVMC.so");
	19	// gSystem->Load("libPhysics.so");
	20	// gSystem->Load("libSTEERBase.so");
	21	// gSystem->Load("libESD.so");
	22	// gSystem->Load("libAOD.so");
	23	// gSystem->Load("libANALYSIS.so");
	24	// gSystem->Load("libANALYSISalice.so");
	25	// gSystem->Load("libANALYSIS");
	26	// gSystem->Load("libANALYSISalice");
	27	// gROOT->LoadMacro("$ALICE_ROOT/PWG0/CreateESDChain.C");
	28	// gSystem->AddIncludePath("-I$ALICE_ROOT/include");
	29	// gStyle->SetPalette(1);
	30	// gStyle->SetFillColor(kWhite);
	31
	32	// gROOT->LoadMacro("AliSpectraAODTrackCuts.cxx+g");
	33	// gROOT->LoadMacro("AliSpectraAODEventCuts.cxx+g");
	34	// gROOT->LoadMacro("AliSpectraAODHistoManager.cxx+g");
	35	// gROOT->LoadMacro("AliAnalysisTaskSpectraAOD.cxx+g");
	36
	37	const char * histoname[] =
	38	{
	39	// names of histos we want to draw. please observe order of histos, elsewise name will not make sense
	40	"#frac{MC_{#sigma, rec}}{Data}, P^{+}",
	41	"#frac{MC_{#sigma, rec}}{Data}, K^{+}",
	42	"#frac{MC_{#sigma, rec}}{Data}, #pi^{+}",
	43	"#frac{MC_{#sigma, rec}}{Data}, P^{-}",
	44	"#frac{MC_{#sigma, rec}}{Data}, K^{-}",
	45	"#frac{MC_{#sigma, rec}}{Data}, #pi^{-}",
	46
	47	"#frac{MC_{#sigma, rec, prim}}{Data}, P^{+}",
	48	"#frac{MC_{#sigma, rec, prim}}{Data}, K^{+}",
	49	"#frac{MC_{#sigma, rec, prim}}{Data}, #pi^{+}",
	50	"#frac{MC_{#sigma, rec, prim}}{Data}, P^{-}",
	51	"#frac{MC_{#sigma, rec, prim}}{Data}, K^{-}",
	52	"#frac{MC_{#sigma, rec, prim}}{Data}, #pi^{-}",
	53
	54	"#frac{MC_{#sigma, rec, sec}}{Data}, P^{+}",
	55	"#frac{MC_{#sigma, rec, sec}}{Data}, K^{+}",
	56	"#frac{MC_{#sigma, rec, sec}}{Data}, #pi^{+}",
	57	"#frac{MC_{#sigma, rec, sec}}{Data}, P^{-}",
	58	"#frac{MC_{#sigma, rec, sec}}{Data}, K^{-}",
	59	"#frac{MC_{#sigma, rec, sec}}{Data}, #pi^{-}",
	60
	61	"#frac{MC_{true, gen}}{MC_{#sigma, rec}}, P^{+}",
	62	"#frac{MC_{true, gen}}{MC_{#sigma, rec}}, K^{+}",
	63	"#frac{MC_{true, gen}}{MC_{#sigma, rec}}, #pi^{+}",
	64	"#frac{MC_{true, gen}}{MC_{#sigma, rec}}, P^{-}",
65	"#frac{MC_{true, gen}}{MC_{#sigma, rec}}, K^{-}",
66	"#frac{MC_{true, gen}}{MC_{#sigma, rec}}, #pi^{-}",
67
68	"#frac{MC_{true, gen}}{MC_{#sigma, rec}} x Data, P^{+}",
69	"#frac{MC_{true, gen}}{MC_{#sigma, rec}} x Data, K^{+}",
70	"#frac{MC_{true, gen}}{MC_{#sigma, rec}} x Data, #pi^{+}",
71	"#frac{MC_{true, gen}}{MC_{#sigma, rec}} x Data, P^{-}",
72	"#frac{MC_{true, gen}}{MC_{#sigma, rec}} x Data, K^{-}",
73	"#frac{MC_{true, gen}}{MC_{#sigma, rec}} x Data, #pi^{-}",
74
75	"Data, P^{+}",
76	"Data, K^{+}",
77	"Data, #pi^{+}",
78	"Data, P^{-}",
79	"Data, K^{-}",
80	"Data, #pi^{-}",
81
82	"#frac{MC_{rec #sigma}}{MC_{true, rec}}, P^{+}",
83	"#frac{MC_{rec #sigma}}{MC_{true, rec}}, K^{+}",
84	"#frac{MC_{rec #sigma}}{MC_{true, rec}}, #pi^{+}",
85	"#frac{MC_{rec #sigma}}{MC_{true, rec}}, P^{-}",
86	"#frac{MC_{rec #sigma}}{MC_{true, rec}}, K^{-}",
87	"#frac{MC_{rec #sigma}}{MC_{true, rec}}, #pi^{-}",
88
89	"#frac{MC_{rec #sigma sec}}{MC_{#sigma, rec}}, P^{+}",
90	"#frac{MC_{rec #sigma sec}}{MC_{#sigma, rec}}, K^{+}",
91	"#frac{MC_{rec #sigma sec}}{MC_{#sigma, rec}}, #pi^{+}",
92	"#frac{MC_{rec #sigma sec}}{MC_{#sigma, rec}}, P^{-}",
93	"#frac{MC_{rec #sigma sec}}{MC_{#sigma, rec}}, K^{-}",
94	"#frac{MC_{rec #sigma sec}}{MC_{#sigma, rec}}, #pi^{-}",
95	};
96
97	const char * histotitle[] =
98	{
99	// titles of histos we want to draw. please observe order of histos, elsewise name will not make sense
100	"MC #sigma, rec / Data, P^{+}",
101	"MC #sigma, rec / Data, K^{+}",
102	"MC #sigma, rec / Data, #pi^{+}",
103	"MC #sigma, rec / Data, P^{-}",
104	"MC #sigma, rec / Data, K^{-}",
105	"MC #sigma, rec / Data, #pi^{-}",
106
107	"MC #sigma, rec, prim / Data, P^{+}",
108	"MC #sigma, rec, prim / Data, K^{+}",
109	"MC #sigma, rec, prim / Data, #pi^{+}",
110	"MC #sigma, rec, prim / Data, P^{-}",
111	"MC #sigma, rec, prim / Data, K^{-}",
112	"MC #sigma, rec, prim / Data, #pi^{-}",
113
114	"MC #sigma, rec, sec / Data, P^{+}",
115	"MC #sigma, rec, sec / Data, K^{+}",
116	"MC #sigma, rec, sec / Data, #pi^{+}",
117	"MC #sigma, rec, sec / Data, P^{-}",
118	"MC #sigma, rec, sec / Data, K^{-}",
119	"MC #sigma, rec, sec / Data, #pi^{-}",
120
121	"Correction factor, P^{+}",
122	"Correction factor, K^{+}",
123	"Correction factor, #pi^{+}",
124	"Correction factor, P^{-}",
125	"Correction factor, K^{-}",
126	"Correction factor, #pi^{-}",
127
128	"True Data, P^{+}",
129	"True Data, K^{+}",
130	"True Data, #pi^{+}",
131	"True Data, P^{-}",
132	"True Data, K^{-}",
133	"True Data, #pi^{-}",
134
135	"Raw Data, P^{+}",
136	"Raw Data, K^{+}",
137	"Raw Data, #pi^{+}",
138	"Raw Data, P^{-}",
139	"Raw Data, K^{-}",
140	"Raw Data, #pi^{-}",
141
142	"Identificaiton quality, P^{+}",
143	"Identificaiton quality, K^{+}",
144	"Identificaiton quality, #pi^{+}",
145	"Identificaiton quality, P^{-}",
146	"Identificaiton quality, K^{-}",
147	"Identificaiton quality, #pi^{-}",
148
149	"Contribution of secondaries, P^{+}",
150	"Contribution of secondaries, K^{+}",
151	"Contribution of secondaries, #pi^{+}",
152	"Contribution of secondaries, P^{-}",
153	"Contribution of secondaries, K^{-}",
154	"Contribution of secondaries, #pi^{-}",
155	};
156
157	// Open root MC file and get classes
158	cout << "Analysis Macro" << endl;
159	cout << " > Reading MC data" << endl;
160	TFile *_mc = TFile::Open(mcFile);
161	AliSpectraAODHistoManager* hman_mc = (AliSpectraAODHistoManager*) _mc->Get("SpectraHistos");
162	AliSpectraAODEventCuts* ecuts_mc = (AliSpectraAODEventCuts*) _mc->Get("Event Cuts");
163	AliSpectraAODTrackCuts* tcuts_mc = (AliSpectraAODTrackCuts*) _mc->Get("Track Cuts");
164	// print info about mc track and Event cuts
165	cout << " -- Info about MC -- "<< endl;
166	ecuts_mc->PrintCuts();
167	tcuts_mc->PrintCuts();
168	// get the mc histo's necessary to write the contamination histos from the rootfile
169	TH1F* spectrahistos_mc[30];
170	cout << " -- Reading and normalizing histograms -- " << endl;
171	for( Int_t i = 0 ; i <= AliSpectraNameSpace::kNPtSpecies ; i++ ) {
172	spectrahistos_mc[i] = dynamic_cast<TH1F*>(hman_mc->GetHistogram((AliSpectraNameSpace::AODPtHist_t)i));
173	Double_t events_mc = 1. / ecuts_mc->NumberOfEvents();
174	if(!spectrahistos_mc[i]) {
175	continue;
176	}
177	spectrahistos_mc[i]->Scale(events_mc, "width");
178	}
179	cout << " > Reading real data " << endl;
180	// proceed likewise for data
181	TFile *_data = TFile::Open(dataFile);
182	AliSpectraAODHistoManager* hman_data = (AliSpectraAODHistoManager*) _data->Get("SpectraHistos");
183	AliSpectraAODEventCuts* ecuts_data = (AliSpectraAODEventCuts*) _data->Get("Event Cuts");
184	AliSpectraAODTrackCuts* tcuts_data = (AliSpectraAODTrackCuts*) _data->Get("Track Cuts");
185	// print info about track and Event cuts
186	cout << " -- Info about data -- " << endl;
187	ecuts_data->PrintCuts();
188	tcuts_data->PrintCuts();
189	// get the histo's necessary to write the contamination histos from the rootfile
7df5f6a3	190	// MF 22/02/2012
	191	// Get the raw data distribution:
	192	// 1. Get DCA vs pt histo for each particle
	193	// 2. Loop over all pt bins, projecting on DCA
	194	// 3. Fit the DCA data distribution with MC templates for primaries, sec. w., sec. m. (TFractionFitter)
	195	// 4. Extract from the fit raw primaries vs pt in data
	196	// 5. Correct for the primaries efficiency estimated in MC
c88234ad	197	TH1F* spectrahistos_data[30];
	198	cout << " -- Reading and normalizing histograms -- " << endl;
	199	for( Int_t i = 0 ; i <= AliSpectraNameSpace::kNPtSpecies ; i++ ) {
	200	spectrahistos_data[i] = dynamic_cast<TH1F*>(hman_data->GetHistogram((AliSpectraNameSpace::AODPtHist_t)i));
	201	if(!spectrahistos_data[i]) {
	202	continue;
	203	}
	204	Double_t events_data = 1. / ecuts_data->NumberOfEvents();
	205	//normalize the histos
	206	spectrahistos_data[i]->Scale(events_data, "width");
	207	}
	208	//create output file and select it as the active file
	209	TFile * output = new TFile("analysis_output.root", "RECREATE");
	210	output->cd();
	211	// Write the scaled data and MC histos in the output file for convenience
	212	output->mkdir("Data");
	213	output->cd("Data");
	214	for(Int_t ihist = 0; ihist <= AliSpectraNameSpace::kNPtSpecies; ihist++){
	215	spectrahistos_data[ihist]->Write();
	216	}
	217	output->mkdir("MC");
	218	output->cd("MC");
	219	for(Int_t ihist = 0; ihist <= AliSpectraNameSpace::kNPtSpecies; ihist++){
	220	spectrahistos_mc[ihist]->Write();
	221	}
	222	output->mkdir("ANALYSIS");
	223	output->cd("ANALYSIS");
	224	// perform some analysis!
	225	// next loop: MC_Sigma_Rec over Rec_Data
	226	TCanvas* c1 = new TCanvas("MC reconstructed over Data","MC reconstructed over Data");
	227	c1->Divide(3,2);
	228	TH1F* kRatioPtMCSigmaRecSigmaRec[6];
	229	for (Int_t i = 0 ; i < 6 ; i ++ )
	230	{
	231	kRatioPtMCSigmaRecSigmaRec[i] = (TH1F*)(spectrahistos_mc[i])->Clone(histotitle[i]);
	232	kRatioPtMCSigmaRecSigmaRec[i]->Divide(spectrahistos_data[i]);
	233	c1->cd(1+i);
	234	kRatioPtMCSigmaRecSigmaRec[i]->SetTitle(histoname[i]);
	235	kRatioPtMCSigmaRecSigmaRec[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
	236	kRatioPtMCSigmaRecSigmaRec[i]->Draw();
	237	kRatioPtMCSigmaRecSigmaRec[i]->Write();
	238	}
	239	// next loop: MC_sigma_primaries over Rec_data
	240	TCanvas* c2 = new TCanvas("MC reconstructed primaries over Data", "MC reconstructed primaries over Data");
	241	c2->Divide(3,2);
	242	TH1F* kRatioPtMCSigmaRecPrimarySigmaRec[6];
	243	for (Int_t i = 0 ; i < 6 ; i ++ )
	244	{
	245	kRatioPtMCSigmaRecPrimarySigmaRec[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaPrimaryProtonPlus])->Clone(histotitle[i+6]);
	246	kRatioPtMCSigmaRecPrimarySigmaRec[i]->Divide(spectrahistos_data[i]);
	247	c2->cd(1+i);
	248	kRatioPtMCSigmaRecPrimarySigmaRec[i]->SetTitle(histoname[i+6]);
	249	kRatioPtMCSigmaRecPrimarySigmaRec[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
	250	kRatioPtMCSigmaRecPrimarySigmaRec[i]->Draw();
	251	kRatioPtMCSigmaRecPrimarySigmaRec[i]->Write();
	252	}
	253	// next loop: MC_Sigma_secondaries over Rec_data
	254	TCanvas* c3 = new TCanvas("MC reconstructed secondaries over Data","MC reconstructed secondaries over Data");
	255	c3->Divide(3,2);
	256	TH1F* kRatioPtMCSigmaRecSigmaRecSecondary[6];
	257	for (Int_t i = 0 ; i < 6 ; i ++ )
	258	{
	259	kRatioPtMCSigmaRecSigmaRecSecondary[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaSecondaryProtonPlus])->Clone(histotitle[i+12]);
	260	kRatioPtMCSigmaRecSigmaRecSecondary[i]->Divide(spectrahistos_data[i]);
261	c3->cd(1+i);
262	kRatioPtMCSigmaRecSigmaRecSecondary[i]->SetTitle(histoname[12+i]);
263	kRatioPtMCSigmaRecSigmaRecSecondary[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
264	kRatioPtMCSigmaRecSigmaRecSecondary[i]->Draw();
265	kRatioPtMCSigmaRecSigmaRecSecondary[i]->Write();
266	}
267	//
268	TCanvas* c4 = new TCanvas("MC true generated over MC true reconstructed","MC true generated over MC true reconstructed");
269	c4->Divide(3,2);
270	TH1F* kRatioPtMCTrueMCRec[6];
271	for (Int_t i = 0 ; i < 6 ; i ++ )
272	{
273	kRatioPtMCTrueMCRec[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus])->Clone(histotitle[i+18]);
274	kRatioPtMCTrueMCRec[i]->Divide(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaPrimaryProtonPlus]);
275	c4->cd(1+i);
276	kRatioPtMCTrueMCRec[i]->SetTitle(histoname[18+i]);
277	kRatioPtMCTrueMCRec[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
278	kRatioPtMCTrueMCRec[i]->Draw();
279	kRatioPtMCTrueMCRec[i]->Write();
280	}
281
282	TCanvas* c5 = new TCanvas("True Data","True Data");
283	c5->Divide(3,2);
284	TH1F* kRatioPtTrueData[6];
285	for (Int_t i = 0 ; i < 6 ; i ++ )
286	{
287	kRatioPtTrueData[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus])->Clone(histotitle[i+24]);
288	kRatioPtTrueData[i]->Divide(spectrahistos_mc[i]); // Michele: Since we are correcting all contributions at once, this should be exaclty the same as in the data!
289	kRatioPtTrueData[i]->Multiply(spectrahistos_data[i]);
290	c5->cd(1+i);
291	kRatioPtTrueData[i]->SetTitle(histoname[24+i]);
292	kRatioPtTrueData[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
293	kRatioPtTrueData[i]->Draw();
294	kRatioPtTrueData[i]->Write();
295	}
296
297	TCanvas* c6 = new TCanvas("Raw data", "Raw data");
298	c6->Divide(3,2);
299	TH1F* kRawData[6];
300	for (Int_t i = 0 ; i < 6 ; i ++ )
301	{
302	kRawData[i] = (TH1F*)(spectrahistos_data[i+AliSpectraNameSpace::kHistPtRecSigmaProtonPlus])->Clone(histotitle[i+30]);
303	c6->cd(1+i);
304	kRawData[i]->SetTitle(histoname[30+i]);
305	kRawData[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
306	kRawData[i]->Draw();
307	kRawData[i]->Write();
308	}
309
310	TCanvas* c7 = new TCanvas("contamination", "contamination");
311	c7->Divide(3,2);
312	TH1F* kContamination[6];
313	for (Int_t i = 0 ; i < 6 ; i ++ )
314	{
315	kContamination[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaProtonPlus])->Clone(histotitle[i+36]);
316	kContamination[i]->Divide(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecTrueProtonPlus]);
317	c7->cd(1+i);
318	kContamination[i]->SetTitle(histoname[36+i]);
319	kContamination[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
320	kContamination[i]->Draw();
321	kContamination[i]->Write();
322	}
323
324	TCanvas* c8 = new TCanvas("contamination secondary", "contamination secondary");
325	c8->Divide(3,2);
326	TH1F* kSecContamination[6];
327	for (Int_t i = 0 ; i < 6 ; i ++ )
328	{
329	kSecContamination[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaSecondaryProtonPlus])->Clone(histotitle[i+42]);
330	kSecContamination[i]->Divide(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaProtonPlus]);
331	c8->cd(1+i);
332	kSecContamination[i]->SetTitle(histoname[42+i]);
333	kSecContamination[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
334	kSecContamination[i]->Draw();
335	kSecContamination[i]->Write();
336	}
337
338	TCanvas* c9 = new TCanvas("true data and generated", "true data and generated");
339	c9->Divide(3,2);
340	TH1F* ktrueandgenerated[6];
341	for (Int_t i = 0 ; i < 6 ; i ++ )
342	{
343	// kSecContamination[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaSecondaryProtonPlus])->Clone(histotitle[i+42]);
344	// kSecContamination[i]->Divide(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaProtonPlus]);
345	c9->cd(1+i);
346	// kSecContamination[i]->SetTitle(histoname[42+i]);
347	// kSecContamination[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
348	kRatioPtTrueData[i]->Draw();
349	spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus]->SetMarkerColor(kRed);
350	spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus]->SetMarkerStyle(kOpenCircle);
351	spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus]->SetLineColor(kRed);
352	spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus]->Draw("same");
353
354	// legend
355	TLegend * leg = new TLegend(0.6,0.64,0.89,0.79); //coordinates are fractions
356	leg->AddEntry(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus],"Generated","p");
357	leg->AddEntry(kRatioPtTrueData[i],"Corrected","p");
358	leg->Draw();
359	}
360
361	cout << "Analysis complete" << endl;
362	cout << " - Please exit root to close and save output file " << endl;
363	}