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