PWGJE/macros/PID/extractGenMCtruth.C

   1 #include "TH1D.h"
   2 #include "TH2D.h"
   3 #include "TCanvas.h"
   4 #include "TFile.h"
   5 #include "TLegend.h"
   6
   7 #include "AliPID.h"
   8
   9 #include <iostream>
  10
  11 #include "THnSparseDefinitions.h"
  12
  13 enum type { kTrackPt = 0, kZ = 1, kXi = 2, kNtypes = 3 };
  14
  15 //___________________________________________________________________
  16 void normaliseHist(TH1* h, Double_t scale = 1.0)
  17 {
  18   if (h->GetSumw2N() <= 0)
  19     h->Sumw2();
  20
  21   h->Scale(scale);
  22
  23   for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
  24     Double_t normFactor = h->GetBinWidth(i);
  25     h->SetBinContent(i, h->GetBinContent(i) / normFactor);
  26     h->SetBinError(i, h->GetBinError(i) / normFactor);
  27   }
  28 }
  29
  30
  31 //___________________________________________________________________
  32 Int_t extractGenMCtruth(TString pathNameData, TString listName,
  33                         Int_t chargeMode /*kNegCharge = -1, kAllCharged = 0, kPosCharge = 1*/,
  34                         Double_t lowerCentrality /*= -2*/, Double_t upperCentrality /*= -2*/,
  35                         Double_t lowerJetPt /*= -1*/ , Double_t upperJetPt/* = -1*/)
  36 {
  37   if (listName == "") {
  38     listName = pathNameData;
  39     listName.Replace(0, listName.Last('/') + 1, "");
  40     listName.ReplaceAll(".root", "");
  41   }
  42
  43   TString titles[kNtypes];
  44   titles[kTrackPt] = "trackPt";
  45   titles[kZ] = "z";
  46   titles[kXi] = "xi";
  47
  48   TString pathData = pathNameData;
  49   pathData.Replace(pathData.Last('/'), pathData.Length(), "");
  50
  51
  52   TH2D* hNjetsGenData = 0x0;
  53   TH2D* hNjetsRecData = 0x0;
  54
  55   TH1D* hMCgenPrimYieldPt[AliPID::kSPECIES] = {0x0, };
  56   TH1D* hMCgenPrimYieldZ[AliPID::kSPECIES] = {0x0, };
  57   TH1D* hMCgenPrimYieldXi[AliPID::kSPECIES] = {0x0, };
  58
  59   TFile* fileData = TFile::Open(pathNameData.Data());
  60   if (!fileData) {
  61     printf("Failed to open data file \"%s\"\n", pathNameData.Data());
  62     return -1;
  63   }
  64
  65   TObjArray* histList = (TObjArray*)(fileData->Get(listName.Data()));
  66   THnSparse* hMCgeneratedYieldsPrimaries = (THnSparse*)histList->FindObject("fhMCgeneratedYieldsPrimaries");
  67
  68   if (!hMCgeneratedYieldsPrimaries) {
  69     printf("Failed to load generated primary yields!\n");
  70     return -1;
  71   }
  72
  73   // Set proper errors, if not yet calculated
  74   if (!hMCgeneratedYieldsPrimaries->GetCalculateErrors()) {
  75     std::cout << "Re-calculating errors of " << hMCgeneratedYieldsPrimaries->GetName() << "..." << std::endl;
  76
  77     hMCgeneratedYieldsPrimaries->Sumw2();
  78
  79     Long64_t nBinsTHnSparseGenYield = hMCgeneratedYieldsPrimaries->GetNbins();
  80     Double_t binContentGenYield = 0;
  81     for (Long64_t bin = 0; bin < nBinsTHnSparseGenYield; bin++) {
  82       binContentGenYield = hMCgeneratedYieldsPrimaries->GetBinContent(bin);
  83       hMCgeneratedYieldsPrimaries->SetBinError(bin, TMath::Sqrt(binContentGenYield));
  84     }
  85   }
  86
  87   hNjetsGenData = (TH2D*)histList->FindObject("fh2FFJetPtGen");
  88
  89   hNjetsRecData = (TH2D*)histList->FindObject("fh2FFJetPtRec");
  90
  91   Bool_t restrictJetPtAxis = (lowerJetPt >= 0 && upperJetPt >= 0);
  92
  93   if (restrictJetPtAxis && (!hNjetsRecData || !hNjetsGenData)) {
  94     printf("Failed to load numJet histo for data!\n");
  95
  96     // For backward compatibility (TODO REMOVE IN FUTURE): Load info from fixed AnalysisResults file (might be wrong, if other
  97     // period is considered; also: No multiplicity information)
  98
  99     TString pathBackward = Form("%s/AnalysisResults.root", pathData.Data());
 100     TFile* fBackward = TFile::Open(pathBackward.Data());
 101
 102     TString dirDataInFile = "";
 103     TDirectory* dirData = fBackward ? (TDirectory*)fBackward->Get(fBackward->GetListOfKeys()->At(0)->GetName()) : 0x0;
 104
 105     TList* list = dirData ? (TList*)dirData->Get(dirData->GetListOfKeys()->At(0)->GetName()) : 0x0;
 106
 107     TH1D* hFFJetPtRec = list ? (TH1D*)list->FindObject("fh1FFJetPtRecCutsInc") : 0x0;
 108     TH1D* hFFJetPtGen = list ? (TH1D*)list->FindObject("fh1FFJetPtGenInc") : 0x0;
 109
 110     if (hFFJetPtRec && hFFJetPtGen) {
 111       printf("***WARNING: For backward compatibility, using file \"%s\" to get number of jets. BUT: Might be wrong period and has no mult info!***\n",
 112         pathBackward.Data());
 113       printf("ALSO: Using Njets for inclusive jets!!!!\n");
 114
 115       hNjetsRecData = new TH2D("fh2FFJetPtRec", "", 1, -1, 1, hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldJetPt)->GetNbins(),
 116                                hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldJetPt)->GetXbins()->GetArray());
 117
 118       for (Int_t iJet = 1; iJet <= hNjetsRecData->GetNbinsY(); iJet++) {
 119         Int_t lowerBin = hFFJetPtRec->FindBin(hNjetsRecData->GetYaxis()->GetBinLowEdge(iJet) + 1e-3);
 120         Int_t upperBin = hFFJetPtRec->FindBin(hNjetsRecData->GetYaxis()->GetBinUpEdge(iJet) - 1e-3);
 121         hNjetsRecData->SetBinContent(1, iJet, hFFJetPtRec->Integral(lowerBin, upperBin));
 122       }
 123
 124       hNjetsGenData = new TH2D("fh2FFJetPtGen", "", 1, -1, 1,  hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldJetPt)->GetNbins(),
 125                                hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldJetPt)->GetXbins()->GetArray());
 126
 127       for (Int_t iJet = 1; iJet <= hNjetsGenData->GetNbinsY(); iJet++) {
 128         Int_t lowerBin = hFFJetPtGen->FindBin(hNjetsGenData->GetYaxis()->GetBinLowEdge(iJet) + 1e-3);
 129         Int_t upperBin = hFFJetPtGen->FindBin(hNjetsGenData->GetYaxis()->GetBinUpEdge(iJet) - 1e-3);
 130         hNjetsGenData->SetBinContent(1, iJet, hFFJetPtGen->Integral(lowerBin, upperBin));
 131       }
 132     }
 133
 134     if (!hNjetsRecData || ! hNjetsGenData)
 135       return -1;
 136   }
 137
 138   const Bool_t restrictCentralityData = ((lowerCentrality >= -1) && (upperCentrality >= -1));
 139   // Integral(lowerCentBinLimitData, uppCentBinLimitData) will not be restricted if these values are kept
 140   const Int_t lowerCentralityBinLimitData = restrictCentralityData ? hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldCentrality)->FindBin(lowerCentrality + 0.001)
 141                                                                    : -1;
 142   const Int_t upperCentralityBinLimitData = restrictCentralityData ? hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldCentrality)->FindBin(upperCentrality - 0.001)
 143                                                                    : -2;
 144
 145   const Double_t actualLowerCentralityData = hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldCentrality)->GetBinLowEdge(lowerCentralityBinLimitData);
 146
 147   const Double_t actualUpperCentralityData = hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldCentrality)->GetBinLowEdge(upperCentralityBinLimitData);
 148
 149   // Integral(lowerJetBinLimitData, uppJetBinLimitData) will not be restricted if these values are kept
 150   const Int_t lowerJetPtBinLimitData = restrictJetPtAxis
 151                                        ? hNjetsRecData->GetYaxis()->FindBin(lowerJetPt + 0.001) : -1;
 152   const Int_t upperJetPtBinLimitData  = restrictJetPtAxis
 153                                        ? hNjetsRecData->GetYaxis()->FindBin(upperJetPt - 0.001) : -2;
 154
 155   const Double_t actualLowerJetPtData = restrictJetPtAxis ? hNjetsRecData->GetYaxis()->GetBinLowEdge(lowerJetPtBinLimitData) : -1;
 156   const Double_t actualUpperJetPtData = restrictJetPtAxis ? hNjetsRecData->GetYaxis()->GetBinUpEdge(upperJetPtBinLimitData)  : -1;
 157
 158   if (restrictCentralityData)
 159     hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldCentrality)->SetRange(lowerCentralityBinLimitData, upperCentralityBinLimitData);
 160
 161   const Bool_t restrictCharge = (chargeMode != kAllCharged);
 162
 163   Int_t lowerChargeBinLimitGenYield = -1;
 164   Int_t upperChargeBinLimitGenYield = -2;
 165
 166   if (restrictCharge) {
 167     const Int_t indexChargeAxisGenYield = GetAxisByTitle(hMCgeneratedYieldsPrimaries, "Charge (e_{0})");
 168     if (indexChargeAxisGenYield < 0) {
 169       std::cout << "Error: Charge axis not found for gen yield histogram!" << std::endl;
 170       return -1;
 171     }
 172
 173     // Add subtract a very small number to avoid problems with values right on the border between to bins
 174     if (chargeMode == kNegCharge) {
 175       lowerChargeBinLimitGenYield = hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->FindBin(-1. + 0.001);
 176       upperChargeBinLimitGenYield = hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->FindBin(0. - 0.001);
 177     }
 178     else if (chargeMode == kPosCharge) {
 179       lowerChargeBinLimitGenYield = hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->FindBin(0. + 0.001);
 180       upperChargeBinLimitGenYield = hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->FindBin(1. - 0.001);
 181     }
 182
 183     // Check if the values look reasonable
 184     if (lowerChargeBinLimitGenYield <= upperChargeBinLimitGenYield && lowerChargeBinLimitGenYield >= 1
 185         && upperChargeBinLimitGenYield <= hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->GetNbins()) {
 186       // OK
 187     }
 188     else {
 189       std::cout << std::endl;
 190       std::cout << "Requested charge range (gen yield) out of limits or upper and lower limit are switched!" << std::endl;
 191       return -1;
 192     }
 193
 194     hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->SetRange(lowerChargeBinLimitGenYield, upperChargeBinLimitGenYield);
 195   }
 196
 197   // If desired, restrict jetPt axis
 198   Int_t lowerJetPtBinLimit = -1;
 199   Int_t upperJetPtBinLimit = -1;
 200   Double_t actualLowerJetPt = -1.;
 201   Double_t actualUpperJetPt = -1.;
 202
 203   if (restrictJetPtAxis) {
 204     // Add subtract a very small number to avoid problems with values right on the border between to bins
 205     lowerJetPtBinLimit = hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldJetPt)->FindBin(lowerJetPt + 0.001);
 206     upperJetPtBinLimit = hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldJetPt)->FindBin(upperJetPt - 0.001);
 207
 208     // Check if the values look reasonable
 209     if (lowerJetPtBinLimit <= upperJetPtBinLimit && lowerJetPtBinLimit >= 1 &&
 210         upperJetPtBinLimit <= hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldJetPt)->GetNbins()) {
 211       actualLowerJetPt = hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldJetPt)->GetBinLowEdge(lowerJetPtBinLimit);
 212       actualUpperJetPt = hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldJetPt)->GetBinUpEdge(upperJetPtBinLimit);
 213
 214       restrictJetPtAxis = kTRUE;
 215
 216       if (TMath::Abs(actualLowerJetPt - actualLowerJetPtData) > 1e-3 || TMath::Abs(actualUpperJetPt - actualUpperJetPtData) > 1e-3) {
 217         std::cout << std::endl;
 218         std::cout << "Mismatch between jet pT range of gen histo and num jet histo!" << std::endl;
 219         return -1;
 220       }
 221     }
 222     else {
 223       std::cout << std::endl;
 224       std::cout << "Requested jet pT range out of limits or upper and lower limit are switched!" << std::endl;
 225       return -1;
 226     }
 227   }
 228
 229   std::cout << "jet pT: ";
 230   if (restrictJetPtAxis) {
 231     std::cout << actualLowerJetPt << " - " << actualUpperJetPt << std::endl;
 232     hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldJetPt)->SetRange(lowerJetPtBinLimit, upperJetPtBinLimit);
 233   }
 234   else {
 235     std::cout << "All" << std::endl;
 236   }
 237
 238
 239
 240   for (Int_t MCid = 0; MCid < AliPID::kSPECIES; MCid++) {
 241     hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldMCpid)->SetRange(MCid + 1, MCid + 1);
 242
 243     hMCgenPrimYieldPt[MCid] = hMCgeneratedYieldsPrimaries->Projection(kPidGenYieldPt, "e");
 244     hMCgenPrimYieldPt[MCid]->SetName(Form("hMCgenYieldsPrimPt_%s", AliPID::ParticleShortName(MCid)));
 245     hMCgenPrimYieldPt[MCid]->SetTitle(Form("MC truth generated primary yield, %s", AliPID::ParticleName(MCid)));
 246     hMCgenPrimYieldPt[MCid]->GetYaxis()->SetTitle("1/N_{Jets} dN/dP_{T} (GeV/c)^{-1}");
 247     hMCgenPrimYieldPt[MCid]->SetStats(kFALSE);
 248
 249     if (hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldZ)) { // For backward compatibility
 250       hMCgenPrimYieldZ[MCid] = hMCgeneratedYieldsPrimaries->Projection(kPidGenYieldZ, "e");
 251       hMCgenPrimYieldZ[MCid]->SetName(Form("hMCgenYieldsPrimZ_%s", AliPID::ParticleShortName(MCid)));
 252       hMCgenPrimYieldZ[MCid]->SetTitle(Form("MC truth generated primary yield, %s", AliPID::ParticleName(MCid)));
 253       hMCgenPrimYieldZ[MCid]->GetYaxis()->SetTitle("1/N_{Jets} dN/dz");
 254       hMCgenPrimYieldZ[MCid]->SetStats(kFALSE);
 255     }
 256
 257     if (hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldXi)) {
 258       hMCgenPrimYieldXi[MCid] = hMCgeneratedYieldsPrimaries->Projection(kPidGenYieldXi, "e");
 259       hMCgenPrimYieldXi[MCid]->SetName(Form("hMCgenYieldsPrimXi_%s", AliPID::ParticleShortName(MCid)));
 260       hMCgenPrimYieldXi[MCid]->SetTitle(Form("MC truth generated primary yield, %s", AliPID::ParticleName(MCid)));
 261       hMCgenPrimYieldXi[MCid]->GetYaxis()->SetTitle("1/N_{Jets} dN/d#xi");
 262       hMCgenPrimYieldXi[MCid]->SetStats(kFALSE);
 263     }
 264
 265     hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldMCpid)->SetRange(0, -1);
 266   }
 267
 268   // Obtain uncorrected fractions
 269   THnSparse* hPIDdata = dynamic_cast<THnSparse*>(histList->FindObject("hPIDdataAll"));
 270   if (!hPIDdata) {
 271     std::cout << std::endl;
 272     std::cout << "Failed to load data histo!" << std::endl;
 273     return -1;
 274   }
 275
 276   // Set proper errors, if not yet calculated
 277   if (!hPIDdata->GetCalculateErrors()) {
 278     std::cout << "Re-calculating errors of " << hPIDdata->GetName() << "..." << std::endl;
 279     hPIDdata->Sumw2();
 280     Long64_t nBinsTHnSparse = hPIDdata->GetNbins();
 281     Double_t binContent = 0;
 282
 283     for (Long64_t bin = 0; bin < nBinsTHnSparse; bin++) {
 284       binContent = hPIDdata->GetBinContent(bin);
 285       hPIDdata->SetBinError(bin, TMath::Sqrt(binContent));
 286     }
 287   }
 288
 289   // Bin limits are the same as in gen yield histo
 290   if (restrictCentralityData)
 291     hPIDdata->GetAxis(kPidCentrality)->SetRange(lowerCentralityBinLimitData, upperCentralityBinLimitData);
 292
 293   if (restrictJetPtAxis)
 294     hPIDdata->GetAxis(kPidJetPt)->SetRange(lowerJetPtBinLimit, upperJetPtBinLimit);
 295
 296   const Int_t indexChargeAxisData = GetAxisByTitle(hPIDdata, "Charge (e_{0})");
 297   if (indexChargeAxisData < 0 && restrictCharge) {
 298     std::cout << "Error: Charge axis not found for data histogram!" << std::endl;
 299     return -1;
 300   }
 301
 302   if (restrictCharge)
 303     hPIDdata->GetAxis(indexChargeAxisData)->SetRange(lowerChargeBinLimitGenYield, upperChargeBinLimitGenYield);
 304
 305
 306   hPIDdata->GetAxis(kPidSelectSpecies)->SetRange(1, 1); // Do not count each particle more than once
 307
 308
 309   TH1D* hMCuncorrYieldPt[AliPID::kSPECIES] = {0x0, };
 310   TH1D* hMCuncorrYieldZ[AliPID::kSPECIES] = {0x0, };
 311   TH1D* hMCuncorrYieldXi[AliPID::kSPECIES] = {0x0, };
 312
 313   for (Int_t MCid = 0; MCid < AliPID::kSPECIES; MCid++) {
 314     // Unfortunately, different MC bins than AliPID
 315     Int_t MCbinInHisto = 0;
 316
 317     switch (MCid) {
 318       case AliPID::kElectron:
 319         MCbinInHisto = 1;
 320         break;
 321       case AliPID::kKaon:
 322         MCbinInHisto = 2;
 323         break;
 324       case AliPID::kMuon:
 325         MCbinInHisto = 3;
 326         break;
 327       case AliPID::kPion:
 328         MCbinInHisto = 4;
 329         break;
 330       case AliPID::kProton:
 331         MCbinInHisto = 5;
 332         break;
 333       default:
 334         break;
 335     }
 336
 337     hPIDdata->GetAxis(kPidMCpid)->SetRange(MCbinInHisto, MCbinInHisto);
 338
 339     hMCuncorrYieldPt[MCid] = hPIDdata->Projection(kPidPt, "e");
 340     hMCuncorrYieldPt[MCid]->SetName(Form("hMCuncorrYieldsPt%s", AliPID::ParticleShortName(MCid)));
 341     hMCuncorrYieldPt[MCid]->SetTitle(Form("MC truth uncorrected yield, %s", AliPID::ParticleName(MCid)));
 342     hMCuncorrYieldPt[MCid]->GetYaxis()->SetTitle("1/N_{Jets} dN/dP_{T} (GeV/c)^{-1}");
 343     hMCuncorrYieldPt[MCid]->SetStats(kFALSE);
 344
 345     if (hPIDdata->GetAxis(kPidZ)) { // For backward compatibility
 346       hMCuncorrYieldZ[MCid] = hPIDdata->Projection(kPidZ, "e");
 347       hMCuncorrYieldZ[MCid]->SetName(Form("hMCuncorrYieldsPrimZ_%s", AliPID::ParticleShortName(MCid)));
 348       hMCuncorrYieldZ[MCid]->SetTitle(Form("MC truth uncorrected yield, %s", AliPID::ParticleName(MCid)));
 349       hMCuncorrYieldZ[MCid]->GetYaxis()->SetTitle("1/N_{Jets} dN/dz");
 350       hMCuncorrYieldZ[MCid]->SetStats(kFALSE);
 351     }
 352
 353     if (hPIDdata->GetAxis(kPidXi)) {
 354       hMCuncorrYieldXi[MCid] = hPIDdata->Projection(kPidXi, "e");
 355       hMCuncorrYieldXi[MCid]->SetName(Form("hMCuncorrYieldsPrimXi_%s", AliPID::ParticleShortName(MCid)));
 356       hMCuncorrYieldXi[MCid]->SetTitle(Form("MC truth uncorrected yield, %s", AliPID::ParticleName(MCid)));
 357       hMCuncorrYieldXi[MCid]->GetYaxis()->SetTitle("1/N_{Jets} dN/d#xi");
 358       hMCuncorrYieldXi[MCid]->SetStats(kFALSE);
 359     }
 360
 361     hPIDdata->GetAxis(kPidMCpid)->SetRange(0, -1);
 362   }
 363
 364   // Save results to file
 365   TString chargeString = "";
 366   if (chargeMode == kPosCharge)
 367     chargeString = "_posCharge";
 368   else if (chargeMode == kNegCharge)
 369     chargeString = "_negCharge";
 370
 371   TString saveFileName = pathNameData;
 372   saveFileName.Replace(0, pathNameData.Last('/') + 1, "");
 373
 374   TString savePath = pathNameData;
 375   savePath.ReplaceAll(Form("/%s", saveFileName.Data()), "");
 376
 377   saveFileName.Prepend("output_extractedGenPrimYields_");
 378   TString centralityString = restrictCentralityData ? Form("_centrality_%.0f_%.0f.root", actualLowerCentralityData,
 379                                                            actualUpperCentralityData)
 380                                                     : "_centrality_all";
 381   TString jetPtString = restrictJetPtAxis ? Form("_jetPt_%.0f_%.0f.root", actualLowerJetPt, actualUpperJetPt)
 382                                           : "";
 383   saveFileName.ReplaceAll(".root", Form("%s%s%s.root", centralityString.Data(), jetPtString.Data(), chargeString.Data()));
 384
 385   TString saveFilePathName = Form("%s/%s", savePath.Data(), saveFileName.Data());
 386   TFile* saveFile = TFile::Open(saveFilePathName.Data(), "RECREATE");
 387
 388   if (!saveFile) {
 389     printf("Failed to save results to file \"%s\"!\n", saveFilePathName.Data());
 390     return -1;
 391   }
 392
 393   saveFile->cd();
 394
 395
 396
 397   // Calculate the fractions and do the normalisation
 398
 399   const Double_t nJetsGenData = hNjetsGenData ? hNjetsGenData->Integral(lowerCentralityBinLimitData, upperCentralityBinLimitData,
 400                                                                         lowerJetPtBinLimitData, upperJetPtBinLimitData) : 1.;
 401   const Double_t nJetsRecData = hNjetsRecData ? hNjetsRecData->Integral(lowerCentralityBinLimitData, upperCentralityBinLimitData,
 402                                                                         lowerJetPtBinLimitData, upperJetPtBinLimitData) : 1.;
 403
 404   for (Int_t type = 0; type < kNtypes; type++) {
 405     TH1D** hMCgenPrimYield;
 406     TH1D** hMCuncorrYield;
 407
 408     if (type == kTrackPt) {
 409       hMCgenPrimYield = &hMCgenPrimYieldPt[0];
 410       hMCuncorrYield = &hMCuncorrYieldPt[0];
 411     }
 412     else if (type == kZ) {
 413       hMCgenPrimYield = &hMCgenPrimYieldZ[0];
 414       hMCuncorrYield = &hMCuncorrYieldZ[0];
 415     }
 416     else if (type == kXi) {
 417       hMCgenPrimYield = &hMCgenPrimYieldXi[0];
 418       hMCuncorrYield = &hMCuncorrYieldXi[0];
 419     }
 420     else
 421       continue;
 422
 423     // Gen yields
 424     if (hMCgenPrimYield[0]) {
 425       for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
 426         normaliseHist(hMCgenPrimYield[species], (nJetsGenData > 0) ? 1. / nJetsGenData : 0.);
 427         hMCgenPrimYield[species]->SetStats(kFALSE);
 428
 429         hMCgenPrimYield[species]->SetTitle(Form("%s (MC)", AliPID::ParticleLatexName(species)));
 430         hMCgenPrimYield[species]->SetMarkerStyle(24);
 431         hMCgenPrimYield[species]->SetLineColor(getLineColorAliPID(species));
 432         hMCgenPrimYield[species]->SetMarkerColor(getLineColorAliPID(species));
 433       }
 434
 435       TH1D* hMCgenPrimYieldTotal = 0x0;
 436       TH1D* hMCgenPrimFraction[AliPID::kSPECIES];
 437       for (Int_t i = 0; i < AliPID::kSPECIES; i++)
 438         hMCgenPrimFraction[i] = 0x0;
 439
 440       hMCgenPrimYieldTotal = new TH1D(*hMCgenPrimYield[0]);
 441       hMCgenPrimYieldTotal->SetLineColor(kBlack);
 442       hMCgenPrimYieldTotal->SetMarkerColor(kBlack);
 443       hMCgenPrimYieldTotal->SetMarkerStyle(24);
 444       hMCgenPrimYieldTotal->SetName(Form("hMCgenPrimYieldTotal%s", titles[type].Data()));
 445       hMCgenPrimYieldTotal->SetTitle("Total (MC)");
 446       //hMCgenPrimYieldTotal->SetTitle("Total generated primary yield (MC truth)");
 447
 448       for (Int_t i = 1; i < AliPID::kSPECIES; i++)
 449         hMCgenPrimYieldTotal->Add(hMCgenPrimYield[i], 1.);
 450
 451       // Calculate the MC fractions
 452       for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
 453         hMCgenPrimYield[species]->SetLineColor(getLineColorAliPID(species));
 454         hMCgenPrimYield[species]->SetMarkerColor(getLineColorAliPID(species));
 455         hMCgenPrimYield[species]->SetMarkerStyle(24);
 456
 457         hMCgenPrimYield[species]->SetTitle(Form("%s (MC)", AliPID::ParticleLatexName(species)));
 458
 459         hMCgenPrimFraction[species] = new TH1D(*hMCgenPrimYield[species]);
 460         TString oldName = hMCgenPrimYield[species]->GetName();
 461         TString newName = oldName.ReplaceAll("Yield", "Fraction");
 462         hMCgenPrimFraction[species]->SetName(newName.Data());
 463         hMCgenPrimFraction[species]->GetYaxis()->SetTitle("Particle Fraction");
 464
 465         // Binomial error as for efficiencies (numerator and denominator are statistically not independent) for correct error calculation
 466         // (numerator is a "selected" subset of the denominator).
 467         hMCgenPrimFraction[species]->Divide(hMCgenPrimFraction[species], hMCgenPrimYieldTotal, 1., 1., "B");
 468         hMCgenPrimFraction[species]->GetYaxis()->SetRangeUser(0., 1.);
 469       }
 470
 471       if (hMCgenPrimYieldTotal) {
 472         hMCgenPrimYieldTotal->GetYaxis()->SetTitleOffset(1.4);
 473         hMCgenPrimYieldTotal->GetXaxis()->SetMoreLogLabels(kTRUE);
 474         hMCgenPrimYieldTotal->GetXaxis()->SetNoExponent(kTRUE);
 475       }
 476
 477       TCanvas* cGenYield = new TCanvas(Form("cGenYield_%s", titles[type].Data()), "Generated Yield", 0, 300, 900, 900);
 478       cGenYield->SetLogx(1);
 479       cGenYield->SetLogy(1);
 480
 481       hMCgenPrimYieldTotal->Draw("E1");
 482
 483       for (Int_t i = 0; i < AliPID::kSPECIES; i++)
 484         hMCgenPrimYield[i]->Draw("E1 same");
 485
 486       TLegend* legTemp = cGenYield->BuildLegend();//0.25, 0.16, 0.65, 0.51);
 487       legTemp->SetFillColor(kWhite);
 488
 489       ClearTitleFromHistoInCanvas(cGenYield);
 490
 491
 492       TCanvas* cGenFractions = new TCanvas(Form("cGenFractions_%s", titles[type].Data()), "Generated Fractions", 0, 300, 900, 900);
 493       cGenFractions->SetLogx(1);
 494
 495       hMCgenPrimFraction[0]->Draw("E1");
 496
 497       for (Int_t i = 1; i < AliPID::kSPECIES; i++)
 498         hMCgenPrimFraction[i]->Draw("E1 same");
 499
 500       cGenFractions->BuildLegend()->SetFillColor(kWhite);
 501
 502       ClearTitleFromHistoInCanvas(cGenFractions);
 503
 504       // Save results
 505       saveFile->cd();
 506
 507       for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
 508         if (hMCgenPrimYield[i])
 509           hMCgenPrimYield[i]->Write();
 510       }
 511
 512       if (hMCgenPrimYieldTotal)
 513           hMCgenPrimYieldTotal->Write();
 514
 515       for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
 516         if (hMCgenPrimFraction[i])
 517           hMCgenPrimFraction[i]->Write();
 518       }
 519
 520       if (cGenYield)
 521         cGenYield->Write();
 522
 523       if (cGenFractions)
 524         cGenFractions->Write();
 525     }
 526
 527
 528     // Uncorr gen yields
 529     if (hMCuncorrYield[0]) {
 530       for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
 531         normaliseHist(hMCuncorrYield[species], (nJetsGenData > 0) ? 1. / nJetsGenData : 0.);
 532         hMCuncorrYield[species]->SetStats(kFALSE);
 533
 534         hMCuncorrYield[species]->SetTitle(Form("%s (MC)", AliPID::ParticleLatexName(species)));
 535         hMCuncorrYield[species]->SetMarkerStyle(24);
 536         hMCuncorrYield[species]->SetLineColor(getLineColorAliPID(species));
 537         hMCuncorrYield[species]->SetMarkerColor(getLineColorAliPID(species));
 538       }
 539
 540       TH1D* hMCuncorrYieldTotal = 0x0;
 541       TH1D* hMCuncorrFraction[AliPID::kSPECIES];
 542       for (Int_t i = 0; i < AliPID::kSPECIES; i++)
 543         hMCuncorrFraction[i] = 0x0;
 544
 545       hMCuncorrYieldTotal = new TH1D(*hMCuncorrYield[0]);
 546       hMCuncorrYieldTotal->SetLineColor(kBlack);
 547       hMCuncorrYieldTotal->SetMarkerColor(kBlack);
 548       hMCuncorrYieldTotal->SetMarkerStyle(24);
 549       hMCuncorrYieldTotal->SetName(Form("hMCuncorrYieldTotal%s", titles[type].Data()));
 550       hMCuncorrYieldTotal->SetTitle("Total (MC)");
 551       //hMCuncorrYieldTotal->SetTitle("Total uncorrected yield (MC truth)");
 552
 553       for (Int_t i = 1; i < AliPID::kSPECIES; i++)
 554         hMCuncorrYieldTotal->Add(hMCuncorrYield[i], 1.);
 555
 556       // Calculate the MC fractions
 557       for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
 558         hMCuncorrYield[species]->SetLineColor(getLineColorAliPID(species));
 559         hMCuncorrYield[species]->SetMarkerColor(getLineColorAliPID(species));
 560         hMCuncorrYield[species]->SetMarkerStyle(24);
 561
 562         hMCuncorrYield[species]->SetTitle(Form("%s (MC)", AliPID::ParticleLatexName(species)));
 563
 564         hMCuncorrFraction[species] = new TH1D(*hMCuncorrYield[species]);
 565         TString oldName = hMCuncorrYield[species]->GetName();
 566         TString newName = oldName.ReplaceAll("Yield", "Fraction");
 567         hMCuncorrFraction[species]->SetName(newName.Data());
 568         hMCuncorrFraction[species]->GetYaxis()->SetTitle("Particle Fraction");
 569
 570         // Binomial error as for efficiencies (numerator and denominator are statistically not independent) for correct error calculation
 571         // (numerator is a "selected" subset of the denominator).
 572         hMCuncorrFraction[species]->Divide(hMCuncorrFraction[species], hMCuncorrYieldTotal, 1., 1., "B");
 573         hMCuncorrFraction[species]->GetYaxis()->SetRangeUser(0., 1.);
 574       }
 575
 576       if (hMCuncorrYieldTotal) {
 577         hMCuncorrYieldTotal->GetYaxis()->SetTitleOffset(1.4);
 578         hMCuncorrYieldTotal->GetXaxis()->SetMoreLogLabels(kTRUE);
 579         hMCuncorrYieldTotal->GetXaxis()->SetNoExponent(kTRUE);
 580       }
 581
 582       TCanvas* cGenUncorrYield = new TCanvas(Form("cGenUncorrYield_%s", titles[type].Data()), "Generated Uncorrected Yield",
 583                                              0, 300, 900, 900);
 584       cGenUncorrYield->SetLogx(1);
 585       cGenUncorrYield->SetLogy(1);
 586
 587       hMCuncorrYieldTotal->Draw("E1");
 588
 589       for (Int_t i = 0; i < AliPID::kSPECIES; i++)
 590         hMCuncorrYield[i]->Draw("E1 same");
 591
 592       TLegend* legTemp = cGenUncorrYield->BuildLegend();//0.25, 0.16, 0.65, 0.51);
 593       legTemp->SetFillColor(kWhite);
 594
 595       ClearTitleFromHistoInCanvas(cGenUncorrYield);
 596
 597
 598       TCanvas* cGenUncorrFractions = new TCanvas(Form("cGenUncorrFractions_%s", titles[type].Data()), "Generated Uncorrected Fractions",
 599                                            0, 300, 900, 900);
 600       cGenUncorrFractions->SetLogx(1);
 601
 602       hMCuncorrFraction[0]->Draw("E1");
 603
 604       for (Int_t i = 1; i < AliPID::kSPECIES; i++)
 605         hMCuncorrFraction[i]->Draw("E1 same");
 606
 607       cGenUncorrFractions->BuildLegend()->SetFillColor(kWhite);
 608
 609       ClearTitleFromHistoInCanvas(cGenUncorrFractions);
 610
 611       // Save results
 612       saveFile->cd();
 613
 614       for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
 615         if (hMCuncorrYield[i])
 616           hMCuncorrYield[i]->Write();
 617       }
 618
 619       if (hMCuncorrYieldTotal)
 620           hMCuncorrYieldTotal->Write();
 621
 622       for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
 623         if (hMCuncorrFraction[i])
 624           hMCuncorrFraction[i]->Write();
 625       }
 626
 627       if (cGenUncorrYield)
 628         cGenUncorrYield->Write();
 629
 630       if (cGenUncorrFractions)
 631         cGenUncorrFractions->Write();
 632     }
 633   }
 634
 635   saveFile->cd();
 636
 637   TNamed* settings = new TNamed(
 638       Form("Settings: Data file \"%s\", lowerCentrality %.3f, upperCentrality %.3f, lowerJetPt %.1f, upperJetPt %.1f\n",
 639            pathNameData.Data(), lowerCentrality, upperCentrality, lowerJetPt, upperJetPt), "");
 640   settings->Write();
 641
 642   saveFile->Close();
 643
 644
 645   return 0;
 646 }