PWGLF/SPECTRA/PiKaPr/TestAOD/analysis_macro.C

   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
 190   TH1F* spectrahistos_data[30];
 191   cout << " -- Reading and normalizing histograms -- " << endl;
 192   for( Int_t i = 0 ; i <= AliSpectraNameSpace::kNPtSpecies ; i++ ) {
 193     spectrahistos_data[i] = dynamic_cast<TH1F*>(hman_data->GetHistogram((AliSpectraNameSpace::AODPtHist_t)i));
 194     if(!spectrahistos_data[i]) {
 195       continue;
 196     }
 197     Double_t events_data =  1. / ecuts_data->NumberOfEvents();
 198     //normalize the histos
 199     spectrahistos_data[i]->Scale(events_data, "width");
 200   }
 201   //create output file and select it as the active file
 202   TFile * output = new TFile("analysis_output.root", "RECREATE");
 203   output->cd();
 204   // Write the scaled data and MC histos in the output file for convenience
 205   output->mkdir("Data");
 206   output->cd("Data");
 207   for(Int_t ihist = 0; ihist <= AliSpectraNameSpace::kNPtSpecies; ihist++){
 208       spectrahistos_data[ihist]->Write();
 209   }
 210   output->mkdir("MC");
 211   output->cd("MC");
 212   for(Int_t ihist = 0; ihist <= AliSpectraNameSpace::kNPtSpecies; ihist++){
 213       spectrahistos_mc[ihist]->Write();
 214   }
 215   output->mkdir("ANALYSIS");
 216   output->cd("ANALYSIS");
 217   // perform some analysis!
 218   //  next loop: MC_Sigma_Rec over Rec_Data
 219   TCanvas* c1 = new TCanvas("MC reconstructed over Data","MC reconstructed over Data");
 220   c1->Divide(3,2);
 221   TH1F* kRatioPtMCSigmaRecSigmaRec[6];
 222   for (Int_t i = 0 ; i < 6 ; i ++ )
 223     {
 224       kRatioPtMCSigmaRecSigmaRec[i] = (TH1F*)(spectrahistos_mc[i])->Clone(histotitle[i]);
 225       kRatioPtMCSigmaRecSigmaRec[i]->Divide(spectrahistos_data[i]);
 226       c1->cd(1+i);
 227       kRatioPtMCSigmaRecSigmaRec[i]->SetTitle(histoname[i]);
 228       kRatioPtMCSigmaRecSigmaRec[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
 229       kRatioPtMCSigmaRecSigmaRec[i]->Draw();
 230       kRatioPtMCSigmaRecSigmaRec[i]->Write();
 231     }
 232   // next loop: MC_sigma_primaries over Rec_data
 233   TCanvas* c2 = new TCanvas("MC reconstructed primaries over Data", "MC reconstructed primaries over Data");
 234   c2->Divide(3,2);
 235   TH1F* kRatioPtMCSigmaRecPrimarySigmaRec[6];
 236   for (Int_t i = 0 ; i < 6 ; i ++ )
 237     {
 238       kRatioPtMCSigmaRecPrimarySigmaRec[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaPrimaryProtonPlus])->Clone(histotitle[i+6]);
 239       kRatioPtMCSigmaRecPrimarySigmaRec[i]->Divide(spectrahistos_data[i]);
 240       c2->cd(1+i);
 241       kRatioPtMCSigmaRecPrimarySigmaRec[i]->SetTitle(histoname[i+6]);
 242       kRatioPtMCSigmaRecPrimarySigmaRec[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
 243       kRatioPtMCSigmaRecPrimarySigmaRec[i]->Draw();
 244       kRatioPtMCSigmaRecPrimarySigmaRec[i]->Write();
 245     }
 246   // next loop: MC_Sigma_secondaries over Rec_data
 247   TCanvas* c3 = new TCanvas("MC reconstructed secondaries over Data","MC reconstructed secondaries over Data");
 248   c3->Divide(3,2);
 249   TH1F* kRatioPtMCSigmaRecSigmaRecSecondary[6];
 250   for (Int_t i = 0 ; i < 6 ; i ++ )
 251     {
 252       kRatioPtMCSigmaRecSigmaRecSecondary[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaSecondaryProtonPlus])->Clone(histotitle[i+12]);
 253       kRatioPtMCSigmaRecSigmaRecSecondary[i]->Divide(spectrahistos_data[i]);
 254       c3->cd(1+i);
 255       kRatioPtMCSigmaRecSigmaRecSecondary[i]->SetTitle(histoname[12+i]);
 256       kRatioPtMCSigmaRecSigmaRecSecondary[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
 257       kRatioPtMCSigmaRecSigmaRecSecondary[i]->Draw();
 258       kRatioPtMCSigmaRecSigmaRecSecondary[i]->Write();
 259     }
 260   //
 261   TCanvas* c4 = new TCanvas("MC true generated over MC true reconstructed","MC true generated over MC true reconstructed");
 262   c4->Divide(3,2);
 263   TH1F* kRatioPtMCTrueMCRec[6];
 264   for (Int_t i = 0 ; i < 6 ; i ++ )
 265     {
 266       kRatioPtMCTrueMCRec[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus])->Clone(histotitle[i+18]);
 267       kRatioPtMCTrueMCRec[i]->Divide(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaPrimaryProtonPlus]);
 268       c4->cd(1+i);
 269       kRatioPtMCTrueMCRec[i]->SetTitle(histoname[18+i]);
 270       kRatioPtMCTrueMCRec[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
 271       kRatioPtMCTrueMCRec[i]->Draw();
 272       kRatioPtMCTrueMCRec[i]->Write();
 273     }
 274
 275   TCanvas* c5 = new TCanvas("True Data","True Data");
 276   c5->Divide(3,2);
 277   TH1F* kRatioPtTrueData[6];
 278   for (Int_t i = 0 ; i < 6 ; i ++ )
 279     {
 280       kRatioPtTrueData[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus])->Clone(histotitle[i+24]);
 281       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!
 282       kRatioPtTrueData[i]->Multiply(spectrahistos_data[i]);
 283       c5->cd(1+i);
 284       kRatioPtTrueData[i]->SetTitle(histoname[24+i]);
 285       kRatioPtTrueData[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
 286       kRatioPtTrueData[i]->Draw();
 287       kRatioPtTrueData[i]->Write();
 288     }
 289
 290   TCanvas* c6 = new TCanvas("Raw data", "Raw data");
 291   c6->Divide(3,2);
 292   TH1F* kRawData[6];
 293   for (Int_t i = 0 ; i < 6 ; i ++ )
 294     {
 295       kRawData[i] = (TH1F*)(spectrahistos_data[i+AliSpectraNameSpace::kHistPtRecSigmaProtonPlus])->Clone(histotitle[i+30]);
 296       c6->cd(1+i);
 297       kRawData[i]->SetTitle(histoname[30+i]);
 298       kRawData[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
 299       kRawData[i]->Draw();
 300       kRawData[i]->Write();
 301     }
 302
 303   TCanvas* c7 = new TCanvas("contamination", "contamination");
 304   c7->Divide(3,2);
 305   TH1F* kContamination[6];
 306   for (Int_t i = 0 ; i < 6 ; i ++ )
 307     {
 308       kContamination[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaProtonPlus])->Clone(histotitle[i+36]);
 309       kContamination[i]->Divide(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecTrueProtonPlus]);
 310       c7->cd(1+i);
 311       kContamination[i]->SetTitle(histoname[36+i]);
 312       kContamination[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
 313       kContamination[i]->Draw();
 314       kContamination[i]->Write();
 315     }
 316
 317   TCanvas* c8 = new TCanvas("contamination secondary", "contamination secondary");
 318   c8->Divide(3,2);
 319   TH1F* kSecContamination[6];
 320   for (Int_t i = 0 ; i < 6 ; i ++ )
 321     {
 322       kSecContamination[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaSecondaryProtonPlus])->Clone(histotitle[i+42]);
 323       kSecContamination[i]->Divide(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaProtonPlus]);
 324       c8->cd(1+i);
 325       kSecContamination[i]->SetTitle(histoname[42+i]);
 326       kSecContamination[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
 327       kSecContamination[i]->Draw();
 328       kSecContamination[i]->Write();
 329     }
 330
 331   TCanvas* c9 = new TCanvas("true data and generated", "true data and generated");
 332   c9->Divide(3,2);
 333   TH1F* ktrueandgenerated[6];
 334   for (Int_t i = 0 ; i < 6 ; i ++ )
 335     {
 336       // kSecContamination[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaSecondaryProtonPlus])->Clone(histotitle[i+42]);
 337       // kSecContamination[i]->Divide(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaProtonPlus]);
 338       c9->cd(1+i);
 339       // kSecContamination[i]->SetTitle(histoname[42+i]);
 340       // kSecContamination[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
 341       kRatioPtTrueData[i]->Draw();
 342       spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus]->SetMarkerColor(kRed);
 343       spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus]->SetMarkerStyle(kOpenCircle);
 344       spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus]->SetLineColor(kRed);
 345       spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus]->Draw("same");
 346
 347       // legend
 348       TLegend * leg = new TLegend(0.6,0.64,0.89,0.79);  //coordinates are fractions
 349       leg->AddEntry(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus],"Generated","p");
 350       leg->AddEntry(kRatioPtTrueData[i],"Corrected","p");
 351       leg->Draw();
 352     }
 353
 354   cout << "Analysis complete" << endl;
 355   cout << " - Please exit root to close and save output file " << endl;
 356 }