PWGJE/macros/PID/calcEfficiency.C

   1 #include "TH1D.h"
   2 #include "TH2D.h"
   3 #include "TCanvas.h"
   4 #include "TFile.h"
   5 #include "TGraphAsymmErrors.h"
   6 #include "TLegend.h"
   7
   8 #include "AliCFContainer.h"
   9 #include "AliCFEffGrid.h"
  10 #include "AliCFDataGrid.h"
  11 #include "AliPID.h"
  12
  13 #include <iostream>
  14
  15 #include "THnSparseDefinitions.h"
  16
  17 enum type { kTrackPt = 0, kZ = 1, kXi = 2, kNtypes = 3 };
  18
  19
  20 const Int_t nPtBinsType2 = 44;
  21 const Double_t binsPtType2[nPtBinsType2+1] = {0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45,
  22           0.5, 0.55, 0.6, 0.65, 0.7, 0.8, 0.9, 1.0, 1.2, 1.4,
  23           1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.4, 3.8,
  24           4.5, 5.5, 6.5, 8.0, 10.0, 12.0, 14.0, 16.0, 20.0, 24.0,
  25           28.0, 32.0, 36.0, 45.0, 50.0 };
  26
  27
  28 //___________________________________________________________________
  29 const Double_t* GetBins(Int_t type, Int_t& nBins)
  30 {
  31   if (type == 2) {
  32     nBins = nPtBinsType2;
  33
  34     return binsPtType2;
  35   }
  36
  37   return 0x0;
  38 }
  39
  40
  41 //___________________________________________________________________
  42 Double_t trackingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
  43 {
  44   return (1. - 0.129758 * TMath::Exp(-pTmc * 0.679612));
  45 }
  46
  47
  48 //___________________________________________________________________
  49 Double_t trackingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
  50 {
  51   return TMath::Min((0.972865 + 0.0117093 * pTmc), 1.);
  52 }
  53
  54
  55 //___________________________________________________________________
  56 Bool_t geantFlukaCorrection(AliCFContainer* data, Int_t genStepToDownscale)
  57 {
  58   // Issue: GEANT/FLUKA correction factor is for MC_pT.
  59   // Finally the effeciency should be DIVIDED by this correction factor.
  60   // To include resolution effects, it is therefore the best way to just
  61   // multiply the generated step with the correction factor.
  62
  63   if (!data) {
  64     printf("No CFContainer for GEANT/FLUKA correction!\n");
  65     return kFALSE;
  66   }
  67
  68   const Int_t iPt     = data->GetVar("P_{T} (GeV/c)");
  69   const Int_t iMCid   = data->GetVar("MC ID");
  70   const Int_t iCharge = data->GetVar("Charge (e_{0})");
  71
  72   if (iPt < 0 || iMCid < 0 || iCharge < 0) {
  73     printf("Data axis for GEANT/FLUKA correction not found!\n");
  74     return kFALSE;
  75   }
  76
  77   if (!data->GetGrid(genStepToDownscale)) {
  78     printf("Step for downscaling (GEANT/FLUKA) not found!\n");
  79     return kFALSE;
  80   }
  81
  82   const Int_t nDim = data->GetNVar();
  83   Int_t coord[nDim];
  84   Double_t binCenterCoord[nDim];
  85
  86   Long64_t nBinsGrid = data->GetGrid(genStepToDownscale)->GetGrid()->GetNbins();
  87
  88   for (Long64_t iBin = 0; iBin < nBinsGrid; iBin++) {
  89     Double_t binContent = data->GetGrid(genStepToDownscale)->GetGrid()->GetBinContent(iBin, coord);
  90     Double_t binError  = data->GetGrid(genStepToDownscale)->GetGrid()->GetBinError(iBin);
  91
  92     for (Int_t iDim = 0; iDim < nDim; iDim++)
  93       binCenterCoord[iDim] = data->GetBinCenter(iDim, coord[iDim]);
  94
  95     if (binCenterCoord[iCharge] < 0) {
  96       Double_t corrFactor = 1.;
  97
  98       if (binCenterCoord[iMCid] - 0.5 == AliPID::kProton)
  99         corrFactor = trackingPtGeantFlukaCorrectionPrMinus(binCenterCoord[iPt]);
 100       else if (binCenterCoord[iMCid] - 0.5 == AliPID::kKaon)
 101         corrFactor = trackingPtGeantFlukaCorrectionKaMinus(binCenterCoord[iPt]);
 102       else
 103         continue;
 104
 105       data->GetGrid(genStepToDownscale)->GetGrid()->SetBinContent(iBin, binContent * corrFactor);
 106       data->GetGrid(genStepToDownscale)->GetGrid()->SetBinError(iBin, binError * corrFactor);
 107     }
 108   }
 109
 110   return kTRUE;
 111 }
 112
 113
 114
 115 //___________________________________________________________________
 116 void setupHist(TH1* h, TString histName, TString histTitle, TString xAxisTitle, TString yAxisTitle, Int_t color)
 117 {
 118   if (histName != "")
 119     h->SetName(histName.Data());
 120   h->SetTitle(histTitle.Data());
 121
 122   if (xAxisTitle != "")
 123     h->GetXaxis()->SetTitle(xAxisTitle.Data());
 124   if (yAxisTitle != "")
 125     h->GetYaxis()->SetTitle(yAxisTitle.Data());
 126
 127   h->SetMarkerStyle(24);
 128   h->SetLineColor(color);
 129   h->SetMarkerColor(color);
 130
 131   h->SetStats(kFALSE);
 132 }
 133
 134
 135 //____________________________________________________________________________________________________________________
 136 void undoJetNormalisationHist2D(TH2* hData, TH2* hNumJets, const Int_t lowerCentralityBinLimit, const Int_t upperCentralityBinLimit)
 137 {
 138   // Undo normalisation to 1/numJets. NOTE: jetPt binning of hData and hNumJets assumed to be the same!
 139
 140   if (!hData || !hNumJets)
 141     return;
 142
 143   for (Int_t binJetPt = 0; binJetPt <= hData->GetNbinsY() + 1; binJetPt++) {
 144     const Double_t numJets = hNumJets->Integral(lowerCentralityBinLimit, upperCentralityBinLimit, binJetPt, binJetPt);
 145     Bool_t noJets = numJets < 1e-13;
 146
 147     for (Int_t binObs = 0; binObs <= hData->GetNbinsX() + 1; binObs++) {
 148       if (noJets)
 149         continue;
 150
 151       hData->SetBinContent(binObs, binJetPt, hData->GetBinContent(binObs, binJetPt) * numJets);
 152       hData->SetBinError(binObs, binJetPt, hData->GetBinError(binObs, binJetPt) * numJets);
 153     }
 154   }
 155 }
 156
 157
 158 //___________________________________________________________________
 159 void normaliseHist(TH1* h, Double_t scale = 1.0)
 160 {
 161   if (h->GetSumw2N() <= 0)
 162     h->Sumw2();
 163
 164   h->Scale(scale);
 165
 166   for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
 167     Double_t normFactor = h->GetBinWidth(i);
 168     h->SetBinContent(i, h->GetBinContent(i) / normFactor);
 169     h->SetBinError(i, h->GetBinError(i) / normFactor);
 170   }
 171 }
 172
 173
 174 //___________________________________________________________________
 175 void multiplyHistsDifferentBinning(TH1* h1, const TH1* h2, Double_t c1, Double_t c2, Bool_t ignoreErrorOfSecondHist = kFALSE)
 176 {
 177    Int_t nbinsx = h1->GetNbinsX();
 178    Int_t nbinsy = h1->GetNbinsY();
 179    Int_t nbinsz = h1->GetNbinsZ();
 180
 181    if (h1->GetDimension() < 2)
 182      nbinsy -= 1;
 183    if (h1->GetDimension() < 3)
 184      nbinsz -= 1;
 185
 186    if (h1->GetSumw2N() == 0)
 187      h1->Sumw2();
 188
 189    Int_t bin, bin2, binx, biny, binz;
 190    Double_t b1,b2,w,d1,d2;
 191    d1 = c1*c1;
 192    d2 = c2*c2;
 193
 194    // Loop over bins (including underflows/overflows)
 195    for (binz = 0; binz <= nbinsz + 1; binz++) {
 196       for (biny = 0; biny <= nbinsy + 1; biny++) {
 197          for (binx = 0; binx <= nbinsx + 1;binx++) {
 198             bin = binx + (nbinsx + 2) * (biny + (nbinsy + 2) * binz);
 199             bin2 = h2->FindFixBin(h1->GetXaxis()->GetBinCenter(binx), h1->GetYaxis()->GetBinCenter(biny),
 200                                   h1->GetZaxis()->GetBinCenter(binz));
 201             b1  = h1->GetBinContent(bin);
 202             b2  = h2->GetBinContent(bin2);
 203             w   = (c1*b1)*(c2*b2);
 204             h1->SetBinContent(bin, w);
 205             Double_t e1 = h1->GetBinError(bin);
 206             Double_t e2 = ignoreErrorOfSecondHist ? 0. : h2->GetBinError(bin2);
 207             h1->SetBinError(bin, TMath::Sqrt(d1*d2*(e1*e1*b2*b2 + e2*e2*b1*b1)));
 208          }
 209       }
 210    }
 211 }
 212
 213
 214 //___________________________________________________________________
 215 void divideHistsDifferentBinning(TH1* h1, const TH1* h2, Double_t c1, Double_t c2, Bool_t ignoreErrorOfSecondHist = kFALSE)
 216 {
 217    Int_t nbinsx = h1->GetNbinsX();
 218    Int_t nbinsy = h1->GetNbinsY();
 219    Int_t nbinsz = h1->GetNbinsZ();
 220
 221    if (h1->GetDimension() < 2)
 222      nbinsy -= 1;
 223    if (h1->GetDimension() < 3)
 224      nbinsz -= 1;
 225
 226    if (h1->GetSumw2N() == 0)
 227      h1->Sumw2();
 228
 229    Int_t bin, bin2, binx, biny, binz;
 230    Double_t b1,b2,w,d1,d2;
 231    d1 = c1*c1;
 232    d2 = c2*c2;
 233
 234    // Loop over bins (including underflows/overflows)
 235    for (binz = 0; binz <= nbinsz + 1; binz++) {
 236       for (biny = 0; biny <= nbinsy + 1; biny++) {
 237          for (binx = 0; binx <= nbinsx + 1;binx++) {
 238             bin = binx + (nbinsx + 2) * (biny + (nbinsy + 2) * binz);
 239             bin2 = h2->FindFixBin(h1->GetXaxis()->GetBinCenter(binx), h1->GetYaxis()->GetBinCenter(biny),
 240                                   h1->GetZaxis()->GetBinCenter(binz));
 241             b1  = h1->GetBinContent(bin);
 242             b2  = h2->GetBinContent(bin2);
 243             if (b2)
 244               w = (c1*b1)/(c2*b2);
 245             else
 246               w = 0;
 247             h1->SetBinContent(bin, w);
 248
 249             if (!b2) {
 250               h1->SetBinError(bin, 0);
 251               continue;
 252             }
 253
 254             Double_t b22 = b2*b2*d2;
 255             Double_t e1 = h1->GetBinError(bin);
 256             Double_t e2 = ignoreErrorOfSecondHist ? 0. : h2->GetBinError(bin2);
 257             h1->SetBinError(bin, TMath::Sqrt(d1*d2*(e1*e1*b2*b2 + e2*e2*b1*b1)/(b22*b22)));
 258          }
 259       }
 260    }
 261 }
 262
 263
 264 //___________________________________________________________________
 265 // Efficiency for inclusive spectra vs. pT (or also jets, but without using the generation)
 266 // E.g. a 'calcEfficiency.C+("finalCuts/MC_pp/7TeV/LHC10f6a/corrected/finalisedSplines/analytical/efficiency_noClCut_preliminary/bhess_PID_Jets_Inclusive_PureGauss_efficiency.root", "finalCuts/pp/7TeV/LHC10e.pass2/corrected/finalisedSplines/finalMapsAndTail/Jets/noCutOn_ncl_or_liav/outputSystematicsTotal_SummedSystematicErrors__2013_10_21.root", kTRUE, kFALSE, 0, -2, -2, -1, -1, 8, 2)' -b -q
 267 Int_t calcEfficiency(TString pathNameEfficiency, TString pathNameData, Bool_t correctGeantFluka, Bool_t scaleStrangeness,
 268                      Int_t chargeMode /*kNegCharge = -1, kAllCharged = 0, kPosCharge = 1*/,
 269                      Double_t lowerCentrality /*= -2*/, Double_t upperCentrality /*= -2*/,
 270                      Double_t lowerJetPt /*= -1*/ , Double_t upperJetPt/* = -1*/,
 271                      Double_t constantCorrectionAbovePtThreshold,
 272                      Int_t rebinEfficiency)
 273 {
 274   TString pathData = pathNameData;
 275   pathData.Replace(pathData.Last('/'), pathData.Length(), "");
 276
 277   TFile* fileEff = TFile::Open(pathNameEfficiency.Data());
 278   if (!fileEff) {
 279     printf("Failed to open efficiency file \"%s\"\n", pathNameEfficiency.Data());
 280     return -1;
 281   }
 282
 283   AliCFContainer* data = (AliCFContainer*)(fileEff->Get("containerEff"));
 284   if (!data) {
 285     printf("Failed to load efficiency container!\n");
 286     return -1;
 287   }
 288
 289   const Int_t iMCid   = data->GetVar("MC ID");
 290   const Int_t iPt     = data->GetVar("P_{T} (GeV/c)");
 291   const Int_t iEta    = data->GetVar("#eta");
 292   const Int_t iCharge = data->GetVar("Charge (e_{0})");
 293   const Int_t iMult   = data->GetVar("Centrality Percentile");
 294   const Int_t iJetPt  = data->GetVar("P_{T}^{jet} (GeV/c)");
 295   //const Int_t iZ      = data->GetVar("z = P_{T}^{track} / P_{T}^{jet}");
 296   //const Int_t iXi     = data->GetVar("#xi = ln(P_{T}^{jet} / P_{T}^{track})");
 297
 298   TFile* fileData = TFile::Open(pathNameData.Data());
 299   if (!fileData) {
 300     printf("Failed to open data file \"%s\"\n", pathNameData.Data());
 301     return -1;
 302   }
 303
 304   TH1D* hYield[AliPID::kSPECIES] = { 0x0, };
 305   TH1D* hYieldSysError[AliPID::kSPECIES] = { 0x0, };
 306   TGraphAsymmErrors* gYieldSysError[AliPID::kSPECIES] = { 0x0, };
 307   TH1D* hMCgenPrimYield[AliPID::kSPECIES] = { 0x0, };
 308   Int_t numMCgenPrimYieldHistsFound = 0;
 309
 310   for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
 311     TString speciesName = AliPID::ParticleName(species);
 312     TString firstLetter = speciesName(0);
 313     firstLetter.ToUpper();
 314     speciesName.Replace(0, 1, firstLetter.Data());
 315     TString histName = Form("hYield%ss", speciesName.Data());
 316     hYield[species] = (TH1D*)fileData->Get(histName.Data());
 317     if (!hYield[species]) {
 318       printf("Failed to load hist \"%s\"\n", histName.Data());
 319       return -1;
 320     }
 321     hYield[species]->SetFillStyle(0);
 322
 323     TString graphName = Form("systematicErrorYields_%s", AliPID::ParticleName(species));
 324     gYieldSysError[species] = (TGraphAsymmErrors*)fileData->Get(graphName.Data());
 325
 326     // In case of MC also retrieve the MC truth generated yields
 327     TString histNameMCgenYields = Form("hMCgenYieldsPrimSpecies_%s", AliPID::ParticleShortName(species));
 328     hMCgenPrimYield[species] = (TH1D*)fileData->Get(histNameMCgenYields.Data());
 329     if (hMCgenPrimYield[species])
 330       numMCgenPrimYieldHistsFound++;
 331   }
 332
 333   if (numMCgenPrimYieldHistsFound > 0 && numMCgenPrimYieldHistsFound != AliPID::kSPECIES) {
 334     printf("Error: Unable to retrieve all MC generated prim yield histos! Got %d.\n", numMCgenPrimYieldHistsFound);
 335     return -1;
 336   }
 337
 338   TCanvas* cFractions = (TCanvas*)fileData->Get("cFractionsWithTotalSystematicError");
 339   if (!cFractions)
 340     cFractions = (TCanvas*)fileData->Get("cFractions");
 341
 342
 343   // Convert graphs with systematic errors into histograms (assume symmetric error, which is currently true)
 344   for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
 345     if (gYieldSysError[species]) {
 346       hYieldSysError[species] = new TH1D(*hYield[species]);
 347       hYieldSysError[species]->SetName(Form("%s_sysError", hYieldSysError[species]->GetName()));
 348       hYieldSysError[species]->SetFillStyle(0);
 349
 350       for (Int_t binX = 1; binX <= hYieldSysError[species]->GetNbinsX(); binX++) {
 351         hYieldSysError[species]->SetBinContent(binX, gYieldSysError[species]->GetY()[binX - 1]);
 352         hYieldSysError[species]->SetBinError(binX, gYieldSysError[species]->GetErrorY(binX - 1));
 353       }
 354     }
 355   }
 356
 357   // Take pT binning from pion yield (binning for all species the same) and create a new AliCFContainer with this new binning for pT
 358   TH1D hDummy(*hYield[AliPID::kPion]);
 359   hDummy.SetName("hDummy");
 360   TAxis* axis = 0x0;
 361
 362   if (rebinEfficiency > 1) {
 363     Int_t nBinsNew = 0;
 364     const Double_t* newBins = GetBins(rebinEfficiency, nBinsNew);
 365
 366     hDummy.SetBins(nBinsNew, newBins);
 367
 368     axis = hDummy.GetXaxis();
 369
 370     //axis = hDummy.Rebin(rebinEfficiency, "", 0)->GetXaxis();
 371   }
 372   else
 373     axis = hDummy.GetXaxis();
 374
 375   const TArrayD* binsPtCurrent = axis->GetXbins();
 376   TArrayD* binsPtNew = new TArrayD(*binsPtCurrent);
 377
 378   const Int_t nEffDims = data->GetNVar();
 379   Int_t nEffBins[nEffDims];
 380
 381   for (Int_t iDim = 0; iDim < nEffDims; iDim++) {
 382     if (iDim == iPt)
 383       nEffBins[iDim] = axis->GetNbins();
 384     else
 385       nEffBins[iDim] = data->GetNBins(iDim);
 386   }
 387
 388
 389   // Just make one large pT bin above some threshold, if desired
 390   if (binsPtNew->fN != 0 && constantCorrectionAbovePtThreshold > 0) {
 391     for (Int_t iBin = 0; iBin < nEffBins[iPt]; iBin++) {
 392       // Find the first bin edged really larger than the threshold.
 393       // If the bin edge before equals the threshold, just set the
 394       // current bin edge to the right end of the spectrum -> Done.
 395       // If the bin edge before is different, set the bin edge to the
 396       // threshold
 397       if (binsPtNew->fArray[iBin] > constantCorrectionAbovePtThreshold) {
 398         if (binsPtNew->fArray[iBin - 1] == constantCorrectionAbovePtThreshold) {
 399           binsPtNew->fArray[iBin] = binsPtNew->fArray[nEffBins[iPt]];
 400           nEffBins[iPt] = iBin;
 401           break;
 402         }
 403         else {
 404           binsPtNew->fArray[iBin] = constantCorrectionAbovePtThreshold;
 405         }
 406       }
 407     }
 408   }
 409
 410
 411   AliCFContainer *dataRebinned = new AliCFContainer(Form("%s_rebinned", data->GetName()), Form("%s (rebinned)", data->GetTitle()),
 412                                                     data->GetNStep(), nEffDims, nEffBins);
 413
 414   for (Int_t iDim = 0; iDim < nEffDims; iDim++) {
 415     dataRebinned->SetVarTitle(iDim, data->GetVarTitle(iDim));
 416
 417     if (iDim == iPt) {
 418       if (binsPtNew->fN == 0)
 419         dataRebinned->SetBinLimits(iDim, axis->GetXmin(), axis->GetXmax());
 420       else
 421         dataRebinned->SetBinLimits(iDim, binsPtNew->fArray);
 422     }
 423     else {
 424       dataRebinned->SetBinLimits(iDim, data->GetBinLimits(iDim));
 425     }
 426   }
 427
 428   for (Int_t iStep = 0; iStep < data->GetNStep(); iStep++)
 429     dataRebinned->SetStepTitle(iStep, data->GetStepTitle(iStep));
 430
 431   Int_t coord[nEffDims];
 432   Double_t binCenterCoord[nEffDims];
 433
 434   // Fill content from old grid into the new grid with proper binning
 435   for (Int_t iStep = 0; iStep < data->GetNStep(); iStep++) {
 436     Long64_t nBinsGrid = data->GetGrid(iStep)->GetGrid()->GetNbins();
 437
 438     for (Long64_t iBin = 0; iBin < nBinsGrid; iBin++) {
 439       Double_t binContent = data->GetGrid(iStep)->GetGrid()->GetBinContent(iBin, coord);
 440       Double_t binError2  = data->GetGrid(iStep)->GetGrid()->GetBinError2(iBin);
 441
 442       for (Int_t iDim = 0; iDim < nEffDims; iDim++) {
 443         binCenterCoord[iDim] = data->GetBinCenter(iDim, coord[iDim]);
 444       }
 445
 446       Long64_t iBinRebinned = dataRebinned->GetGrid(iStep)->GetGrid()->GetBin(binCenterCoord);
 447       dataRebinned->GetGrid(iStep)->GetGrid()->AddBinContent(iBinRebinned, binContent);
 448       dataRebinned->GetGrid(iStep)->GetGrid()->AddBinError2(iBinRebinned, binError2);
 449     }
 450   }
 451
 452
 453   // If desired, restrict centrality axis
 454   Int_t lowerCentralityBinLimit = -1;
 455   Int_t upperCentralityBinLimit = -2; // Integral(lowerCentBinLimit, uppCentBinLimit) will not be restricted if these values are kept
 456   Bool_t restrictCentralityAxis = kFALSE;
 457   Double_t actualLowerCentrality = -1.;
 458   Double_t actualUpperCentrality = -1.;
 459
 460   if (lowerCentrality >= -1 && upperCentrality >= -1) {
 461     // Add subtract a very small number to avoid problems with values right on the border between to bins
 462     lowerCentralityBinLimit = dataRebinned->GetAxis(iMult, 0)->FindBin(lowerCentrality + 0.001);
 463     upperCentralityBinLimit = dataRebinned->GetAxis(iMult, 0)->FindBin(upperCentrality - 0.001);
 464
 465     // Check if the values look reasonable
 466     if (lowerCentralityBinLimit <= upperCentralityBinLimit && lowerCentralityBinLimit >= 1
 467         && upperCentralityBinLimit <= dataRebinned->GetAxis(iMult, 0)->GetNbins()) {
 468       actualLowerCentrality = dataRebinned->GetAxis(iMult, 0)->GetBinLowEdge(lowerCentralityBinLimit);
 469       actualUpperCentrality = dataRebinned->GetAxis(iMult, 0)->GetBinUpEdge(upperCentralityBinLimit);
 470
 471       restrictCentralityAxis = kTRUE;
 472     }
 473     else {
 474       std::cout << std::endl;
 475       std::cout << "Requested centrality range out of limits or upper and lower limit are switched!" << std::endl;
 476       return -1;
 477     }
 478   }
 479
 480   std::cout << "centrality: ";
 481   if (restrictCentralityAxis) {
 482     std::cout << actualLowerCentrality << " - " << actualUpperCentrality << std::endl;
 483     dataRebinned->SetRangeUser(iMult, lowerCentralityBinLimit, upperCentralityBinLimit, kTRUE);
 484   }
 485   else {
 486     std::cout << "All" << std::endl;
 487   }
 488
 489
 490
 491   // If desired, restrict jetPt axis
 492   Int_t lowerJetPtBinLimit = -1;
 493   Int_t upperJetPtBinLimit = -1;
 494   Bool_t restrictJetPtAxis = kFALSE;
 495   Double_t actualLowerJetPt = -1.;
 496   Double_t actualUpperJetPt = -1.;
 497
 498   if (lowerJetPt >= 0 && upperJetPt >= 0) {
 499     // Add subtract a very small number to avoid problems with values right on the border between to bins
 500     lowerJetPtBinLimit = dataRebinned->GetAxis(iJetPt, 0)->FindBin(lowerJetPt + 0.001);
 501     upperJetPtBinLimit = dataRebinned->GetAxis(iJetPt, 0)->FindBin(upperJetPt - 0.001);
 502
 503     // Check if the values look reasonable
 504     if (lowerJetPtBinLimit <= upperJetPtBinLimit && lowerJetPtBinLimit >= 1 &&
 505         upperJetPtBinLimit <= dataRebinned->GetAxis(iJetPt, 0)->GetNbins()) {
 506       actualLowerJetPt = dataRebinned->GetAxis(iJetPt, 0)->GetBinLowEdge(lowerJetPtBinLimit);
 507       actualUpperJetPt = dataRebinned->GetAxis(iJetPt, 0)->GetBinUpEdge(upperJetPtBinLimit);
 508
 509       restrictJetPtAxis = kTRUE;
 510     }
 511     else {
 512       std::cout << std::endl;
 513       std::cout << "Requested jet pT range out of limits or upper and lower limit are switched!" << std::endl;
 514       return -1;
 515     }
 516   }
 517
 518   std::cout << "jet pT: ";
 519   if (restrictJetPtAxis) {
 520     std::cout << actualLowerJetPt << " - " << actualUpperJetPt << std::endl;
 521     dataRebinned->SetRangeUser(iJetPt, lowerJetPtBinLimit, upperJetPtBinLimit, kTRUE);
 522   }
 523   else {
 524     std::cout << "All" << std::endl;
 525   }
 526
 527   // If desired, restrict charge axis
 528   std::cout << "Charge selection (efficiency): ";
 529   if (chargeMode == kAllCharged)
 530     std::cout << "All charged particles" << std::endl;
 531   else if (chargeMode == kNegCharge)
 532     std::cout << "Negative particles only" << std::endl;
 533   else if (chargeMode == kPosCharge)
 534     std::cout << "Positive particles only" << std::endl;
 535   else {
 536     std::cout << "Unknown -> ERROR" << std::endl;
 537     return -1;
 538   }
 539
 540   const Bool_t restrictCharge = (chargeMode != kAllCharged);
 541
 542   Int_t lowerChargeBinLimit = -1;
 543   Int_t upperChargeBinLimit = -2;
 544   Double_t actualLowerCharge = -999;
 545   Double_t actualUpperCharge = -999;
 546
 547   if (restrictCharge) {
 548     // Add subtract a very small number to avoid problems with values right on the border between to bins
 549     if (chargeMode == kNegCharge) {
 550       lowerChargeBinLimit = dataRebinned->GetAxis(iCharge, 0)->FindBin(-1. + 0.001);
 551       upperChargeBinLimit = dataRebinned->GetAxis(iCharge, 0)->FindBin(0. - 0.001);
 552     }
 553     else if (chargeMode == kPosCharge) {
 554       lowerChargeBinLimit = dataRebinned->GetAxis(iCharge, 0)->FindBin(0. + 0.001);
 555       upperChargeBinLimit = dataRebinned->GetAxis(iCharge, 0)->FindBin(1. - 0.001);
 556     }
 557
 558     // Check if the values look reasonable
 559     if (lowerChargeBinLimit <= upperChargeBinLimit && lowerChargeBinLimit >= 1
 560         && upperChargeBinLimit <= dataRebinned->GetAxis(iCharge, 0)->GetNbins()) {
 561       actualLowerCharge = dataRebinned->GetAxis(iCharge, 0)->GetBinLowEdge(lowerChargeBinLimit);
 562       actualUpperCharge = dataRebinned->GetAxis(iCharge, 0)->GetBinUpEdge(upperChargeBinLimit);
 563
 564       std::cout << "Charge range (efficiency): " << actualLowerCharge << " - " << actualUpperCharge << std::endl;
 565     }
 566     else {
 567       std::cout << std::endl;
 568       std::cout << "Requested charge range (efficiency) out of limits or upper and lower limit are switched!" << std::endl;
 569       return -1;
 570     }
 571
 572     dataRebinned->SetRangeUser(iCharge, lowerChargeBinLimit, upperChargeBinLimit, kTRUE);
 573     data->SetRangeUser(iCharge, lowerChargeBinLimit, upperChargeBinLimit, kTRUE);
 574   }
 575
 576   std::cout << std::endl;
 577
 578
 579   // If jet axis is restricted (i.e. jet input), also need num jet histos for proper normalisation.
 580   // Load numJet histos from bhess_PID*.root file related to efficiency file.
 581   TH2D* hNjetsGen = 0x0;
 582   TH2D* hNjetsRec = 0x0;
 583
 584   if (restrictJetPtAxis) {
 585     TString pathNameDataMC = pathNameEfficiency;
 586     pathNameDataMC.ReplaceAll("_efficiency", "");
 587
 588     TFile* fDataMC = TFile::Open(pathNameDataMC.Data());
 589     if (!fDataMC)  {
 590       std::cout << std::endl;
 591       std::cout << "Failed to open file \"" << pathNameData.Data() << "\" to obtain num of rec/gen jets!" << std::endl;
 592
 593       return -1;
 594     }
 595
 596     TString listName = pathNameDataMC;
 597     listName.Replace(0, listName.Last('/') + 1, "");
 598     listName.ReplaceAll(".root", "");
 599
 600     TObjArray* histList = (TObjArray*)(fDataMC->Get(listName.Data()));
 601     if (!histList) {
 602       std::cout << std::endl;
 603       std::cout << "Failed to load list \"" << listName.Data() << "\" to obtain num of rec/gen jets!" << std::endl;
 604       return -1;
 605     }
 606
 607     hNjetsGen = (TH2D*)histList->FindObject("fh2FFJetPtGen");
 608     hNjetsRec = (TH2D*)histList->FindObject("fh2FFJetPtRec");
 609
 610     if (!hNjetsRec || !hNjetsGen) {
 611       std::cout << "Failed to load number of jets histos!" << std::endl;
 612
 613
 614       // For backward compatibility (TODO REMOVE IN FUTURE): Load info from fixed AnalysisResults file (might be wrong, if other
 615       // period is considered; also: No multiplicity information)
 616       TString pathEfficiency = pathNameEfficiency;
 617       pathEfficiency.Replace(pathEfficiency.Last('/'), pathEfficiency.Length(), "");
 618       TString pathBackward = Form("%s/AnalysisResults.root", pathEfficiency.Data());
 619       TFile* fBackward = TFile::Open(pathBackward.Data());
 620
 621       TString dirDataInFile = "";
 622       TDirectory* dirData = fBackward ? (TDirectory*)fBackward->Get(fBackward->GetListOfKeys()->At(0)->GetName()) : 0x0;
 623
 624       TList* list = dirData ? (TList*)dirData->Get(dirData->GetListOfKeys()->At(0)->GetName()) : 0x0;
 625
 626       TH1D* hFFJetPtRec = list ? (TH1D*)list->FindObject("fh1FFJetPtRecCuts") : 0x0;
 627       TH1D* hFFJetPtGen = list ? (TH1D*)list->FindObject("fh1FFJetPtGen") : 0x0;
 628
 629       if (hFFJetPtRec && hFFJetPtGen) {
 630         printf("***WARNING: For backward compatibility, using file \"%s\" to get number of jets. BUT: Might be wrong period and has no mult info!***\n",
 631           pathBackward.Data());
 632
 633         hNjetsRec = new TH2D("fh2FFJetPtRec", "", 1, -1, 1, dataRebinned->GetAxis(iJetPt, 0)->GetNbins(),
 634                             dataRebinned->GetAxis(iJetPt, 0)->GetXbins()->GetArray());
 635
 636         for (Int_t iJet = 1; iJet <= hNjetsRec->GetNbinsY(); iJet++) {
 637           Int_t lowerBin = hFFJetPtRec->FindBin(hNjetsRec->GetYaxis()->GetBinLowEdge(iJet) + 1e-3);
 638           Int_t upperBin = hFFJetPtRec->FindBin(hNjetsRec->GetYaxis()->GetBinUpEdge(iJet) - 1e-3);
 639           hNjetsRec->SetBinContent(1, iJet, hFFJetPtRec->Integral(lowerBin, upperBin));
 640         }
 641
 642         hNjetsGen = new TH2D("fh2FFJetPtGen", "", 1, -1, 1,  dataRebinned->GetAxis(iJetPt, 0)->GetNbins(),
 643                             dataRebinned->GetAxis(iJetPt, 0)->GetXbins()->GetArray());
 644
 645         for (Int_t iJet = 1; iJet <= hNjetsGen->GetNbinsY(); iJet++) {
 646           Int_t lowerBin = hFFJetPtGen->FindBin(hNjetsGen->GetYaxis()->GetBinLowEdge(iJet) + 1e-3);
 647           Int_t upperBin = hFFJetPtGen->FindBin(hNjetsGen->GetYaxis()->GetBinUpEdge(iJet) - 1e-3);
 648           hNjetsGen->SetBinContent(1, iJet, hFFJetPtGen->Integral(lowerBin, upperBin));
 649         }
 650       }
 651
 652       if (!hNjetsRec || ! hNjetsGen)
 653         return -1;
 654     }
 655   }
 656
 657   // For normalisation to number of jets
 658   // NOTE: These numbers are for the efficiency only! The data will be normalised to its own number!!!
 659   const Double_t nJetsGen = hNjetsGen ? hNjetsGen->Integral(lowerCentralityBinLimit, upperCentralityBinLimit, lowerJetPtBinLimit,
 660                                                             upperJetPtBinLimit) : 1.;
 661   const Double_t nJetsRec = hNjetsRec ? hNjetsRec->Integral(lowerCentralityBinLimit, upperCentralityBinLimit, lowerJetPtBinLimit,
 662                                                             upperJetPtBinLimit) : 1.;
 663
 664   // Secondary correction
 665   AliCFEffGrid* sec = new AliCFEffGrid("sec", "Secondary Contamination", *dataRebinned);
 666   AliCFEffGrid* secStrangeScale = new AliCFEffGrid("secStrangeScale", "Secondary Contamination with Strangeness Scaling",
 667                                                    *dataRebinned);
 668
 669   // Either one can take kStepRecWithGenCutsMeasuredObs or, what I prefer, one can take
 670   // kStepRecWithRecCutsMeasuredObsPrimaries => The difference is only the eta cut, which is on the rec level
 671   // in the latter case, i.e. one corrects for eta resolution (although the effect should be very small)
 672   // => TESTED: There is NO difference (only 1 bin vs. pt shows a deviation by one entry), so this cut has,
 673   // in principle, more or less no effect
 674
 675   // For testing with data set w/o strangeness secStrangeScale->CalculateEfficiency(kStepRecWithRecCutsMeasuredObsPrimaries, kStepRecWithRecCutsMeasuredObs);
 676   secStrangeScale->CalculateEfficiency(kStepRecWithRecCutsMeasuredObsPrimaries, kStepRecWithRecCutsMeasuredObsStrangenessScaled);
 677   sec->CalculateEfficiency(kStepRecWithRecCutsMeasuredObsPrimaries, kStepRecWithRecCutsMeasuredObs);
 678
 679   // QA plots for secondaries
 680   TCanvas* cSec = new TCanvas("cSec", "Secondary Contamination", 100, 10, 1200, 800);
 681   cSec->Divide(2, 1);
 682   cSec->GetPad(1)->SetLogx(kTRUE);
 683   cSec->cd(1);
 684   TH1* hSecPt = sec->Project(iPt);
 685   hSecPt->SetName("hSecPt");
 686   hSecPt->SetStats(kFALSE);
 687   hSecPt->GetXaxis()->SetRangeUser(0.15, 50.);
 688   hSecPt->GetXaxis()->SetMoreLogLabels(kTRUE);
 689   hSecPt->GetXaxis()->SetNoExponent(kTRUE);
 690   hSecPt->GetYaxis()->SetTitle("Primary Fraction");
 691   hSecPt->Draw("E1");
 692   cSec->cd(2);
 693   TH1* hEtaSec = sec->Project(iEta);
 694   hEtaSec->SetName("hEtaSec");
 695   hEtaSec->SetStats(kFALSE);
 696   hEtaSec->GetYaxis()->SetTitle("Primary Fraction");
 697   hEtaSec->Draw("E1");
 698   TH2D* hSecID2Pt = (TH2D*)sec->Project(iPt, iMCid);
 699   hSecID2Pt->SetName("hSecID2Pt");
 700
 701   // Get the secondary contamination vs. pT for each species
 702   TH1D* hSec[AliPID::kSPECIES];
 703   for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
 704     hSec[species] = hSecID2Pt->ProjectionX(Form("hSec_%s", AliPID::ParticleShortName(species)), species + 1, species + 1, "e");
 705     hSec[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species)));
 706     hSec[species]->SetLineColor(hYield[species]->GetLineColor());
 707     hSec[species]->SetMarkerColor(hYield[species]->GetLineColor());
 708     hSec[species]->SetLineWidth(2.);
 709     hSec[species]->GetXaxis()->SetRangeUser(0.15, 50);
 710     hSec[species]->GetYaxis()->SetRangeUser(0., 1.01);
 711     hSec[species]->GetXaxis()->SetMoreLogLabels(kTRUE);
 712     hSec[species]->GetXaxis()->SetNoExponent(kTRUE);
 713     hSec[species]->SetStats(kFALSE);
 714     hSec[species]->GetYaxis()->SetTitle("Primary Fraction");
 715   }
 716
 717   TCanvas* cSec2 = new TCanvas("cSec2", "Primary fraction for different species", 100, 10, 1200, 800);
 718   cSec2->SetGridx(1);
 719   cSec2->SetGridy(1);
 720   cSec2->SetLogx(1);
 721
 722   hSec[0]->Draw("E1");
 723
 724   for (Int_t i = 1; i < AliPID::kSPECIES; i++) {
 725     hSec[i]->Draw("E1 same");
 726   }
 727   cSec2->BuildLegend()->SetFillColor(kWhite);
 728
 729   ClearTitleFromHistoInCanvas(cSec2);
 730
 731   // QA plots for secondaries with strangeness scaling
 732   TCanvas* cSecSS = new TCanvas("cSecSS", "Secondary Contamination (strangeness scaled)", 100, 10, 1200, 800);
 733   cSecSS->Divide(2, 1);
 734   cSecSS->GetPad(1)->SetLogx(kTRUE);
 735   cSecSS->cd(1);
 736   TH1* hSecSSPt = secStrangeScale->Project(iPt);
 737   hSecSSPt->SetName("hSecSSPt");
 738   hSecSSPt->SetStats(kFALSE);
 739   hSecSSPt->GetXaxis()->SetRangeUser(0.15, 50.);
 740   hSecSSPt->GetXaxis()->SetMoreLogLabels(kTRUE);
 741   hSecSSPt->GetXaxis()->SetNoExponent(kTRUE);
 742   hSecSSPt->GetYaxis()->SetTitle("Primary Fraction");
 743   hSecSSPt->Draw("E1");
 744   cSecSS->cd(2);
 745   TH1* hEtaSecSS = secStrangeScale->Project(iEta);
 746   hEtaSecSS->SetName("hEtaSecSS");
 747   hEtaSecSS->SetStats(kFALSE);
 748   hEtaSecSS->GetYaxis()->SetTitle("Primary Fraction");
 749   hEtaSecSS->Draw("E1");
 750   TH2D* hSecSSID2Pt = (TH2D*)secStrangeScale->Project(iPt, iMCid);
 751   hSecSSID2Pt->SetName("hSecSSID2Pt");
 752
 753   // Get the secondary contamination vs. pT for each species
 754   TH1D* hSecSS[AliPID::kSPECIES];
 755   for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
 756     hSecSS[species] = hSecSSID2Pt->ProjectionX(Form("hSecSS_%s", AliPID::ParticleShortName(species)), species + 1, species + 1, "e");
 757     hSecSS[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species)));
 758     hSecSS[species]->SetLineColor(hYield[species]->GetLineColor());
 759     hSecSS[species]->SetMarkerColor(hYield[species]->GetLineColor());
 760     hSecSS[species]->SetLineWidth(2.);
 761     hSecSS[species]->GetXaxis()->SetRangeUser(0.15, 50);
 762     hSecSS[species]->GetYaxis()->SetRangeUser(0., 1.01);
 763     hSecSS[species]->GetXaxis()->SetMoreLogLabels(kTRUE);
 764     hSecSS[species]->GetXaxis()->SetNoExponent(kTRUE);
 765     hSecSS[species]->SetStats(kFALSE);
 766     hSecSS[species]->GetYaxis()->SetTitle("Primary Fraction");
 767   }
 768
 769   TCanvas* cSecSS2 = new TCanvas("cSecSS2", "Primary fraction for different species (strangeness scaled)", 100, 10, 1200, 800);
 770   cSecSS2->SetGridx(1);
 771   cSecSS2->SetGridy(1);
 772   cSecSS2->SetLogx(1);
 773
 774   hSecSS[0]->Draw("E1");
 775
 776   for (Int_t i = 1; i < AliPID::kSPECIES; i++) {
 777     hSecSS[i]->Draw("E1 same");
 778   }
 779   cSecSS2->BuildLegend()->SetFillColor(kWhite);
 780
 781   ClearTitleFromHistoInCanvas(cSecSS2);
 782
 783
 784   // Secondary correction for to-pi-ratios
 785   TH1D* hSecToPiRatio[AliPID::kSPECIES] = { 0x0, };
 786   TH1D* hSecSSToPiRatio[AliPID::kSPECIES] = { 0x0, };
 787
 788   for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
 789     if (species == AliPID::kPion)
 790       continue; // Do not consider pion-to-pion ratio
 791
 792     hSecToPiRatio[species] = new TH1D(*hSec[species]);
 793     hSecToPiRatio[species]->Reset();
 794     hSecToPiRatio[species]->SetName(Form("hSecSSToPionRatio_%s", AliPID::ParticleShortName(species)));
 795     hSecToPiRatio[species]->SetTitle(Form("%s/#pi", AliPID::ParticleLatexName(species)));
 796     hSecToPiRatio[species]->SetLineColor(getLineColorAliPID(species));
 797     hSecToPiRatio[species]->SetMarkerColor(getLineColorAliPID(species));
 798     hSecToPiRatio[species]->SetMarkerStyle(24);
 799     hSecToPiRatio[species]->SetLineWidth(2.);
 800     hSecToPiRatio[species]->SetStats(kFALSE);
 801     hSecToPiRatio[species]->GetYaxis()->SetTitle("Primary Fraction of Ratio");
 802
 803     hSecSSToPiRatio[species] = new TH1D(*hSecToPiRatio[species]);
 804     hSecSSToPiRatio[species]->SetName(Form("hSecSSToPionRatio_%s", AliPID::ParticleShortName(species)));
 805
 806     // Samples for different species are independent, so just divide correction factors
 807     hSecToPiRatio[species]->Divide(hSec[species], hSec[AliPID::kPion], 1., 1., "");
 808     hSecToPiRatio[species]->GetYaxis()->SetRangeUser(0., 1.1);
 809
 810     hSecSSToPiRatio[species]->Divide(hSecSS[species], hSecSS[AliPID::kPion], 1., 1., "");
 811     hSecSSToPiRatio[species]->GetYaxis()->SetRangeUser(0., 1.1);
 812   }
 813
 814   TCanvas* cSecToPiRatio = new TCanvas("cSecToPiRatio", "Primary fraction of to-#pi-ratio for different species", 100, 10, 1200, 800);
 815   cSecToPiRatio->SetGridx(1);
 816   cSecToPiRatio->SetGridy(1);
 817   cSecToPiRatio->SetLogx(1);
 818
 819   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
 820     if (i == AliPID::kPion)
 821       continue;
 822
 823     hSecToPiRatio[i]->Draw(Form("E1%s", i == 0 ? "" : " same"));
 824   }
 825
 826   cSecToPiRatio->BuildLegend()->SetFillColor(kWhite);
 827
 828   ClearTitleFromHistoInCanvas(cSecToPiRatio);
 829
 830
 831
 832   TCanvas* cSecSSToPiRatio = new TCanvas("cSecSSToPiRatio",
 833                                          "Primary fraction of to-#pi-ratio for different species (strangeness scaled)",
 834                                          100, 10, 1200, 800);
 835   cSecSSToPiRatio->SetGridx(1);
 836   cSecSSToPiRatio->SetGridy(1);
 837   cSecSSToPiRatio->SetLogx(1);
 838
 839   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
 840     if (i == AliPID::kPion)
 841       continue;
 842
 843     hSecSSToPiRatio[i]->Draw(Form("E1%s", i == 0 ? "" : " same"));
 844   }
 845
 846   cSecSSToPiRatio->BuildLegend()->SetFillColor(kWhite);
 847
 848   ClearTitleFromHistoInCanvas(cSecSSToPiRatio);
 849
 850
 851
 852
 853   // If desired, apply the GEANT/FLUKA correction first
 854   // NOTE: This will change dataRebinned! If anything else should also be calculated, it must be done before.
 855   // Otherwise, the GEANT/FLUKA correction factor for the efficiency will also affect it!
 856   const Int_t genStepEff = kStepGenWithGenCuts;
 857   if (correctGeantFluka) {
 858     printf("Applying GEANT/FLUKA correction...\n");
 859     if (!geantFlukaCorrection(dataRebinned, genStepEff)) {
 860       printf("GEANT/FLUKA correction could not be applied!\n");
 861       return kFALSE;
 862     }
 863   }
 864
 865   // Construct the efficiency grid from the data container
 866   AliCFEffGrid* eff = new AliCFEffGrid("eff", "Efficiency x Acceptance x pT Resolution", *dataRebinned);
 867
 868   // Either one can take kStepRecWithGenCutsMeasuredObs or, what I prefer, one can take
 869   // kStepRecWithRecCutsMeasuredObsPrimaries => The difference is only the eta cut, which is on the rec level
 870   // in the latter case, i.e. one corrects for eta resolution (although the effect should be very small)
 871   eff->CalculateEfficiency(kStepRecWithRecCutsMeasuredObsPrimaries, genStepEff);
 872
 873   // If the jet axis is restricted (i.e. jet input), scale with the corresponding number of jets.
 874   // Note: Since this is supposed to be a real scaling, set the error of the scale factor to zero
 875   // (second element of factor array)
 876   if (restrictJetPtAxis) {
 877     Double_t factor_Numerator[2] = { nJetsRec > 0 ? 1. / nJetsRec : 0., 0.  };
 878     Double_t factor_Denominator[2] = { nJetsGen > 0 ? 1. / nJetsGen : 0., 0.  };
 879     eff->GetNum()->Scale(factor_Numerator);
 880     eff->GetDen()->Scale(factor_Denominator);
 881   }
 882
 883   //The efficiency along pt and vertex, and 2-D projection
 884   TCanvas* cEff = new TCanvas("cEff", "Efficiency x Acceptance x pT Resolution", 100, 10, 1200, 800);
 885   cEff->Divide(2, 1);
 886   cEff->GetPad(1)->SetLogx(kTRUE);
 887   cEff->cd(1);
 888   TH1* hEffPt = eff->Project(iPt); //the efficiency vs pt
 889   hEffPt->SetStats(kFALSE);
 890   hEffPt->GetXaxis()->SetRangeUser(0.15, 50.);
 891   hEffPt->GetXaxis()->SetMoreLogLabels(kTRUE);
 892   hEffPt->GetXaxis()->SetNoExponent(kTRUE);
 893   hEffPt->GetYaxis()->SetTitle("Efficiency x Acceptance x pT Resolution");
 894   hEffPt->Draw("E1");
 895   cEff->cd(2);
 896   TH1* hEffEta = eff->Project(iEta); //the efficiency vs eta
 897   hEffEta->SetStats(kFALSE);
 898   hEffEta->GetYaxis()->SetTitle("Efficiency x Acceptance x pT Resolution");
 899   hEffEta->GetYaxis()->SetTitleSize(0.05);
 900   hEffEta->Draw("E1");
 901   TH2D* hEffID2Pt = (TH2D*)eff->Project(iPt, iMCid);
 902
 903   // Get the efficiencies vs. pT for each species
 904   TH1D* hEfficiency[AliPID::kSPECIES];
 905   for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
 906     hEfficiency[species] = hEffID2Pt->ProjectionX(Form("hEfficiency_%s", AliPID::ParticleShortName(species)), species + 1, species + 1, "e");
 907     hEfficiency[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species)));
 908     hEfficiency[species]->SetLineColor(hYield[species]->GetLineColor());
 909     hEfficiency[species]->SetMarkerColor(hYield[species]->GetLineColor());
 910     hEfficiency[species]->SetLineWidth(2.);
 911     hEfficiency[species]->GetXaxis()->SetRangeUser(0.15, 50);
 912     hEfficiency[species]->GetYaxis()->SetRangeUser(0., 1.01);
 913     hEfficiency[species]->GetXaxis()->SetMoreLogLabels(kTRUE);
 914     hEfficiency[species]->GetXaxis()->SetNoExponent(kTRUE);
 915     hEfficiency[species]->SetStats(kFALSE);
 916     hEfficiency[species]->GetYaxis()->SetTitle("Efficiency x Acceptance x pT Resolution");
 917     hEfficiency[species]->GetYaxis()->SetTitleSize(0.05);
 918   }
 919
 920   TCanvas* cEff2 = new TCanvas("cEff2", "Efficiency x Acceptance x pT Resolution for different species", 0, 300, 900, 900);
 921   cEff2->SetGridx(1);
 922   cEff2->SetGridy(1);
 923   cEff2->SetLogx(1);
 924
 925   hEfficiency[0]->Draw("E1");
 926
 927   for (Int_t i = 1; i < AliPID::kSPECIES; i++) {
 928     hEfficiency[i]->Draw("E1 same");
 929   }
 930   cEff2->BuildLegend()->SetFillColor(kWhite);
 931
 932   ClearTitleFromHistoInCanvas(cEff2);
 933
 934
 935
 936
 937   // Efficiency correction for to-pi-ratios
 938   TH1D* hEfficiencyToPiRatio[AliPID::kSPECIES] = { 0x0, };
 939
 940   for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
 941     if (species == AliPID::kPion)
 942       continue; // Do not consider pion-to-pion ratio
 943
 944     hEfficiencyToPiRatio[species] = new TH1D(*hEfficiency[species]);
 945     hEfficiencyToPiRatio[species]->Reset();
 946     hEfficiencyToPiRatio[species]->SetName(Form("hEfficiencyToPionRatio_%s", AliPID::ParticleShortName(species)));
 947     hEfficiencyToPiRatio[species]->SetTitle(Form("%s/#pi", AliPID::ParticleLatexName(species)));
 948
 949     // Samples for different species are independent, so just divide correction factors
 950     hEfficiencyToPiRatio[species]->Divide(hEfficiency[species], hEfficiency[AliPID::kPion], 1., 1., "");
 951     hEfficiencyToPiRatio[species]->GetYaxis()->SetRangeUser(0., 2.0);
 952   }
 953
 954   TCanvas* cEffToPiRatio = new TCanvas("cEffToPiRatio", "Efficiency x Acceptance x pT Resolution of to-#pi-ratio for different species",
 955                                        100, 10, 1200, 800);
 956   cEffToPiRatio->SetGridx(1);
 957   cEffToPiRatio->SetGridy(1);
 958   cEffToPiRatio->SetLogx(1);
 959
 960   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
 961     if (i == AliPID::kPion)
 962       continue;
 963
 964     hEfficiencyToPiRatio[i]->Draw(Form("E1%s", i == 0 ? "" : " same"));
 965   }
 966
 967   cEffToPiRatio->BuildLegend()->SetFillColor(kWhite);
 968
 969   ClearTitleFromHistoInCanvas(cEffToPiRatio);
 970
 971
 972   // Correct the yields with the efficiencies and primary fractions
 973   TH1D* hYieldCorrected[AliPID::kSPECIES] = { 0x0, };
 974   TH1D* hYieldCorrectedSysError[AliPID::kSPECIES] = { 0x0, };
 975   for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
 976     hYieldCorrected[species] = new TH1D(*hYield[species]);
 977     hYieldCorrected[species]->SetName(Form("%s_corrected", hYield[species]->GetName()));
 978     hYieldCorrected[species]->SetTitle(Form("%s", hYield[species]->GetTitle()));
 979     //hYieldCorrected[species]->SetTitle(Form("%s, secondary and efficiency x acceptance x pT resolution corrected",
 980     //                                        hYield[species]->GetTitle()));
 981
 982     // Correction factor histos can have different binning than data histos (constant factor above some threshold)
 983     // -> Need special functions to multiply and divide such histos
 984     multiplyHistsDifferentBinning(hYieldCorrected[species], scaleStrangeness ? hSecSS[species] : hSec[species], 1., 1.);
 985     divideHistsDifferentBinning(hYieldCorrected[species], hEfficiency[species], 1., 1.);
 986
 987     if (hYieldSysError[species]) {
 988       hYieldCorrectedSysError[species] = new TH1D(*hYieldSysError[species]);
 989       hYieldCorrectedSysError[species]->SetName(Form("%s_corrected", hYieldSysError[species]->GetName()));
 990       hYieldCorrectedSysError[species]->SetTitle(Form("%s", hYieldSysError[species]->GetTitle()));
 991
 992       multiplyHistsDifferentBinning(hYieldCorrectedSysError[species], scaleStrangeness ? hSecSS[species] : hSec[species], 1., 1.,
 993                                     kTRUE);
 994       divideHistsDifferentBinning(hYieldCorrectedSysError[species], hEfficiency[species], 1., 1., kTRUE);
 995     }
 996   }
 997
 998   // Calculate the total corrected yield. The original total yield had no error (just the number of detected tracks in a pT bin),
 999   // but due to the correction there is some error for the total yield. Also the error of the fractions introduces uncertainties
1000   // for the yields of individual species
1001   TH1D* hYieldCorrectedTotal = new TH1D(*hYieldCorrected[0]);
1002   hYieldCorrectedTotal->SetLineColor(kBlack);
1003   hYieldCorrectedTotal->SetMarkerColor(kBlack);
1004   hYieldCorrectedTotal->SetName("hYieldCorrectedTotal");
1005   hYieldCorrectedTotal->SetTitle("Total");
1006   //hYieldCorrectedTotal->SetTitle("Total yield, secondary and efficiency x acceptance x pT resolution corrected");
1007
1008   for (Int_t i = 1; i < AliPID::kSPECIES; i++)
1009     hYieldCorrectedTotal->Add(hYieldCorrected[i], 1.);
1010
1011   // Calculate the corrected fractions
1012   TH1D* hFractionCorrected[AliPID::kSPECIES] = { 0x0, };
1013   TH1D* hFractionCorrectedSysError[AliPID::kSPECIES] = { 0x0, };
1014   for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
1015     hFractionCorrected[species] = new TH1D(*hYield[species]);
1016     TString oldName = hYield[species]->GetName();
1017     TString newName = oldName.ReplaceAll("Yield", "Fraction");
1018     TString oldTitle = hYield[species]->GetTitle();
1019     TString newTitle = oldTitle.ReplaceAll("Yield", "Fraction");
1020     hFractionCorrected[species]->SetName(Form("%s_corrected", newName.Data()));
1021     hFractionCorrected[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species)));
1022     hFractionCorrected[species]->GetYaxis()->SetTitle("Corrected Fraction");
1023     hFractionCorrected[species]->GetXaxis()->SetMoreLogLabels(kTRUE);
1024     hFractionCorrected[species]->GetXaxis()->SetNoExponent(kTRUE);
1025
1026     // Binomial error as for efficiencies (numerator and denominator are statistically not independent) for correct error calculation
1027     // (numerator is a "selected" subset of the denominator). It doesn't matter that the error of the histos is not just "sqrt(content)"
1028     // because the error formula also works for weighted histograms (which means that the error can be more or less anything)
1029     hFractionCorrected[species]->Divide(hYieldCorrected[species], hYieldCorrectedTotal, 1., 1., "B");
1030
1031
1032     //  The systematic errors just need to be scaled in the same way as the fractions.
1033     // So, just take the ratio to the uncorrected fraction and scale the sys. error accordingly
1034     // or, in this case, just divide by the same total yield as for yield -> fractions
1035     if (hYieldCorrectedSysError[species]) {
1036       hFractionCorrectedSysError[species] = new TH1D(*hFractionCorrected[species]);
1037       hFractionCorrectedSysError[species]->SetName(Form("%s_sysError", hFractionCorrected[species]->GetName()));
1038       hFractionCorrectedSysError[species]->SetTitle(Form("%s (sys. error)", hFractionCorrected[species]->GetTitle()));
1039
1040       for (Int_t binX = 1; binX <= hFractionCorrectedSysError[species]->GetNbinsX(); binX++) {
1041         const Double_t corrTotalYield = hYieldCorrectedTotal->GetBinContent(binX);
1042         const Double_t scaleFactor = corrTotalYield > 0 ? 1.0 / corrTotalYield : 1.;
1043         hFractionCorrectedSysError[species]->SetBinError(binX,   hYieldCorrectedSysError[species]->GetBinError(binX) * scaleFactor);
1044       }
1045     }
1046   }
1047
1048   // If MC is available, calculate the generated fractions
1049   TH1D* hMCgenPrimYieldTotal = 0x0;
1050   TH1D* hMCgenPrimFraction[AliPID::kSPECIES];
1051   for (Int_t i = 0; i < AliPID::kSPECIES; i++)
1052     hMCgenPrimFraction[i] = 0x0;
1053
1054   if (numMCgenPrimYieldHistsFound > 0) {
1055     hMCgenPrimYieldTotal = new TH1D(*hMCgenPrimYield[0]);
1056     hMCgenPrimYieldTotal->SetLineColor(kBlack);
1057     hMCgenPrimYieldTotal->SetMarkerColor(kBlack);
1058     hMCgenPrimYieldTotal->SetName("hMCgenPrimYieldTotal");
1059     hMCgenPrimYieldTotal->SetTitle("Total (MC truth)");
1060     //hMCgenPrimYieldTotal->SetTitle("Total generated primary yield (MC truth)");
1061
1062     for (Int_t i = 1; i < AliPID::kSPECIES; i++)
1063       hMCgenPrimYieldTotal->Add(hMCgenPrimYield[i], 1.);
1064
1065     // Calculate the MC fractions
1066     for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
1067       hMCgenPrimYield[species]->SetTitle(Form("%s (MC truth)", AliPID::ParticleLatexName(species)));
1068
1069       hMCgenPrimFraction[species] = new TH1D(*hMCgenPrimYield[species]);
1070       TString oldName = hMCgenPrimYield[species]->GetName();
1071       TString newName = oldName.ReplaceAll("Yield", "Fraction");
1072       TString oldTitle = hMCgenPrimYield[species]->GetTitle();
1073       TString newTitle = oldTitle.ReplaceAll("yield", "fraction");
1074       hMCgenPrimFraction[species]->SetName(newName.Data());
1075       hMCgenPrimFraction[species]->SetTitle(newTitle.Data());
1076
1077       // Binomial error as for efficiencies (numerator and denominator are statistically not independent) for correct error calculation
1078       // (numerator is a "selected" subset of the denominator).
1079       hMCgenPrimFraction[species]->Divide(hMCgenPrimFraction[species], hMCgenPrimYieldTotal, 1., 1., "B");
1080     }
1081   }
1082
1083   TCanvas* cCorrData = new TCanvas("cCorrData", "Corrected data", 0, 300, 900, 900);
1084   cCorrData->Divide(2, 1);//, 0., 0.01);
1085   cCorrData->GetPad(1)->SetLogx(1);
1086   cCorrData->GetPad(1)->SetLogy(1);
1087   cCorrData->GetPad(2)->SetLogx(1);
1088   cCorrData->GetPad(2)->SetLogy(1);
1089
1090   cCorrData->GetPad(1)->SetRightMargin(0.001);
1091   cCorrData->GetPad(2)->SetRightMargin(0.001);
1092
1093   cCorrData->GetPad(1)->SetLeftMargin(0.2);
1094   cCorrData->GetPad(2)->SetLeftMargin(0.2);
1095
1096   cCorrData->cd(1); // uncorrected
1097   hYield[AliPID::kPion]->GetYaxis()->SetTitleOffset(1.4);
1098   hYield[AliPID::kPion]->Draw("E1");
1099
1100   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1101     if (hYieldSysError[i])
1102       hYieldSysError[i]->Draw("E2 same");
1103
1104     if (i == AliPID::kPion) continue;
1105     hYield[i]->Draw("E1 same");
1106   }
1107
1108   ClearTitleFromHistoInCanvas(cCorrData, 1);
1109
1110   cCorrData->cd(2); // corrected
1111   hYieldCorrectedTotal->GetYaxis()->SetTitleOffset(1.4);
1112   hYieldCorrectedTotal->GetXaxis()->SetMoreLogLabels(kTRUE);
1113   hYieldCorrectedTotal->GetXaxis()->SetNoExponent(kTRUE);
1114   hYieldCorrectedTotal->GetYaxis()->SetRangeUser(hYieldCorrected[AliPID::kMuon]->GetBinContent(hYieldCorrected[AliPID::kMuon]->FindLastBinAbove(0.)) * 0.1,
1115                                                  hYieldCorrectedTotal->GetBinContent(hYieldCorrectedTotal->GetMaximumBin()) * 10.);
1116
1117   if (hMCgenPrimYieldTotal) {
1118     hMCgenPrimYieldTotal->GetYaxis()->SetTitleOffset(1.4);
1119     hMCgenPrimYieldTotal->GetXaxis()->SetMoreLogLabels(kTRUE);
1120     hMCgenPrimYieldTotal->GetXaxis()->SetNoExponent(kTRUE);
1121     hMCgenPrimYieldTotal->GetXaxis()->SetTitle(hYieldCorrectedTotal->GetXaxis()->GetTitle());
1122     hMCgenPrimYieldTotal->GetYaxis()->SetRangeUser(hYieldCorrected[AliPID::kMuon]->GetBinContent(hYieldCorrected[AliPID::kMuon]->FindLastBinAbove(0.)) * 0.1,
1123                                                    hYieldCorrectedTotal->GetBinContent(hYieldCorrectedTotal->GetMaximumBin()) * 10.);
1124   }
1125
1126   if (hMCgenPrimYieldTotal) {
1127     hMCgenPrimYieldTotal->Draw("E1");
1128     hYieldCorrectedTotal->Draw("E1 same");
1129   }
1130   else
1131     hYieldCorrectedTotal->Draw("E1");
1132
1133   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1134     if (hMCgenPrimYield[i])
1135       hMCgenPrimYield[i]->Draw("E1 same");
1136     hYieldCorrected[i]->Draw("E1 same");
1137   }
1138
1139   TLegend* legTemp = cCorrData->cd(2)->BuildLegend(0.25, 0.16, 0.65, 0.51);
1140   legTemp->SetNColumns(2);
1141   legTemp->SetFillColor(kWhite);
1142
1143   // Do not include in legend
1144   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1145     if (hYieldCorrectedSysError[i])
1146       hYieldCorrectedSysError[i]->Draw("E2 same");
1147   }
1148
1149   ClearTitleFromHistoInCanvas(cCorrData, 2);
1150
1151   TCanvas* cCorrYieldsRatio = 0x0;
1152
1153   TH1D* hYieldCorrectedRatioToMC[AliPID::kSPECIES];
1154   for (Int_t i = 0; i < AliPID::kSPECIES; i++)
1155     hYieldCorrectedRatioToMC[i] = 0x0;
1156
1157   TH1D* hYieldCorrectedTotalRatioToMC = 0x0;
1158
1159   if (numMCgenPrimYieldHistsFound > 0) {
1160     // Compare with MC truth
1161     for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
1162       hYieldCorrectedRatioToMC[species] = new TH1D(*hYieldCorrected[species]);
1163       hYieldCorrectedRatioToMC[species]->SetName(Form("hYieldCorrectedRatioToMC_%s", AliPID::ParticleShortName(species)));
1164       hYieldCorrectedRatioToMC[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species)));
1165       hYieldCorrectedRatioToMC[species]->GetYaxis()->SetTitle("Corrected Yield: Fit / MC Truth");
1166       hYieldCorrectedRatioToMC[species]->Divide(hMCgenPrimYield[species]);
1167     }
1168
1169     hYieldCorrectedTotalRatioToMC = new TH1D(*hYieldCorrectedTotal);
1170     hYieldCorrectedTotalRatioToMC->SetName("hYieldCorrectedTotalRatioToMC");
1171     hYieldCorrectedTotalRatioToMC->SetTitle("Total");
1172     hYieldCorrectedTotalRatioToMC->GetYaxis()->SetTitle("Corrected Yield: Fit / MC Truth");
1173     hYieldCorrectedTotalRatioToMC->Divide(hMCgenPrimYieldTotal);
1174
1175     cCorrYieldsRatio = new TCanvas("cCorrYieldsRatio", "Corrected Yields Comparison to MC", 0, 300, 900, 900);
1176     cCorrYieldsRatio->SetGridx(1);
1177     cCorrYieldsRatio->SetGridy(1);
1178     cCorrYieldsRatio->SetLogx(1);
1179
1180     hYieldCorrectedTotalRatioToMC->GetYaxis()->SetRangeUser(0.6, 1.6);
1181     hYieldCorrectedTotalRatioToMC->GetYaxis()->SetTitleOffset(0.85);
1182     hYieldCorrectedTotalRatioToMC->Draw("E1");
1183
1184     for (Int_t species = 0; species < AliPID::kSPECIES; species++)
1185       hYieldCorrectedRatioToMC[species]->Draw("E1 same");
1186
1187     cCorrYieldsRatio->BuildLegend()->SetFillColor(kWhite);
1188
1189     ClearTitleFromHistoInCanvas(cCorrYieldsRatio);
1190   }
1191
1192
1193   if (cFractions)
1194     cFractions->Draw();
1195
1196   TCanvas* cCorrFractions = new TCanvas("cCorrFractions", "Corrected particleFractions", 0, 300, 900, 900);
1197   cCorrFractions->SetLogx(1);
1198   hFractionCorrected[0]->GetYaxis()->SetRangeUser(0., 1.);
1199   hFractionCorrected[0]->Draw("E1");
1200   if (hMCgenPrimFraction[0])
1201     hMCgenPrimFraction[0]->Draw("E1 same");
1202
1203   for (Int_t i = 1; i < AliPID::kSPECIES; i++) {
1204     hFractionCorrected[i]->Draw("E1 same");
1205     if (hMCgenPrimFraction[i])
1206       hMCgenPrimFraction[i]->Draw("E1 same");
1207   }
1208
1209   cCorrFractions->BuildLegend()->SetFillColor(kWhite);
1210
1211   // Do not include in legend!!
1212   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1213     if (hFractionCorrectedSysError[i])
1214       hFractionCorrectedSysError[i]->Draw("E2 same");
1215   }
1216
1217
1218   ClearTitleFromHistoInCanvas(cCorrFractions);
1219
1220   // Save results to file
1221   TString chargeString = "";
1222   if (chargeMode == kPosCharge)
1223     chargeString = "_posCharge";
1224   else if (chargeMode == kNegCharge)
1225     chargeString = "_negCharge";
1226
1227   TString saveFileName = pathNameData;
1228   saveFileName.ReplaceAll(Form("%s/", pathData.Data()), "");
1229   saveFileName.Prepend("output_EfficiencyCorrection_");
1230   saveFileName.ReplaceAll(".root", Form("%s.root", chargeString.Data()));
1231
1232   TString saveFilePathName = Form("%s/%s", pathData.Data(), saveFileName.Data());
1233   TFile* saveFile = TFile::Open(saveFilePathName.Data(), "RECREATE");
1234   saveFile->cd();
1235
1236   if (cSec)
1237     cSec->Write();
1238
1239   if (cSecSS)
1240     cSecSS->Write();
1241
1242   if (cSec2)
1243     cSec2->Write();
1244
1245   if (cSecSS2)
1246     cSecSS2->Write();
1247
1248   if (cSecToPiRatio)
1249     cSecToPiRatio->Write();
1250
1251   if (cSecSSToPiRatio)
1252     cSecSSToPiRatio->Write();
1253
1254   if (cEff)
1255     cEff->Write();
1256
1257   if (cEff2)
1258     cEff2->Write();
1259
1260   if (cEffToPiRatio)
1261     cEffToPiRatio->Write();
1262
1263   if (cCorrData)
1264     cCorrData->Write();
1265
1266   if (cCorrYieldsRatio)
1267     cCorrYieldsRatio->Write();
1268
1269   if (cFractions)
1270     cFractions->Write();
1271
1272   if (cCorrFractions)
1273     cCorrFractions->Write();
1274
1275   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1276     if (hSec[i])
1277       hSec[i]->Write();
1278
1279     if (hSecSS[i])
1280       hSecSS[i]->Write();
1281   }
1282
1283   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1284     if (hSecToPiRatio[i])
1285       hSecToPiRatio[i]->Write();
1286
1287     if (hSecSSToPiRatio[i])
1288       hSecSSToPiRatio[i]->Write();
1289   }
1290
1291   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1292     if (hEfficiency[i])
1293       hEfficiency[i]->Write();
1294   }
1295
1296   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1297     if (hEfficiencyToPiRatio[i])
1298       hEfficiencyToPiRatio[i]->Write();
1299   }
1300
1301   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1302     if (hYield[i])
1303       hYield[i]->Write();
1304
1305     if (hYieldSysError[i])
1306       hYieldSysError[i]->Write();
1307   }
1308
1309   if (hYieldCorrectedTotal)
1310     hYieldCorrectedTotal->Write();
1311
1312   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1313     if (hYieldCorrected[i])
1314       hYieldCorrected[i]->Write();
1315
1316     if (hYieldCorrectedSysError[i])
1317       hYieldCorrectedSysError[i]->Write();
1318   }
1319
1320   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1321     if (hMCgenPrimYield[i])
1322       hMCgenPrimYield[i]->Write();
1323   }
1324
1325   if (hMCgenPrimYieldTotal)
1326       hMCgenPrimYieldTotal->Write();
1327
1328   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1329     if (hYieldCorrectedRatioToMC[i])
1330       hYieldCorrectedRatioToMC[i]->Write();
1331   }
1332
1333   if (hYieldCorrectedTotalRatioToMC)
1334       hYieldCorrectedTotalRatioToMC->Write();
1335
1336   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1337     if (hFractionCorrected[i])
1338       hFractionCorrected[i]->Write();
1339
1340     if (hFractionCorrectedSysError[i])
1341       hFractionCorrectedSysError[i]->Write();
1342   }
1343
1344   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1345     if (hMCgenPrimFraction[i])
1346       hMCgenPrimFraction[i]->Write();
1347   }
1348
1349   TNamed* settings = new TNamed(
1350       Form("Settings: Efficiency file \"%s\", Data file \"%s\", lowerCentrality %.3f, upperCentrality %.3f, lowerJetPt %.1f, upperJetPt %.1f, constantCorrectionAbovePtThreshold %.3f\n",
1351            pathNameEfficiency.Data(), pathNameData.Data(), lowerCentrality, upperCentrality, lowerJetPt,
1352            upperJetPt, constantCorrectionAbovePtThreshold), "");
1353   settings->Write();
1354
1355   saveFile->Close();
1356
1357   return 0;
1358 }
1359
1360
1361 //___________________________________________________________________
1362 // Efficiency for jet spectra (pT, z and xi)
1363 // E.g. a 'calcEfficiency.C+("finalCuts/MC_pp/7TeV/LHC11a1/corrected/allJetPtMergedWeighted/bhess_PID_Jets_PureGauss_efficiency.root", "finalCuts/pp/7TeV/LHC10e.pass2/corrected/finalisedSplines/finalMapsAndTail/Jets/noCutOn_ncl_or_liav/output_extractedFFs_bhess_PID_Jets_centrality_all_highN_rougherXi.root", kTRUE, kFALSE, 0, -2, -2, -2, -2, 40, 80, -100, -100, 2, 1, 2)' -b -q
1364 Int_t calcEfficiency(TString pathNameEfficiency, TString pathNameData, Bool_t correctGeantFluka, Bool_t scaleStrangeness,
1365                      Int_t chargeMode /*kNegCharge = -1, kAllCharged = 0, kPosCharge = 1*/,
1366                      Double_t lowerCentralityData /*= -2*/, Double_t upperCentralityData /*= -2*/,
1367                      Double_t lowerCentrality /*= -2*/, Double_t upperCentrality /*= -2*/,
1368                      Double_t lowerJetPt /*= -1*/ , Double_t upperJetPt/* = -1*/,
1369                      Double_t constantCorrectionAbovePtThreshold,
1370                      Double_t constantCorrectionAboveXiThreshold,
1371                      Int_t rebinEfficiencyPt,
1372                      Int_t rebinEfficiencyZ,
1373                      Int_t rebinEfficiencyXi)
1374 {
1375   TString titles[kNtypes];
1376   titles[kTrackPt] = "trackPt";
1377   titles[kZ] = "z";
1378   titles[kXi] = "xi";
1379
1380   TString pathData = pathNameData;
1381   pathData.Replace(pathData.Last('/'), pathData.Length(), "");
1382
1383   TFile* fileEff = TFile::Open(pathNameEfficiency.Data());
1384   if (!fileEff) {
1385     printf("Failed to open efficiency file \"%s\"\n", pathNameEfficiency.Data());
1386     return -1;
1387   }
1388
1389   AliCFContainer* data = (AliCFContainer*)(fileEff->Get("containerEff"));
1390   if (!data) {
1391     printf("Failed to load efficiency container!\n");
1392     return -1;
1393   }
1394
1395   const Int_t iMCid   = data->GetVar("MC ID");
1396   const Int_t iPt     = data->GetVar("P_{T} (GeV/c)");
1397   //const Int_t iEta    = data->GetVar("#eta");
1398   const Int_t iCharge = data->GetVar("Charge (e_{0})");
1399   const Int_t iMult   = data->GetVar("Centrality Percentile");
1400   const Int_t iJetPt  = data->GetVar("P_{T}^{jet} (GeV/c)");
1401   const Int_t iZ      = data->GetVar("z = P_{T}^{track} / P_{T}^{jet}");
1402   const Int_t iXi     = data->GetVar("#xi = ln(P_{T}^{jet} / P_{T}^{track})");
1403
1404   TH2D* hNjetsGen = 0x0;
1405   TH2D* hNjetsRec = 0x0;
1406
1407   TH2D* hNjetsGenData = 0x0;
1408   TH2D* hNjetsRecData = 0x0;
1409
1410   TH1D* hYieldPt[AliPID::kSPECIES] = {0x0, };
1411   TH1D* hYieldPtSysError[AliPID::kSPECIES] = {0x0, };
1412   TH1D* hMCgenPrimYieldPt[AliPID::kSPECIES] = {0x0, };
1413
1414   TH1D* hYieldZ[AliPID::kSPECIES] = {0x0, };
1415   TH1D* hYieldZSysError[AliPID::kSPECIES] = {0x0, };
1416   TH1D* hMCgenPrimYieldZ[AliPID::kSPECIES] = {0x0, };
1417
1418   TH1D* hYieldXi[AliPID::kSPECIES] = {0x0, };
1419   TH1D* hYieldXiSysError[AliPID::kSPECIES] = {0x0, };
1420   TH1D* hMCgenPrimYieldXi[AliPID::kSPECIES] = {0x0, };
1421
1422
1423   TH2D* hTemp = 0x0;
1424
1425   Int_t numMCgenPrimYieldHistsFound = 0;
1426
1427   Int_t lowerJetPtBinYields = -1;
1428   Int_t upperJetPtBinYields = -1;
1429
1430   Double_t actualLowerJetPtYields = -10;
1431   Double_t actualUpperJetPtYields = -10;
1432
1433   TFile* fileData = TFile::Open(pathNameData.Data());
1434   if (!fileData) {
1435     printf("Failed to open data file \"%s\"\n", pathNameData.Data());
1436     return -1;
1437   }
1438
1439
1440   hNjetsGenData = (TH2D*)fileData->Get("fh2FFJetPtGen");
1441
1442   hNjetsRecData = (TH2D*)fileData->Get("fh2FFJetPtRec");
1443
1444   if (!hNjetsRecData) {
1445     printf("Failed to load numJet histo for data!\n");
1446     return -1;
1447   }
1448
1449   const Bool_t restrictCentralityData = ((lowerCentralityData >= -1) && (upperCentralityData >= -1));
1450   // Integral(lowerCentBinLimitData, uppCentBinLimitData) will not be restricted if these values are kept
1451   const Int_t lowerCentralityBinLimitData = restrictCentralityData ? hNjetsRecData->GetXaxis()->FindBin(lowerCentralityData + 0.001)
1452                                                                    : -1;
1453   const Int_t upperCentralityBinLimitData = restrictCentralityData ? hNjetsRecData->GetXaxis()->FindBin(upperCentralityData - 0.001)
1454                                                                    : -2;
1455
1456   // Integral(lowerJetBinLimitData, uppJetBinLimitData) will not be restricted if these values are kept
1457   const Int_t lowerJetPtBinLimitData = ((lowerJetPt >= 0) && (upperJetPt >= 0))
1458                                        ? hNjetsRecData->GetYaxis()->FindBin(lowerJetPt + 0.001) : -1;
1459   const Int_t upperJetPtBinLimitData  = ((lowerJetPt >= 0) && (upperJetPt >= 0))
1460                                        ? hNjetsRecData->GetYaxis()->FindBin(upperJetPt - 0.001) : -2;
1461
1462   for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
1463     TString histName = Form("hIDFFz_%s", AliPID::ParticleShortName(species));
1464     hTemp = (TH2D*)fileData->Get(histName.Data());
1465
1466     if (hTemp && (lowerJetPtBinYields < 0 || upperJetPtBinYields < 0)) {
1467       // If the jetPt range is to be restricted
1468       if (lowerJetPt >= 0 && upperJetPt >= 0) {
1469         lowerJetPtBinYields = hTemp->GetYaxis()->FindBin(lowerJetPt + 0.001);
1470         upperJetPtBinYields = hTemp->GetYaxis()->FindBin(upperJetPt - 0.001);
1471
1472         actualLowerJetPtYields = hTemp->GetYaxis()->GetBinLowEdge(lowerJetPtBinYields);
1473         actualUpperJetPtYields = hTemp->GetYaxis()->GetBinUpEdge(upperJetPtBinYields);
1474       }
1475       else {
1476         lowerJetPtBinYields = 0;
1477         upperJetPtBinYields = -1;
1478       }
1479     }
1480
1481     if ((lowerJetPtBinYields < 0 || upperJetPtBinYields < 0) && (lowerJetPt >= 0 && upperJetPt >= 0)) {
1482       printf("Failed to determine jet pt bin limits!\n");
1483       return -1;
1484     }
1485
1486     // NOTE: The yields have already been normalised to nJets in extractFFs. However, the binning in jetPt might be different now.
1487     // Thus, in order to get the correct result, the old normalisation (just the same binnnig as in the nJet histo) has to be
1488     // undone and then it has to be normalised to the number of jets in the new binning again (see below).
1489
1490     if (hTemp) {
1491       undoJetNormalisationHist2D(hTemp, hNjetsRecData, lowerCentralityBinLimitData, upperCentralityBinLimitData);
1492       hYieldZ[species] = hTemp->ProjectionX(histName.ReplaceAll("hIDFF", "hYieldIDFF").Data(), lowerJetPtBinYields, upperJetPtBinYields,
1493                                             "e");
1494     }
1495
1496     if (!hYieldZ[species]) {
1497       printf("Failed to load hist \"%s\"\n", histName.Data());
1498       return -1;
1499     }
1500
1501
1502
1503     histName = Form("hIDFFxi_%s", AliPID::ParticleShortName(species));
1504
1505     hTemp = (TH2D*)fileData->Get(histName.Data());
1506
1507     if (hTemp) {
1508       undoJetNormalisationHist2D(hTemp, hNjetsRecData, lowerCentralityBinLimitData, upperCentralityBinLimitData);
1509       hYieldXi[species] = hTemp->ProjectionX(histName.ReplaceAll("hIDFF", "hYieldIDFF").Data(),
1510                                              lowerJetPtBinYields, upperJetPtBinYields, "e");
1511     }
1512
1513     if (!hYieldXi[species]) {
1514       printf("Failed to load hist \"%s\"\n", histName.Data());
1515       return -1;
1516     }
1517
1518     histName = Form("hIDFFtrackPt_%s", AliPID::ParticleShortName(species));
1519
1520     hTemp = (TH2D*)fileData->Get(histName.Data());
1521
1522     if (hTemp) {
1523       undoJetNormalisationHist2D(hTemp, hNjetsRecData, lowerCentralityBinLimitData, upperCentralityBinLimitData);
1524       hYieldPt[species] = hTemp->ProjectionX(histName.ReplaceAll("hIDFF", "hYieldIDFF").Data(),
1525                                              lowerJetPtBinYields, upperJetPtBinYields, "e");
1526     }
1527
1528     if (!hYieldPt[species]) {
1529       printf("Failed to load hist \"%s\"\n", histName.Data());
1530       return -1;
1531     }
1532
1533
1534
1535     histName = Form("hIDFFz_%s_sysError", AliPID::ParticleShortName(species));
1536
1537     hTemp = (TH2D*)fileData->Get(histName.Data());
1538
1539     if (hTemp) {
1540       undoJetNormalisationHist2D(hTemp, hNjetsRecData, lowerCentralityBinLimitData, upperCentralityBinLimitData);
1541       hYieldZSysError[species] = hTemp->ProjectionX(histName.ReplaceAll("hIDFF", "hYieldIDFF").Data(), lowerJetPtBinYields,
1542                                                     upperJetPtBinYields, "e");
1543     }
1544
1545     if (!hYieldZSysError[species]) {
1546       printf("Failed to load hist \"%s\"\n", histName.Data());
1547       return -1;
1548     }
1549
1550     histName = Form("hIDFFxi_%s_sysError", AliPID::ParticleShortName(species));
1551
1552     hTemp = (TH2D*)fileData->Get(histName.Data());
1553
1554     if (hTemp) {
1555       undoJetNormalisationHist2D(hTemp, hNjetsRecData, lowerCentralityBinLimitData, upperCentralityBinLimitData);
1556       hYieldXiSysError[species] = hTemp->ProjectionX(histName.ReplaceAll("hIDFF", "hYieldIDFF").Data(), lowerJetPtBinYields,
1557                                                      upperJetPtBinYields, "e");
1558     }
1559
1560     if (!hYieldXiSysError[species]) {
1561       printf("Failed to load hist \"%s\"\n", histName.Data());
1562       return -1;
1563     }
1564
1565     histName = Form("hIDFFtrackPt_%s_sysError", AliPID::ParticleShortName(species));
1566
1567     hTemp = (TH2D*)fileData->Get(histName.Data());
1568
1569     if (hTemp) {
1570       undoJetNormalisationHist2D(hTemp, hNjetsRecData, lowerCentralityBinLimitData, upperCentralityBinLimitData);
1571       hYieldPtSysError[species] = hTemp->ProjectionX(histName.ReplaceAll("hIDFF", "hYieldIDFF").Data(), lowerJetPtBinYields,
1572                                                      upperJetPtBinYields, "e");
1573     }
1574
1575     if (!hYieldPtSysError[species]) {
1576       printf("Failed to load hist \"%s\"\n", histName.Data());
1577       return -1;
1578     }
1579
1580     // In case of MC also retrieve the MC truth generated yields
1581     histName = Form("hMCgenYieldsPrimZ_%s", AliPID::ParticleShortName(species));
1582
1583     hTemp = (TH2D*)fileData->Get(histName.Data());
1584
1585     if (hTemp) {
1586       undoJetNormalisationHist2D(hTemp, hNjetsGenData, lowerCentralityBinLimitData, upperCentralityBinLimitData);
1587       hMCgenPrimYieldZ[species] = hTemp->ProjectionX(histName.ReplaceAll("hMCgen", "hYieldMCgen").Data(), lowerJetPtBinYields,
1588                                                      upperJetPtBinYields, "e");
1589     }
1590
1591
1592     histName = Form("hMCgenYieldsPrimXi_%s", AliPID::ParticleShortName(species));
1593
1594     hTemp = (TH2D*)fileData->Get(histName.Data());
1595
1596     if (hTemp) {
1597       undoJetNormalisationHist2D(hTemp, hNjetsGenData, lowerCentralityBinLimitData, upperCentralityBinLimitData);
1598       hMCgenPrimYieldXi[species] = hTemp->ProjectionX(histName.ReplaceAll("hMCgen", "hYieldMCgen").Data(), lowerJetPtBinYields,
1599                                                       upperJetPtBinYields, "e");
1600     }
1601
1602
1603     histName = Form("hMCgenYieldsPrimPt_%s", AliPID::ParticleShortName(species));
1604
1605     hTemp = (TH2D*)fileData->Get(histName.Data());
1606
1607     if (hTemp) {
1608       undoJetNormalisationHist2D(hTemp, hNjetsGenData, lowerCentralityBinLimitData, upperCentralityBinLimitData);
1609       hMCgenPrimYieldPt[species] = hTemp->ProjectionX(histName.ReplaceAll("hMCgen", "hYieldMCgen").Data(), lowerJetPtBinYields,
1610                                                       upperJetPtBinYields, "e");
1611     }
1612
1613     if (hMCgenPrimYieldPt[species] && hMCgenPrimYieldZ[species] && hMCgenPrimYieldXi[species])
1614       numMCgenPrimYieldHistsFound++;
1615   }
1616
1617   if (numMCgenPrimYieldHistsFound > 0 && numMCgenPrimYieldHistsFound != AliPID::kSPECIES) {
1618     printf("Error: Unable to retrieve all MC generated prim yield histos! Got %d.\n", numMCgenPrimYieldHistsFound);
1619     return -1;
1620   }
1621
1622
1623   /*OLD: Load from AnalysisResults.root
1624
1625   File* fileData = TFile::Open(pathNameData.Data());
1626   if (!fileData) {
1627     printf("Failed to open data file \"%s\"\n", pathNameData.Data());
1628     return -1;
1629   }
1630
1631   TString dirDataInFile = "";
1632   TDirectory* dirData = (TDirectory*)fileData->Get(fileData->GetListOfKeys()->At(0)->GetName());
1633   if (!dirData) {
1634     printf("No directory found inside data file \"%s\"\n", pathNameData.Data());
1635     return -1;
1636   }
1637
1638   TList* list = (TList*)dirData->Get(dirData->GetListOfKeys()->At(0)->GetName());
1639   if (!list) {
1640     printf("No list found in directory \"%s\" inside data file \"%s\"\n",
1641            dirData->GetListOfKeys()->At(0)->GetName(),pathNameData.Data());
1642     return -1;
1643   }
1644
1645   hNjetsGen = (TH1D*)list->FindObject("fh1FFJetPtGen");
1646   hNjetsRec = (TH1D*)list->FindObject("fh1FFJetPtRecCuts");
1647
1648   if (!hNjetsRec) {
1649     printf("Failed to load number of jets (rec) histo!\n");
1650     return -1;
1651   }
1652
1653
1654   for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
1655     TString histName = Form("fh2FFZRecCuts_%s", AliPID::ParticleName(species));
1656     hTemp = (TH2D*)list->FindObject(histName.Data());
1657
1658     if (hTemp && (lowerJetPtBinYields < 0 || upperJetPtBinYields < 0)) {
1659       // If the jetPt range is to be restricted
1660       if (lowerJetPt >= 0 && upperJetPt >= 0) {
1661         lowerJetPtBinYields = hTemp->GetXaxis()->FindBin(lowerJetPt + 0.001);
1662         upperJetPtBinYields = hTemp->GetXaxis()->FindBin(upperJetPt - 0.001);
1663
1664         actualLowerJetPtYields = hTemp->GetXaxis()->GetBinLowEdge(lowerJetPtBinYields);
1665         actualUpperJetPtYields = hTemp->GetXaxis()->GetBinUpEdge(upperJetPtBinYields);
1666       }
1667       else {
1668         lowerJetPtBinYields = 0;
1669         upperJetPtBinYields = -1;
1670       }
1671     }
1672
1673     if ((lowerJetPtBinYields < 0 || upperJetPtBinYields < 0) && (lowerJetPt >= 0 && upperJetPt >= 0)) {
1674       printf("Failed to determine jet pt bin limits!\n");
1675       return -1;
1676     }
1677
1678     if (hTemp)
1679       hYieldZ[species] = hTemp->ProjectionY(histName.ReplaceAll("fh2", "fh").Data(), lowerJetPtBinYields, upperJetPtBinYields, "e");
1680
1681     if (!hYieldZ[species]) {
1682       printf("Failed to load hist \"%s\"\n", histName.Data());
1683       return -1;
1684     }
1685
1686
1687
1688     histName = Form("fh2FFXiRecCuts_%s", AliPID::ParticleName(species));
1689
1690     hTemp = (TH2D*)list->FindObject(histName.Data());
1691
1692     if (hTemp)
1693       hYieldXi[species] = hTemp->ProjectionY(histName.ReplaceAll("fh2", "fh").Data(), lowerJetPtBinYields, upperJetPtBinYields, "e");
1694
1695     if (!hYieldXi[species]) {
1696       printf("Failed to load hist \"%s\"\n", histName.Data());
1697       return -1;
1698     }
1699
1700     histName = Form("fh2FFTrackPtRecCuts_%s", AliPID::ParticleName(species));
1701
1702     hTemp = (TH2D*)list->FindObject(histName.Data());
1703
1704     if (hTemp)
1705       hYieldPt[species] = hTemp->ProjectionY(histName.ReplaceAll("fh2", "fh").Data(), lowerJetPtBinYields, upperJetPtBinYields, "e");
1706
1707     if (!hYieldPt[species]) {
1708       printf("Failed to load hist \"%s\"\n", histName.Data());
1709       return -1;
1710     }
1711
1712     // In case of MC also retrieve the MC truth generated yields
1713     histName = Form("fh2FFZGen_%s", AliPID::ParticleName(species));
1714
1715     hTemp = (TH2D*)list->FindObject(histName.Data());
1716     if (hTemp)
1717       hMCgenPrimYieldZ[species] = hTemp->ProjectionY(histName.ReplaceAll("fh2", "fh").Data(), lowerJetPtBinYields, upperJetPtBinYields,
1718                                                      "e");
1719
1720
1721     histName = Form("fh2FFXiGen_%s", AliPID::ParticleName(species));
1722
1723     hTemp = (TH2D*)list->FindObject(histName.Data());
1724     if (hTemp)
1725       hMCgenPrimYieldXi[species] = hTemp->ProjectionY(histName.ReplaceAll("fh2", "fh").Data(), lowerJetPtBinYields, upperJetPtBinYields,
1726                                                       "e");
1727
1728
1729     histName = Form("fh2FFTrackPtGen_%s", AliPID::ParticleName(species));
1730
1731     hTemp = (TH2D*)list->FindObject(histName.Data());
1732     if (hTemp)
1733       hMCgenPrimYieldPt[species] = hTemp->ProjectionY(histName.ReplaceAll("fh2", "fh").Data(), lowerJetPtBinYields, upperJetPtBinYields,
1734                                                       "e");
1735
1736     if (hMCgenPrimYieldPt[species] && hMCgenPrimYieldZ[species] && hMCgenPrimYieldXi[species])
1737       numMCgenPrimYieldHistsFound++;
1738   }
1739
1740   if (numMCgenPrimYieldHistsFound > 0 && numMCgenPrimYieldHistsFound != AliPID::kSPECIES) {
1741     printf("Error: Unable to retrieve all MC generated prim yield histos! Got %d.\n", numMCgenPrimYieldHistsFound);
1742     return -1;
1743   }
1744
1745
1746   const Double_t nJetsGen = hNjetsGen ? hNjetsGen->Integral(lowerJetPtBinYields, upperJetPtBinYields) : -1;
1747   const Double_t nJetsRec = hNjetsRec->Integral(lowerJetPtBinYields, upperJetPtBinYields);
1748   */
1749
1750   // Take binning from pion yield (binning for all species the same) and create a new AliCFContainer with this new binning
1751   TH1D hDummyPt(*hYieldPt[AliPID::kPion]);
1752   hDummyPt.SetName("hDummyPt");
1753   TAxis* axisPt = 0x0;
1754
1755   if (rebinEfficiencyPt > 1) {
1756     Int_t nBinsNew = 0;
1757     const Double_t* newBins = GetBins(rebinEfficiencyPt, nBinsNew);
1758
1759     hDummyPt.SetBins(nBinsNew, newBins);
1760
1761     axisPt = hDummyPt.GetXaxis();
1762
1763     //axisPt = hDummyPt.Rebin(rebinEfficiencyPt, "", 0)->GetXaxis();
1764   }
1765   else
1766     axisPt = hDummyPt.GetXaxis();
1767
1768   const TArrayD* binsPtCurrent = axisPt->GetXbins();
1769   TArrayD* binsPtNew = new TArrayD(*binsPtCurrent);
1770
1771   TH1D hDummyZ(*hYieldZ[AliPID::kPion]);
1772   hDummyZ.SetName("hDummyZ");
1773   TAxis* axisZ = 0x0;
1774
1775   if (rebinEfficiencyZ > 1)
1776     axisZ = hDummyZ.Rebin(rebinEfficiencyZ, "", 0)->GetXaxis();
1777   else
1778     axisZ = hDummyZ.GetXaxis();
1779
1780   const TArrayD* binsZCurrent = axisZ->GetXbins();
1781   TArrayD* binsZNew = new TArrayD(*binsZCurrent);
1782
1783   TH1D hDummyXi(*hYieldXi[AliPID::kPion]);
1784   hDummyXi.SetName("hDummyXi");
1785   TAxis* axisXi = 0x0;
1786
1787   if (rebinEfficiencyXi > 1)
1788     axisXi = hDummyXi.Rebin(rebinEfficiencyXi, "", 0)->GetXaxis();
1789   else
1790     axisXi = hDummyXi.GetXaxis();
1791
1792   const TArrayD* binsXiCurrent = axisXi->GetXbins();
1793   TArrayD* binsXiNew = new TArrayD(*binsXiCurrent);
1794
1795
1796   const Int_t nEffDims = data->GetNVar();
1797   Int_t nEffBins[nEffDims];
1798
1799   for (Int_t iDim = 0; iDim < nEffDims; iDim++) {
1800     if (iDim == iPt)
1801       nEffBins[iDim] = axisPt->GetNbins();
1802     else if (iDim == iZ)
1803       nEffBins[iDim] = axisZ->GetNbins();
1804     else if (iDim == iXi)
1805       nEffBins[iDim] = axisXi->GetNbins();
1806     else
1807       nEffBins[iDim] = data->GetNBins(iDim);
1808   }
1809
1810
1811   // Just make one large pT bin above some threshold, if desired
1812   if (binsPtNew->fN != 0 && constantCorrectionAbovePtThreshold > 0) {
1813     for (Int_t iBin = 0; iBin <= nEffBins[iPt] + 1; iBin++) {
1814       // Find the first bin edged really larger than the threshold.
1815       // If the bin edge before equals the threshold, just set the
1816       // current bin edge to the right end of the spectrum -> Done.
1817       // If the bin edge before is different, set the bin edge to the
1818       // threshold
1819       if (binsPtNew->fArray[iBin] > constantCorrectionAbovePtThreshold) {
1820         if (binsPtNew->fArray[iBin - 1] == constantCorrectionAbovePtThreshold) {
1821           binsPtNew->fArray[iBin] = binsPtNew->fArray[nEffBins[iPt]];
1822           nEffBins[iPt] = iBin;
1823           break;
1824         }
1825         else {
1826           binsPtNew->fArray[iBin] = constantCorrectionAbovePtThreshold;
1827         }
1828       }
1829     }
1830   }
1831
1832   // Just make one large Xi bin above some threshold, if desired
1833   if (binsXiNew->fN != 0 && constantCorrectionAboveXiThreshold > 0) {
1834     for (Int_t iBin = 0; iBin <= nEffBins[iXi] + 1; iBin++) {
1835       // Find the first bin edged really larger than the threshold.
1836       // If the bin edge before equals the threshold, just set the
1837       // current bin edge to the right end of the spectrum -> Done.
1838       // If the bin edge before is different, set the bin edge to the
1839       // threshold
1840       if (binsXiNew->fArray[iBin] > constantCorrectionAboveXiThreshold) {
1841         if (binsXiNew->fArray[iBin - 1] == constantCorrectionAboveXiThreshold) {
1842           binsXiNew->fArray[iBin] = binsXiNew->fArray[nEffBins[iXi]];
1843           nEffBins[iXi] = iBin;
1844           break;
1845         }
1846         else {
1847           binsXiNew->fArray[iBin] = constantCorrectionAboveXiThreshold;
1848         }
1849       }
1850     }
1851   }
1852
1853   AliCFContainer *dataRebinned = new AliCFContainer(Form("%s_rebinned", data->GetName()), Form("%s (rebinned)", data->GetTitle()),
1854                                                     data->GetNStep(), nEffDims, nEffBins);
1855
1856   for (Int_t iDim = 0; iDim < nEffDims; iDim++) {
1857     dataRebinned->SetVarTitle(iDim, data->GetVarTitle(iDim));
1858
1859     if (iDim == iPt) {
1860       if (binsPtNew->fN == 0)
1861         dataRebinned->SetBinLimits(iDim, axisPt->GetXmin(), axisPt->GetXmax());
1862       else
1863         dataRebinned->SetBinLimits(iDim, binsPtNew->fArray);
1864     }
1865     else if (iDim == iZ) {
1866       if (binsZNew->fN == 0)
1867         dataRebinned->SetBinLimits(iDim, axisZ->GetXmin(), axisZ->GetXmax());
1868       else
1869         dataRebinned->SetBinLimits(iDim, binsZNew->fArray);
1870     }
1871     else if (iDim == iXi) {
1872       if (binsXiNew->fN == 0)
1873         dataRebinned->SetBinLimits(iDim, axisXi->GetXmin(), axisXi->GetXmax());
1874       else
1875         dataRebinned->SetBinLimits(iDim, binsXiNew->fArray);
1876     }
1877     else {
1878       dataRebinned->SetBinLimits(iDim, data->GetBinLimits(iDim));
1879     }
1880   }
1881
1882   for (Int_t iStep = 0; iStep < data->GetNStep(); iStep++)
1883     dataRebinned->SetStepTitle(iStep, data->GetStepTitle(iStep));
1884
1885   Int_t coord[nEffDims];
1886   Double_t binCenterCoord[nEffDims];
1887
1888   // Fill content from old grid into the new grid with proper binning
1889   for (Int_t iStep = 0; iStep < data->GetNStep(); iStep++) {
1890     Long64_t nBinsGrid = data->GetGrid(iStep)->GetGrid()->GetNbins();
1891
1892     for (Long64_t iBin = 0; iBin <= nBinsGrid + 1; iBin++) {
1893       Double_t binContent = data->GetGrid(iStep)->GetGrid()->GetBinContent(iBin, coord);
1894       Double_t binError2  = data->GetGrid(iStep)->GetGrid()->GetBinError2(iBin);
1895
1896       for (Int_t iDim = 0; iDim < nEffDims; iDim++) {
1897         binCenterCoord[iDim] = data->GetBinCenter(iDim, coord[iDim]);
1898       }
1899
1900       Long64_t iBinRebinned = dataRebinned->GetGrid(iStep)->GetGrid()->GetBin(binCenterCoord);
1901       dataRebinned->GetGrid(iStep)->GetGrid()->AddBinContent(iBinRebinned, binContent);
1902       dataRebinned->GetGrid(iStep)->GetGrid()->AddBinError2(iBinRebinned, binError2);
1903     }
1904   }
1905
1906
1907   // If desired, restrict centrality axis
1908   Int_t lowerCentralityBinLimit = -1;
1909   Int_t upperCentralityBinLimit = -2; // Integral(lowerCentBinLimit, uppCentBinLimit) will not be restricted if these values are kept
1910   Bool_t restrictCentralityAxis = kFALSE;
1911   Double_t actualLowerCentrality = -1.;
1912   Double_t actualUpperCentrality = -1.;
1913
1914   if (lowerCentrality >= -1 && upperCentrality >= -1) {
1915     // Add subtract a very small number to avoid problems with values right on the border between to bins
1916     lowerCentralityBinLimit = dataRebinned->GetAxis(iMult, 0)->FindBin(lowerCentrality + 0.001);
1917     upperCentralityBinLimit = dataRebinned->GetAxis(iMult, 0)->FindBin(upperCentrality - 0.001);
1918
1919     // Check if the values look reasonable
1920     if (lowerCentralityBinLimit <= upperCentralityBinLimit && lowerCentralityBinLimit >= 1
1921         && upperCentralityBinLimit <= dataRebinned->GetAxis(iMult, 0)->GetNbins()) {
1922       actualLowerCentrality = dataRebinned->GetAxis(iMult, 0)->GetBinLowEdge(lowerCentralityBinLimit);
1923       actualUpperCentrality = dataRebinned->GetAxis(iMult, 0)->GetBinUpEdge(upperCentralityBinLimit);
1924
1925       restrictCentralityAxis = kTRUE;
1926     }
1927     else {
1928       std::cout << std::endl;
1929       std::cout << "Requested centrality range out of limits or upper and lower limit are switched!" << std::endl;
1930       return -1;
1931     }
1932   }
1933
1934   std::cout << "centrality: ";
1935   if (restrictCentralityAxis) {
1936     std::cout << actualLowerCentrality << " - " << actualUpperCentrality << std::endl;
1937     dataRebinned->SetRangeUser(iMult, lowerCentralityBinLimit, upperCentralityBinLimit, kTRUE);
1938     data->SetRangeUser(iMult, lowerCentralityBinLimit, upperCentralityBinLimit, kTRUE);
1939   }
1940   else {
1941     std::cout << "All" << std::endl;
1942   }
1943
1944
1945
1946   // If desired, restrict jetPt axis
1947   Int_t lowerJetPtBinLimit = -1;
1948   Int_t upperJetPtBinLimit = -2; // Integral(lowerJetBinLimit, uppJetBinLimit) will not be restricted if these values are kept
1949   Bool_t restrictJetPtAxis = kFALSE;
1950   Double_t actualLowerJetPt = -1.;
1951   Double_t actualUpperJetPt = -1.;
1952
1953   if (lowerJetPt >= 0 && upperJetPt >= 0) {
1954     // Add subtract a very small number to avoid problems with values right on the border between to bins
1955     lowerJetPtBinLimit = dataRebinned->GetAxis(iJetPt, 0)->FindBin(lowerJetPt + 0.001);
1956     upperJetPtBinLimit = dataRebinned->GetAxis(iJetPt, 0)->FindBin(upperJetPt - 0.001);
1957
1958     // Check if the values look reasonable
1959     if (lowerJetPtBinLimit <= upperJetPtBinLimit && lowerJetPtBinLimit >= 1 &&
1960         upperJetPtBinLimit <= dataRebinned->GetAxis(iJetPt, 0)->GetNbins()) {
1961       actualLowerJetPt = dataRebinned->GetAxis(iJetPt, 0)->GetBinLowEdge(lowerJetPtBinLimit);
1962       actualUpperJetPt = dataRebinned->GetAxis(iJetPt, 0)->GetBinUpEdge(upperJetPtBinLimit);
1963
1964       restrictJetPtAxis = kTRUE;
1965     }
1966     else {
1967       std::cout << std::endl;
1968       std::cout << "Requested jet pT range out of limits or upper and lower limit are switched!" << std::endl;
1969       return -1;
1970     }
1971   }
1972
1973   std::cout << "jet pT: ";
1974   if (restrictJetPtAxis) {
1975     std::cout << actualLowerJetPt << " - " << actualUpperJetPt << std::endl;
1976     dataRebinned->SetRangeUser(iJetPt, lowerJetPtBinLimit, upperJetPtBinLimit, kTRUE);
1977     data->SetRangeUser(iJetPt, lowerJetPtBinLimit, upperJetPtBinLimit, kTRUE);
1978   }
1979   else {
1980     std::cout << "All" << std::endl;
1981   }
1982
1983   if (restrictJetPtAxis &&
1984       (TMath::Abs(actualLowerJetPt - actualLowerJetPtYields) > 1e-3 || TMath::Abs(actualUpperJetPt - actualUpperJetPtYields) > 1e-3)) {
1985     std::cout << "ERROR: Disagreement between jetPt bounds of data and efficiency!" << std::endl;
1986     return -1;
1987   }
1988
1989   // If desired, restrict charge axis
1990   std::cout << "Charge selection (efficiency): ";
1991   if (chargeMode == kAllCharged)
1992     std::cout << "All charged particles" << std::endl;
1993   else if (chargeMode == kNegCharge)
1994     std::cout << "Negative particles only" << std::endl;
1995   else if (chargeMode == kPosCharge)
1996     std::cout << "Positive particles only" << std::endl;
1997   else {
1998     std::cout << "Unknown -> ERROR" << std::endl;
1999     return -1;
2000   }
2001
2002   const Bool_t restrictCharge = (chargeMode != kAllCharged);
2003
2004   Int_t lowerChargeBinLimit = -1;
2005   Int_t upperChargeBinLimit = -2;
2006   Double_t actualLowerCharge = -999;
2007   Double_t actualUpperCharge = -999;
2008
2009   if (restrictCharge) {
2010     // Add subtract a very small number to avoid problems with values right on the border between to bins
2011     if (chargeMode == kNegCharge) {
2012       lowerChargeBinLimit = dataRebinned->GetAxis(iCharge, 0)->FindBin(-1. + 0.001);
2013       upperChargeBinLimit = dataRebinned->GetAxis(iCharge, 0)->FindBin(0. - 0.001);
2014     }
2015     else if (chargeMode == kPosCharge) {
2016       lowerChargeBinLimit = dataRebinned->GetAxis(iCharge, 0)->FindBin(0. + 0.001);
2017       upperChargeBinLimit = dataRebinned->GetAxis(iCharge, 0)->FindBin(1. - 0.001);
2018     }
2019
2020     // Check if the values look reasonable
2021     if (lowerChargeBinLimit <= upperChargeBinLimit && lowerChargeBinLimit >= 1
2022         && upperChargeBinLimit <= dataRebinned->GetAxis(iCharge, 0)->GetNbins()) {
2023       actualLowerCharge = dataRebinned->GetAxis(iCharge, 0)->GetBinLowEdge(lowerChargeBinLimit);
2024       actualUpperCharge = dataRebinned->GetAxis(iCharge, 0)->GetBinUpEdge(upperChargeBinLimit);
2025
2026       std::cout << "Charge range (efficiency): " << actualLowerCharge << " - " << actualUpperCharge << std::endl;
2027     }
2028     else {
2029       std::cout << std::endl;
2030       std::cout << "Requested charge range (efficiency) out of limits or upper and lower limit are switched!" << std::endl;
2031       return -1;
2032     }
2033
2034     dataRebinned->SetRangeUser(iCharge, lowerChargeBinLimit, upperChargeBinLimit, kTRUE);
2035     data->SetRangeUser(iCharge, lowerChargeBinLimit, upperChargeBinLimit, kTRUE);
2036   }
2037
2038   std::cout << std::endl;
2039
2040
2041   // Load numJet histos from bhess_PID*.root file related to efficiency file.
2042   TString pathNameDataMC = pathNameEfficiency;
2043   pathNameDataMC.ReplaceAll("_efficiency", "");
2044
2045   TFile* fDataMC = TFile::Open(pathNameDataMC.Data());
2046   if (!fDataMC)  {
2047     std::cout << std::endl;
2048     std::cout << "Failed to open file \"" << pathNameData.Data() << "\" to obtain num of rec/gen jets!" << std::endl;
2049
2050     return -1;
2051   }
2052
2053   TString listName = pathNameDataMC;
2054   listName.Replace(0, listName.Last('/') + 1, "");
2055   listName.ReplaceAll(".root", "");
2056
2057   TObjArray* histList = (TObjArray*)(fDataMC->Get(listName.Data()));
2058   if (!histList) {
2059     std::cout << std::endl;
2060     std::cout << "Failed to load list \"" << listName.Data() << "\" to obtain num of rec/gen jets!" << std::endl;
2061     return -1;
2062   }
2063
2064   hNjetsGen = (TH2D*)histList->FindObject("fh2FFJetPtGen");
2065   hNjetsRec = (TH2D*)histList->FindObject("fh2FFJetPtRec");
2066
2067   if (!hNjetsRec || !hNjetsGen) {
2068     std::cout << "Failed to load number of jets histos!" << std::endl;
2069
2070
2071     // For backward compatibility (TODO REMOVE IN FUTURE): Load info from fixed AnalysisResults file (might be wrong, if other
2072     // period is considered; also: No multiplicity information)
2073     TString pathEfficiency = pathNameEfficiency;
2074     pathEfficiency.Replace(pathEfficiency.Last('/'), pathEfficiency.Length(), "");
2075     TString pathBackward = Form("%s/AnalysisResults.root", pathEfficiency.Data());
2076     TFile* fBackward = TFile::Open(pathBackward.Data());
2077
2078     TString dirDataInFile = "";
2079     TDirectory* dirData = fBackward ? (TDirectory*)fBackward->Get(fBackward->GetListOfKeys()->At(0)->GetName()) : 0x0;
2080
2081     TList* list = dirData ? (TList*)dirData->Get(dirData->GetListOfKeys()->At(0)->GetName()) : 0x0;
2082
2083     TH1D* hFFJetPtRec = list ? (TH1D*)list->FindObject("fh1FFJetPtRecCutsInc") : 0x0;
2084     TH1D* hFFJetPtGen = list ? (TH1D*)list->FindObject("fh1FFJetPtGenInc") : 0x0;
2085
2086     if (hFFJetPtRec && hFFJetPtGen) {
2087       printf("***WARNING: For backward compatibility, using file \"%s\" to get number of jets. BUT: Might be wrong period and has no mult info!***\n",
2088         pathBackward.Data());
2089       printf("ALSO: Using Njets for inclusive jets!!!!\n");
2090
2091       hNjetsRec = new TH2D("fh2FFJetPtRec", "", 1, -1, 1, dataRebinned->GetAxis(iJetPt, 0)->GetNbins(),
2092                           dataRebinned->GetAxis(iJetPt, 0)->GetXbins()->GetArray());
2093
2094       for (Int_t iJet = 1; iJet <= hNjetsRec->GetNbinsY(); iJet++) {
2095         Int_t lowerBin = hFFJetPtRec->FindBin(hNjetsRec->GetYaxis()->GetBinLowEdge(iJet) + 1e-3);
2096         Int_t upperBin = hFFJetPtRec->FindBin(hNjetsRec->GetYaxis()->GetBinUpEdge(iJet) - 1e-3);
2097         hNjetsRec->SetBinContent(1, iJet, hFFJetPtRec->Integral(lowerBin, upperBin));
2098       }
2099
2100       hNjetsGen = new TH2D("fh2FFJetPtGen", "", 1, -1, 1,  dataRebinned->GetAxis(iJetPt, 0)->GetNbins(),
2101                           dataRebinned->GetAxis(iJetPt, 0)->GetXbins()->GetArray());
2102
2103       for (Int_t iJet = 1; iJet <= hNjetsGen->GetNbinsY(); iJet++) {
2104         Int_t lowerBin = hFFJetPtGen->FindBin(hNjetsGen->GetYaxis()->GetBinLowEdge(iJet) + 1e-3);
2105         Int_t upperBin = hFFJetPtGen->FindBin(hNjetsGen->GetYaxis()->GetBinUpEdge(iJet) - 1e-3);
2106         hNjetsGen->SetBinContent(1, iJet, hFFJetPtGen->Integral(lowerBin, upperBin));
2107       }
2108     }
2109
2110     if (!hNjetsRec || ! hNjetsGen)
2111       return -1;
2112   }
2113
2114   // For normalisation to number of jets
2115   // NOTE: These numbers are for the efficiency only! The data will be normalised to its own number!!!
2116   const Double_t nJetsGen = hNjetsGen ? hNjetsGen->Integral(lowerCentralityBinLimit, upperCentralityBinLimit, lowerJetPtBinLimit,
2117                                                             upperJetPtBinLimit) : 1.;
2118   const Double_t nJetsRec = hNjetsRec ? hNjetsRec->Integral(lowerCentralityBinLimit, upperCentralityBinLimit, lowerJetPtBinLimit,
2119                                                             upperJetPtBinLimit) : 1.;
2120
2121
2122
2123   // Secondary correction
2124   AliCFEffGrid* sec = new AliCFEffGrid("sec", "Secondary Contamination", *dataRebinned);
2125   sec->CalculateEfficiency(kStepRecWithRecCutsMeasuredObsPrimaries, kStepRecWithRecCutsMeasuredObs);
2126
2127   TH1 *hSecAllPt = 0x0, *hSecID2Pt = 0x0;
2128   TH1D* hSecPt[AliPID::kSPECIES] = { 0x0, };
2129   TH1D* hSecToPiRatioPt[AliPID::kSPECIES] = { 0x0, };
2130   TCanvas *cSecPt = 0x0;
2131   TCanvas *cSecToPiRatioPt = 0x0;
2132
2133   TH1 *hSecAllZ = 0x0, *hSecID2Z = 0x0;
2134   TH1D* hSecZ[AliPID::kSPECIES] = { 0x0, };
2135   TH1D* hSecToPiRatioZ[AliPID::kSPECIES] = { 0x0, };
2136   TCanvas *cSecZ = 0x0;
2137   TCanvas *cSecToPiRatioZ = 0x0;
2138
2139   TH1 *hSecAllXi = 0x0, *hSecID2Xi = 0x0;
2140   TH1D* hSecXi[AliPID::kSPECIES] = { 0x0, };
2141   TH1D* hSecToPiRatioXi[AliPID::kSPECIES] = { 0x0, };
2142   TCanvas *cSecXi = 0x0;
2143   TCanvas *cSecToPiRatioXi = 0x0;
2144
2145   // Secondary correction with strangeness scaling
2146   AliCFEffGrid* secStrangeScale = new AliCFEffGrid("secStrangeScale", "Secondary Contamination with strangeness scaling",
2147                                                    *dataRebinned);
2148   secStrangeScale->CalculateEfficiency(kStepRecWithRecCutsMeasuredObsPrimaries, kStepRecWithRecCutsMeasuredObsStrangenessScaled);
2149
2150   TH1 *hSecSSAllPt = 0x0, *hSecSSID2Pt = 0x0;
2151   TH1D* hSecSSPt[AliPID::kSPECIES] = { 0x0, };
2152   TH1D* hSecSSToPiRatioPt[AliPID::kSPECIES] = { 0x0, };
2153   TCanvas *cSecSSPt = 0x0;
2154   TCanvas *cSecSSToPiRatioPt = 0x0;
2155
2156   TH1 *hSecSSAllZ = 0x0, *hSecSSID2Z = 0x0;
2157   TH1D* hSecSSZ[AliPID::kSPECIES] = { 0x0, };
2158   TH1D* hSecSSToPiRatioZ[AliPID::kSPECIES] = { 0x0, };
2159   TCanvas *cSecSSZ = 0x0;
2160   TCanvas *cSecSSToPiRatioZ = 0x0;
2161
2162   TH1 *hSecSSAllXi = 0x0, *hSecSSID2Xi = 0x0;
2163   TH1D* hSecSSXi[AliPID::kSPECIES] = { 0x0, };
2164   TH1D* hSecSSToPiRatioXi[AliPID::kSPECIES] = { 0x0, };
2165   TCanvas *cSecSSXi = 0x0;
2166   TCanvas *cSecSSToPiRatioXi = 0x0;
2167
2168
2169   const Double_t upperTrackPt = restrictJetPtAxis ? actualUpperJetPt : 50.;
2170
2171   TString strangenessString[2] = { "", "SS" };
2172
2173   // Get the secondary contamination vs. pT,z,xi for each species w/ and w/o strangeness scaling
2174
2175   for (Int_t type = 0; type < kNtypes; type++) {
2176     for (Int_t strangenessScaling = 0; strangenessScaling < 2; strangenessScaling++) {
2177       TH1 **hSecAll, **hSecID2;
2178       TH1D **hSec, **hSecToPiRatio;
2179       TCanvas **cSec, **cSecToPiRatio;
2180
2181       Int_t iDimProj = -1;
2182
2183       AliCFEffGrid* secCurr = strangenessScaling ? secStrangeScale : sec;
2184
2185       if (type == kTrackPt) {
2186         hSecAll = strangenessScaling ? &hSecSSAllPt : &hSecAllPt;
2187         hSecID2 = strangenessScaling ? &hSecSSID2Pt : &hSecID2Pt;
2188         hSec = strangenessScaling ? &hSecSSPt[0] : &hSecPt[0];
2189         hSecToPiRatio = strangenessScaling ? &hSecSSToPiRatioPt[0] : &hSecToPiRatioPt[0];
2190         cSec = strangenessScaling ? &cSecSSPt : &cSecPt;
2191         cSecToPiRatio = strangenessScaling ? &cSecSSToPiRatioPt : &cSecToPiRatioPt;
2192         iDimProj = iPt;
2193       }
2194       else if (type == kZ) {
2195         hSecAll = strangenessScaling ? &hSecSSAllZ : &hSecAllZ;
2196         hSecID2 = strangenessScaling ? &hSecSSID2Z : &hSecID2Z;
2197         hSec = strangenessScaling ? &hSecSSZ[0] : &hSecZ[0];
2198         hSecToPiRatio = strangenessScaling ? &hSecSSToPiRatioZ[0] : &hSecToPiRatioZ[0];
2199         cSec = strangenessScaling ? &cSecSSZ : &cSecZ;
2200         cSecToPiRatio = strangenessScaling ? &cSecSSToPiRatioZ : &cSecToPiRatioZ;
2201         iDimProj = iZ;
2202       }
2203       else if (type == kXi) {
2204         hSecAll = strangenessScaling ? &hSecSSAllXi : &hSecAllXi;
2205         hSecID2 = strangenessScaling ? &hSecSSID2Xi : &hSecID2Xi;
2206         hSec = strangenessScaling ? &hSecSSXi[0] : &hSecXi[0];
2207         hSecToPiRatio = strangenessScaling ? &hSecSSToPiRatioXi[0] : &hSecToPiRatioXi[0];
2208         cSec = strangenessScaling ? &cSecSSXi : &cSecXi;
2209         cSecToPiRatio = strangenessScaling ? &cSecSSToPiRatioXi : &cSecToPiRatioXi;
2210         iDimProj = iXi;
2211       }
2212       else
2213         continue;
2214
2215       *hSecAll = secCurr->Project(iDimProj);
2216       (*hSecAll)->SetName(Form("hSec%s_%s_all", strangenessString[strangenessScaling].Data(), titles[type].Data()));
2217       (*hSecAll)->SetTitle("All");
2218       (*hSecAll)->SetMarkerStyle(24);
2219       (*hSecAll)->SetLineWidth(2.);
2220       (*hSecAll)->SetStats(kFALSE);
2221       (*hSecAll)->GetYaxis()->SetTitle("Primary Fraction");
2222       (*hSecID2) = (TH2D*)secCurr->Project(iDimProj, iMCid);
2223       (*hSecID2)->SetName(Form("hSecID2_%s", titles[type].Data()));
2224
2225       for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
2226         hSec[species] = ((TH2D*)*hSecID2)->ProjectionX(Form("hSec%s_%s_%s", strangenessString[strangenessScaling].Data(),
2227                                                             titles[type].Data(), AliPID::ParticleShortName(species)),
2228                                                       species + 1, species + 1, "e");
2229         hSec[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species)));
2230         hSec[species]->SetLineColor(getLineColorAliPID(species));
2231         hSec[species]->SetMarkerColor(getLineColorAliPID(species));
2232         hSec[species]->SetMarkerStyle(24);
2233         hSec[species]->SetLineWidth(2.);
2234         hSec[species]->GetYaxis()->SetRangeUser(0., 1.01);
2235         hSec[species]->SetStats(kFALSE);
2236         hSec[species]->GetYaxis()->SetTitle("Primary Fraction");
2237       }
2238
2239       *cSec = new TCanvas(Form("cSec%s_%s", strangenessString[strangenessScaling].Data(), titles[type].Data()),
2240                           Form("Primary fraction for different species%s", strangenessScaling ? " (strangeness scaled)" : ""),
2241                           100, 10, 1200, 800);
2242       (*cSec)->cd();
2243       (*cSec)->SetGridx(1);
2244       (*cSec)->SetGridy(1);
2245
2246       if (type == kTrackPt)
2247         (*hSecAll)->GetXaxis()->SetRangeUser(0.15, upperTrackPt);
2248
2249       (*hSecAll)->Draw("E1");
2250
2251       for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2252         if (type == kTrackPt)
2253           hSec[i]->GetXaxis()->SetRangeUser(0.15, upperTrackPt);
2254
2255         hSec[i]->Draw("E1 same");
2256       }
2257       (*cSec)->BuildLegend()->SetFillColor(kWhite);
2258
2259       ClearTitleFromHistoInCanvas(*cSec);
2260
2261       // To-pi ratios
2262       for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
2263         if (species == AliPID::kPion)
2264           continue; // Do not consider pion-to-pion ratio
2265
2266         hSecToPiRatio[species] = new TH1D(*hSec[species]);
2267         hSecToPiRatio[species]->Reset();
2268         hSecToPiRatio[species]->SetName(Form("hSec%sToPionRatio_%s_%s", strangenessString[strangenessScaling].Data(),
2269                                              titles[type].Data(), AliPID::ParticleShortName(species)));
2270         hSecToPiRatio[species]->SetTitle(Form("%s/#pi", AliPID::ParticleLatexName(species)));
2271         hSecToPiRatio[species]->SetLineColor(getLineColorAliPID(species));
2272         hSecToPiRatio[species]->SetMarkerColor(getLineColorAliPID(species));
2273         hSecToPiRatio[species]->SetMarkerStyle(24);
2274         hSecToPiRatio[species]->SetLineWidth(2.);
2275         hSecToPiRatio[species]->SetStats(kFALSE);
2276         hSecToPiRatio[species]->GetYaxis()->SetTitle("Primary Fraction of Ratio");
2277
2278         // Samples for different species are independent, so just divide correction factors
2279         hSecToPiRatio[species]->Divide(hSec[species], hSec[AliPID::kPion], 1., 1., "");
2280         hSecToPiRatio[species]->GetYaxis()->SetRangeUser(0., 1.1);
2281       }
2282
2283       *cSecToPiRatio = new TCanvas(Form("cSec%sToPiRatio_%s", strangenessString[strangenessScaling].Data(), titles[type].Data()),
2284                                    Form("Primary fraction of to-#pi-ratio for different species%s",
2285                                         strangenessScaling ? " (strangeness scaled)" : ""),
2286                                    100, 10, 1200, 800);
2287       (*cSecToPiRatio)->cd();
2288       (*cSecToPiRatio)->SetGridx(1);
2289       (*cSecToPiRatio)->SetGridy(1);
2290
2291       for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2292         if (i == AliPID::kPion)
2293           continue;
2294
2295         if (type == kTrackPt)
2296           hSecToPiRatio[i]->GetXaxis()->SetRangeUser(0.15, upperTrackPt);
2297
2298         hSecToPiRatio[i]->Draw(Form("E1%s", i == 0 ? "" : " same"));
2299       }
2300       (*cSecToPiRatio)->BuildLegend()->SetFillColor(kWhite);
2301
2302       ClearTitleFromHistoInCanvas(*cSecToPiRatio);
2303     }
2304   }
2305
2306
2307   // Construct the efficiency grid from the data container
2308   // NOTE: Here, the scale factors 1/nJet are different. This effects the error of the ratio.
2309   // Thus, the histograms must be extracted, scaled and THEN divided to arrive at the correction
2310   // factor and the standard "calculateEfficiency" function CANNOT be used (multiplying the ratio
2311   // of the scale factor afterwards would result in wrong errors)
2312
2313
2314
2315   // If desired, apply the GEANT/FLUKA correction first
2316   // NOTE: This will change dataRebinned! If anything else should also be calculated, it must be done before.
2317   // Otherwise, the GEANT/FLUKA correction factor for the efficiency will also affect it!
2318   const Int_t genStepEff = kStepGenWithGenCuts;
2319   if (correctGeantFluka) {
2320     printf("Applying GEANT/FLUKA correction...\n");
2321     if (!geantFlukaCorrection(dataRebinned, genStepEff)) {
2322       printf("GEANT/FLUKA correction could not be applied!\n");
2323       return kFALSE;
2324     }
2325   }
2326
2327
2328   // Either one can take kStepRecWithGenCutsMeasuredObs or, what I prefer, one can take
2329   // kStepRecWithRecCutsMeasuredObsPrimaries => The difference is only the eta cut, which is on the rec level
2330   // in the latter case, i.e. one corrects for eta resolution (although the effect should be very small)
2331   const Int_t recStepEff = kStepRecWithRecCutsMeasuredObsPrimaries;
2332
2333   // The efficiency vs. pT, z, xi
2334   TH1* hSpecPt_RecPrimaries = 0x0;
2335   TH1* hSpecPt_GenPrimaries = 0x0;
2336   TH1D* hEfficiencyAllPt = 0x0;
2337
2338   TH2* hSpecID2Pt_RecPrimaries = 0x0;
2339   TH2* hSpecID2Pt_GenPrimaries = 0x0;
2340   TH2D* hEffID2Pt = 0x0;
2341
2342   TH1D* hSpecIDPt_GenPrimaries[AliPID::kSPECIES] = { 0x0, };
2343   TH1D* hSpecIDPt_RecPrimaries[AliPID::kSPECIES] = { 0x0, };
2344   TH1D* hEfficiencyPt[AliPID::kSPECIES] = { 0x0, };
2345   TH1D* hEfficiencyToPiRatioPt[AliPID::kSPECIES] = { 0x0, };
2346
2347   TCanvas *cSpecPt = 0x0, *cEffPt = 0x0, *cEffToPiRatioPt = 0x0;
2348
2349
2350   TH1* hSpecZ_RecPrimaries = 0x0;
2351   TH1* hSpecZ_GenPrimaries = 0x0;
2352   TH1D* hEfficiencyAllZ = 0x0;
2353
2354   TH2* hSpecID2Z_RecPrimaries = 0x0;
2355   TH2* hSpecID2Z_GenPrimaries = 0x0;
2356   TH2D* hEffID2Z = 0x0;
2357
2358   TH1D* hSpecIDZ_GenPrimaries[AliPID::kSPECIES] = { 0x0, };
2359   TH1D* hSpecIDZ_RecPrimaries[AliPID::kSPECIES] = { 0x0, };
2360   TH1D* hEfficiencyZ[AliPID::kSPECIES] = { 0x0, };
2361   TH1D* hEfficiencyToPiRatioZ[AliPID::kSPECIES] = { 0x0, };
2362
2363   TCanvas *cSpecZ = 0x0, *cEffZ = 0x0, *cEffToPiRatioZ = 0x0;
2364
2365
2366   TH1* hSpecXi_RecPrimaries = 0x0;
2367   TH1* hSpecXi_GenPrimaries = 0x0;
2368   TH1D* hEfficiencyAllXi = 0x0;
2369
2370   TH2* hSpecID2Xi_RecPrimaries = 0x0;
2371   TH2* hSpecID2Xi_GenPrimaries = 0x0;
2372   TH2D* hEffID2Xi = 0x0;
2373
2374   TH1D* hSpecIDXi_GenPrimaries[AliPID::kSPECIES] = { 0x0, };
2375   TH1D* hSpecIDXi_RecPrimaries[AliPID::kSPECIES] = { 0x0, };
2376   TH1D* hEfficiencyXi[AliPID::kSPECIES] = { 0x0, };
2377   TH1D* hEfficiencyToPiRatioXi[AliPID::kSPECIES] = { 0x0, };
2378
2379   TCanvas *cSpecXi = 0x0, *cEffXi = 0x0, *cEffToPiRatioXi = 0x0;
2380
2381   TString titleYaxis[kNtypes];
2382   titleYaxis[kTrackPt] = "1/N_{Jets} dN/dP_{T} (GeV/c)^{-1}";
2383   titleYaxis[kZ] = "1/N_{Jets} dN/dz";
2384   titleYaxis[kXi] = "1/N_{Jets} dN/d#xi";
2385
2386   for (Int_t type = 0; type < kNtypes; type++) {
2387     TH1** hSpec_RecPrimaries;
2388     TH1** hSpec_GenPrimaries;
2389     TH1D** hEfficiencyAll;
2390
2391     TH2** hSpecID2_RecPrimaries;
2392     TH2** hSpecID2_GenPrimaries;
2393     TH2D** hEffID2;
2394
2395     TH1D** hSpecID_GenPrimaries;
2396     TH1D** hSpecID_RecPrimaries;
2397     TH1D** hEfficiency;
2398     TH1D** hEfficiencyToPiRatio;
2399
2400     TCanvas **cSpec, **cEff, **cEffToPiRatio;
2401
2402     Int_t iDimProj = -1;
2403
2404
2405     if (type == kTrackPt) {
2406       hSpec_RecPrimaries = &hSpecPt_RecPrimaries;
2407       hSpec_GenPrimaries = &hSpecPt_GenPrimaries;
2408       hEfficiencyAll = &hEfficiencyAllPt;
2409
2410       hSpecID2_RecPrimaries = &hSpecID2Pt_RecPrimaries;
2411       hSpecID2_GenPrimaries = &hSpecID2Pt_GenPrimaries;
2412       hEffID2 = &hEffID2Pt;
2413
2414       hSpecID_RecPrimaries = &hSpecIDPt_RecPrimaries[0];
2415       hSpecID_GenPrimaries = &hSpecIDPt_GenPrimaries[0];
2416       hEfficiency = &hEfficiencyPt[0];
2417       hEfficiencyToPiRatio = &hEfficiencyToPiRatioPt[0];
2418
2419       cSpec = &cSpecPt;
2420       cEff = &cEffPt;
2421       cEffToPiRatio = &cEffToPiRatioPt;
2422
2423       iDimProj = iPt;
2424     }
2425     else if (type == kZ) {
2426       hSpec_RecPrimaries = &hSpecZ_RecPrimaries;
2427       hSpec_GenPrimaries = &hSpecZ_GenPrimaries;
2428       hEfficiencyAll = &hEfficiencyAllZ;
2429
2430       hSpecID2_RecPrimaries = &hSpecID2Z_RecPrimaries;
2431       hSpecID2_GenPrimaries = &hSpecID2Z_GenPrimaries;
2432       hEffID2 = &hEffID2Z;
2433
2434       hSpecID_RecPrimaries = &hSpecIDZ_RecPrimaries[0];
2435       hSpecID_GenPrimaries = &hSpecIDZ_GenPrimaries[0];
2436       hEfficiency = &hEfficiencyZ[0];
2437       hEfficiencyToPiRatio = &hEfficiencyToPiRatioZ[0];
2438
2439       cSpec = &cSpecZ;
2440       cEff = &cEffZ;
2441       cEffToPiRatio = &cEffToPiRatioZ;
2442
2443       iDimProj = iZ;
2444     }
2445     else if (type == kXi) {
2446       hSpec_RecPrimaries = &hSpecXi_RecPrimaries;
2447       hSpec_GenPrimaries = &hSpecXi_GenPrimaries;
2448       hEfficiencyAll = &hEfficiencyAllXi;
2449
2450       hSpecID2_RecPrimaries = &hSpecID2Xi_RecPrimaries;
2451       hSpecID2_GenPrimaries = &hSpecID2Xi_GenPrimaries;
2452       hEffID2 = &hEffID2Xi;
2453
2454       hSpecID_RecPrimaries = &hSpecIDXi_RecPrimaries[0];
2455       hSpecID_GenPrimaries = &hSpecIDXi_GenPrimaries[0];
2456       hEfficiency = &hEfficiencyXi[0];
2457       hEfficiencyToPiRatio = &hEfficiencyToPiRatioXi[0];
2458
2459       cSpec = &cSpecXi;
2460       cEff = &cEffXi;
2461       cEffToPiRatio = &cEffToPiRatioXi;
2462
2463       iDimProj = iXi;
2464     }
2465     else
2466       continue;
2467
2468
2469     // Divide steps: recStepEff (=kStepRecWithRecCutsMeasuredObsPrimaries) / genStepEff (=kStepGenWithGenCuts)
2470     *hSpec_RecPrimaries = dataRebinned->GetGrid(recStepEff)->Project(iDimProj);
2471     setupHist(*hSpec_RecPrimaries, Form("hSpec_%s_RecPrimaries", titles[type].Data()), "All, rec", "",
2472               titleYaxis[type].Data(), kBlack);
2473     normaliseHist(*hSpec_RecPrimaries, nJetsRec > 0 ? 1. / nJetsRec : 0.);
2474     (*hSpec_RecPrimaries)->SetLineWidth(2.);
2475     (*hSpec_RecPrimaries)->SetMarkerStyle(20);
2476
2477     *hSpec_GenPrimaries = dataRebinned->GetGrid(genStepEff)->Project(iDimProj);
2478     setupHist(*hSpec_GenPrimaries, Form("hSpec_%s_GenPrimaries", titles[type].Data()), "All, gen", "",
2479               titleYaxis[type].Data(), kBlack);
2480     normaliseHist(*hSpec_GenPrimaries, nJetsGen > 0 ? 1. / nJetsGen : 0.);
2481     (*hSpec_GenPrimaries)->SetLineWidth(2.);
2482     (*hSpec_GenPrimaries)->SetLineStyle(2);
2483
2484     *hEfficiencyAll = new TH1D(*(TH1D*)(*hSpec_RecPrimaries));
2485     (*hEfficiencyAll)->SetName(Form("hEfficiencyAll_%s", titles[type].Data()));
2486     (*hEfficiencyAll)->GetYaxis()->SetTitle("");
2487     (*hEfficiencyAll)->SetTitle("All");
2488     (*hEfficiencyAll)->SetLineWidth(2.0);
2489     (*hEfficiencyAll)->SetStats(kFALSE);
2490     (*hEfficiencyAll)->GetYaxis()->SetTitle("Efficiency x Acceptance x pT Resolution");
2491     (*hEfficiencyAll)->GetYaxis()->SetTitleSize(0.05);
2492     (*hEfficiencyAll)->Divide(*hSpec_RecPrimaries, *hSpec_GenPrimaries, 1., 1., "B");
2493     (*hEfficiencyAll)->GetYaxis()->SetRangeUser(0., 2.0);
2494
2495     *hSpecID2_RecPrimaries = (TH2*)dataRebinned->GetGrid(recStepEff)->Project(iDimProj, iMCid);
2496     (*hSpecID2_RecPrimaries)->SetName(Form("hSpecID2_RecPrimaries_%s", titles[type].Data()));
2497     (*hSpecID2_RecPrimaries)->Scale(nJetsRec > 0 ? 1. / nJetsRec : 0.);
2498
2499     *hSpecID2_GenPrimaries = (TH2*)dataRebinned->GetGrid(genStepEff)->Project(iDimProj, iMCid);
2500     (*hSpecID2_GenPrimaries)->SetName(Form("hSpecID2_GenPrimaries_%s", titles[type].Data()));
2501     (*hSpecID2_GenPrimaries)->Scale(nJetsGen > 0 ? 1. / nJetsGen : 0.);
2502
2503     *hEffID2 = new TH2D(*(TH2D*)(*hSpecID2_RecPrimaries));
2504     (*hEffID2)->SetName(Form("hEff_%s", titles[type].Data()));
2505     (*hEffID2)->SetTitle("All");
2506     (*hEffID2)->Divide(*hSpecID2_RecPrimaries, *hSpecID2_GenPrimaries, 1., 1., "B");
2507
2508     // Get the efficiencies vs. pT,z,xi for each species
2509     for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
2510       hEfficiency[species] = ((TH2D*)*hEffID2)->ProjectionX(Form("hEfficiency_%s_%s", titles[type].Data(),
2511                                                                  AliPID::ParticleShortName(species)),
2512                                                             species + 1, species + 1, "e");
2513       setupHist(hEfficiency[species], "", Form("%s", AliPID::ParticleLatexName(species)), "", "",
2514                 getLineColorAliPID(species));
2515       hEfficiency[species]->SetLineWidth(2.);
2516       hEfficiency[species]->GetYaxis()->SetRangeUser(0., 2.0);
2517       hEfficiency[species]->SetStats(kFALSE);
2518       hEfficiency[species]->GetYaxis()->SetTitle("Efficiency x Acceptance x pT Resolution");
2519       hEfficiency[species]->GetYaxis()->SetTitleSize(0.05);
2520
2521       hSpecID_GenPrimaries[species] = ((TH2D*)*hSpecID2_GenPrimaries)->ProjectionX(Form("hSpecID_%s_gen_%s",
2522                                                                                         titles[type].Data(),
2523                                                                                         AliPID::ParticleShortName(species)),
2524                                                                                         species + 1, species + 1, "e");
2525       setupHist(hSpecID_GenPrimaries[species], "", Form("%s, gen", AliPID::ParticleLatexName(species)), "",
2526                 titleYaxis[type].Data(), getLineColorAliPID(species));
2527       normaliseHist(hSpecID_GenPrimaries[species], 1.);
2528       hSpecID_GenPrimaries[species]->SetLineWidth(2.);
2529       hSpecID_GenPrimaries[species]->SetLineStyle(2.);
2530       hSpecID_GenPrimaries[species]->SetStats(kFALSE);
2531
2532
2533       hSpecID_RecPrimaries[species] = ((TH2D*)*hSpecID2_RecPrimaries)->ProjectionX(Form("hSpecID_%s_rec_%s",
2534                                                                                         titles[type].Data(),
2535                                                                                         AliPID::ParticleShortName(species)),
2536                                                                                         species + 1, species + 1, "e");
2537       setupHist(hSpecID_RecPrimaries[species], "", Form("%s, rec", AliPID::ParticleLatexName(species)), "",
2538                 titleYaxis[type].Data(), getLineColorAliPID(species));
2539       normaliseHist(hSpecID_RecPrimaries[species], 1.);
2540       hSpecID_RecPrimaries[species]->SetLineWidth(2.);
2541       hSpecID_RecPrimaries[species]->SetMarkerStyle(20);
2542       hSpecID_RecPrimaries[species]->SetStats(kFALSE);
2543     }
2544
2545     *cSpec = new TCanvas(Form("cSpec_%s", titles[type].Data()), "Spectra for different species", 0, 300, 900, 900);
2546     (*cSpec)->cd();
2547     (*cSpec)->SetGridx(1);
2548     (*cSpec)->SetGridy(1);
2549     (*cSpec)->SetLogy(1);
2550
2551     if (type == kTrackPt)
2552       (*hSpec_GenPrimaries)->GetYaxis()->SetRangeUser(1e-8, 1e2);
2553     else if (type == kZ)
2554       (*hSpec_GenPrimaries)->GetYaxis()->SetRangeUser(1e-5, 1e3);
2555     else if (type == kXi)
2556       (*hSpec_GenPrimaries)->GetYaxis()->SetRangeUser(1e-5, 1e3);
2557
2558
2559     if (type == kTrackPt)
2560       (*hSpec_GenPrimaries)->GetXaxis()->SetRangeUser(0.15, upperTrackPt);
2561     (*hSpec_GenPrimaries)->Draw("E1");
2562
2563     if (type == kTrackPt)
2564       (*hSpec_RecPrimaries)->GetXaxis()->SetRangeUser(0.15, upperTrackPt);
2565     (*hSpec_RecPrimaries)->Draw("E1 same");
2566
2567     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2568       if (type == kTrackPt) {
2569         hSpecID_GenPrimaries[i]->GetXaxis()->SetRangeUser(0.15, upperTrackPt);
2570         hSpecID_RecPrimaries[i]->GetXaxis()->SetRangeUser(0.15, upperTrackPt);
2571       }
2572
2573       hSpecIDPt_GenPrimaries[i]->Draw("E1 same");
2574       hSpecIDPt_RecPrimaries[i]->Draw("E1 same");
2575     }
2576     TLegend* leg = (*cSpec)->BuildLegend();
2577     leg->SetFillColor(kWhite);
2578     leg->SetNColumns(2);
2579
2580     ClearTitleFromHistoInCanvas(*cSpec);
2581
2582
2583     *cEff = new TCanvas(Form("cEff_%s", titles[type].Data()), "Efficiency x Acceptance x pT Resolution for different species",
2584                         0, 300, 900, 900);
2585     (*cEff)->cd();
2586     (*cEff)->SetGridx(1);
2587     (*cEff)->SetGridy(1);
2588
2589     if (type == kTrackPt)
2590       (*hEfficiencyAll)->GetXaxis()->SetRangeUser(0.15, upperTrackPt);
2591     (*hEfficiencyAll)->Draw("E1");
2592
2593     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2594       if (type == kTrackPt)
2595         hEfficiency[i]->GetXaxis()->SetRangeUser(0.15, upperTrackPt);
2596       hEfficiency[i]->Draw("E1 same");
2597     }
2598     (*cEff)->BuildLegend()->SetFillColor(kWhite);
2599
2600     ClearTitleFromHistoInCanvas(*cEff);
2601
2602
2603     // To-pi ratios
2604     for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
2605       if (species == AliPID::kPion)
2606         continue; // Do not consider pion-to-pion ratio
2607
2608       hEfficiencyToPiRatio[species] = new TH1D(*hEfficiency[species]);
2609       hEfficiencyToPiRatio[species]->Reset();
2610       hEfficiencyToPiRatio[species]->SetName(Form("hEfficiency_ToPionRatio_%s_%s", titles[type].Data(),
2611                                                   AliPID::ParticleShortName(species)));
2612
2613       setupHist(hEfficiencyToPiRatio[species], "", Form("%s/#pi", AliPID::ParticleLatexName(species)), "", "",
2614                 getLineColorAliPID(species));
2615       hEfficiencyToPiRatio[species]->SetLineWidth(2.);
2616       hEfficiencyToPiRatio[species]->SetStats(kFALSE);
2617       hEfficiencyToPiRatio[species]->GetYaxis()->SetTitle("Efficiency x Acceptance x pT Resolution");
2618       hEfficiencyToPiRatio[species]->GetYaxis()->SetTitleSize(0.05);
2619
2620       // Samples for different species are independent, so just divide correction factors
2621       hEfficiencyToPiRatio[species]->Divide(hEfficiency[species], hEfficiency[AliPID::kPion], 1., 1., "");
2622       hEfficiencyToPiRatio[species]->GetYaxis()->SetRangeUser(0., 2.0);
2623     }
2624
2625     *cEffToPiRatio = new TCanvas(Form("cEffToPiRatio_%s", titles[type].Data()),
2626                                  "Efficiency x Acceptance x pT Resolution for different to-#pi-ratios",
2627                                  100, 10, 1200, 800);
2628     (*cEffToPiRatio)->cd();
2629     (*cEffToPiRatio)->SetGridx(1);
2630     (*cEffToPiRatio)->SetGridy(1);
2631
2632     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2633       if (i == AliPID::kPion)
2634         continue;
2635
2636       if (type == kTrackPt)
2637         hEfficiencyToPiRatio[i]->GetXaxis()->SetRangeUser(0.15, upperTrackPt);
2638
2639       hEfficiencyToPiRatio[i]->Draw(Form("E1%s", i == 0 ? "" : " same"));
2640     }
2641     (*cEffToPiRatio)->BuildLegend()->SetFillColor(kWhite);
2642
2643     ClearTitleFromHistoInCanvas(*cEffToPiRatio);
2644   }
2645
2646   // Save results to file
2647   TString chargeString = "";
2648   if (chargeMode == kPosCharge)
2649     chargeString = "_posCharge";
2650   else if (chargeMode == kNegCharge)
2651     chargeString = "_negCharge";
2652
2653   TString saveFileName = pathNameData;
2654   saveFileName.ReplaceAll(Form("%s/", pathData.Data()), "");
2655   saveFileName.Prepend("output_EfficiencyCorrection_");
2656   saveFileName.ReplaceAll(".root", Form("_jetPt%.1f_%.1f%s.root", actualLowerJetPt, actualUpperJetPt, chargeString.Data()));
2657
2658   TString saveFilePathName = Form("%s/%s", pathData.Data(), saveFileName.Data());
2659   TFile* saveFile = TFile::Open(saveFilePathName.Data(), "RECREATE");
2660
2661   if (!saveFile) {
2662     printf("Could not save results!\n");
2663
2664     return -1;
2665   }
2666
2667   saveFile->cd();
2668
2669   if (cSecPt)
2670     cSecPt->Write();
2671
2672   if (cSecSSPt)
2673     cSecSSPt->Write();
2674
2675   if (hSecAllPt)
2676       hSecAllPt->Write();
2677
2678   if (hSecSSAllPt)
2679       hSecSSAllPt->Write();
2680
2681   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2682     if (hSecPt[i])
2683       hSecPt[i]->Write();
2684
2685     if (hSecSSPt[i])
2686       hSecSSPt[i]->Write();
2687   }
2688
2689   if (cSecToPiRatioPt)
2690     cSecToPiRatioPt->Write();
2691
2692   if (cSecSSToPiRatioPt)
2693     cSecSSToPiRatioPt->Write();
2694
2695   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2696     if (hSecToPiRatioPt[i])
2697       hSecToPiRatioPt[i]->Write();
2698
2699     if (hSecSSToPiRatioPt[i])
2700       hSecSSToPiRatioPt[i]->Write();
2701   }
2702
2703
2704
2705   if (cSecZ)
2706     cSecZ->Write();
2707
2708   if (cSecSSZ)
2709     cSecSSZ->Write();
2710
2711   if (hSecAllZ)
2712       hSecAllZ->Write();
2713
2714   if (hSecSSAllZ)
2715       hSecSSAllZ->Write();
2716
2717   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2718     if (hSecZ[i])
2719       hSecZ[i]->Write();
2720
2721     if (hSecSSZ[i])
2722       hSecSSZ[i]->Write();
2723   }
2724
2725   if (cSecToPiRatioZ)
2726     cSecToPiRatioZ->Write();
2727
2728   if (cSecSSToPiRatioZ)
2729     cSecSSToPiRatioZ->Write();
2730
2731   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2732     if (hSecToPiRatioZ[i])
2733       hSecToPiRatioZ[i]->Write();
2734
2735     if (hSecSSToPiRatioZ[i])
2736       hSecSSToPiRatioZ[i]->Write();
2737   }
2738
2739
2740
2741   if (cSecXi)
2742     cSecXi->Write();
2743
2744   if (cSecSSXi)
2745     cSecSSXi->Write();
2746
2747   if (hSecAllXi)
2748       hSecAllXi->Write();
2749
2750   if (hSecSSAllXi)
2751       hSecSSAllXi->Write();
2752
2753   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2754     if (hSecXi[i])
2755       hSecXi[i]->Write();
2756
2757     if (hSecSSXi[i])
2758       hSecSSXi[i]->Write();
2759   }
2760
2761   if (cSecToPiRatioXi)
2762     cSecToPiRatioXi->Write();
2763
2764   if (cSecSSToPiRatioXi)
2765     cSecSSToPiRatioXi->Write();
2766
2767   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2768     if (hSecToPiRatioXi[i])
2769       hSecToPiRatioXi[i]->Write();
2770
2771     if (hSecSSToPiRatioXi[i])
2772       hSecSSToPiRatioXi[i]->Write();
2773   }
2774
2775
2776
2777   if (cSpecPt)
2778     cSpecPt->Write();
2779
2780   if (hSpecPt_RecPrimaries)
2781     hSpecPt_RecPrimaries->Write();
2782
2783   if (hSpecPt_GenPrimaries)
2784     hSpecPt_GenPrimaries->Write();
2785
2786   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2787     if (hSpecIDPt_GenPrimaries[i])
2788       hSpecIDPt_GenPrimaries[i]->Write();
2789
2790     if (hSpecIDPt_RecPrimaries[i])
2791       hSpecIDPt_RecPrimaries[i]->Write();
2792   }
2793
2794
2795   if (cEffPt)
2796     cEffPt->Write();
2797
2798   if (hEfficiencyAllPt)
2799       hEfficiencyAllPt->Write();
2800
2801   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2802     if (hEfficiencyPt[i])
2803       hEfficiencyPt[i]->Write();
2804   }
2805
2806   if (cEffToPiRatioPt)
2807     cEffToPiRatioPt->Write();
2808
2809   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2810     if (hEfficiencyToPiRatioPt[i])
2811       hEfficiencyToPiRatioPt[i]->Write();
2812   }
2813
2814
2815
2816   if (cSpecZ)
2817     cSpecZ->Write();
2818
2819   if (hSpecZ_RecPrimaries)
2820     hSpecZ_RecPrimaries->Write();
2821
2822   if (hSpecZ_GenPrimaries)
2823     hSpecZ_GenPrimaries->Write();
2824
2825   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2826     if (hSpecIDZ_GenPrimaries[i])
2827       hSpecIDZ_GenPrimaries[i]->Write();
2828
2829     if (hSpecIDZ_RecPrimaries[i])
2830       hSpecIDZ_RecPrimaries[i]->Write();
2831   }
2832
2833
2834   if (cEffZ)
2835     cEffZ->Write();
2836
2837   if (hEfficiencyAllZ)
2838       hEfficiencyAllZ->Write();
2839
2840   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2841     if (hEfficiencyZ[i])
2842       hEfficiencyZ[i]->Write();
2843   }
2844
2845   if (cEffToPiRatioZ)
2846     cEffToPiRatioZ->Write();
2847
2848   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2849     if (hEfficiencyToPiRatioZ[i])
2850       hEfficiencyToPiRatioZ[i]->Write();
2851   }
2852
2853
2854
2855   if (cSpecXi)
2856     cSpecXi->Write();
2857
2858   if (hSpecXi_RecPrimaries)
2859     hSpecXi_RecPrimaries->Write();
2860
2861   if (hSpecXi_GenPrimaries)
2862     hSpecXi_GenPrimaries->Write();
2863
2864   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2865     if (hSpecIDXi_GenPrimaries[i])
2866       hSpecIDXi_GenPrimaries[i]->Write();
2867
2868     if (hSpecIDXi_RecPrimaries[i])
2869       hSpecIDXi_RecPrimaries[i]->Write();
2870   }
2871
2872
2873   if (cEffXi)
2874     cEffXi->Write();
2875
2876   if (hEfficiencyAllXi)
2877       hEfficiencyAllXi->Write();
2878
2879   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2880     if (hEfficiencyXi[i])
2881       hEfficiencyXi[i]->Write();
2882   }
2883
2884   if (cEffToPiRatioXi)
2885     cEffToPiRatioXi->Write();
2886
2887   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2888     if (hEfficiencyToPiRatioXi[i])
2889       hEfficiencyToPiRatioXi[i]->Write();
2890   }
2891
2892
2893
2894   // Correct the yields with the efficiencies and primary fractions
2895   for (Int_t type = 0; type < kNtypes; type++) {
2896     TH1D* hYieldCorrected[AliPID::kSPECIES];
2897     TH1D* hYieldCorrectedSysError[AliPID::kSPECIES];
2898
2899     TH1D** hYield;
2900     TH1D** hYieldSysError;
2901     TH1D** hMCgenPrimYield;
2902     TH1D** hSecSS;
2903     TH1D** hSec;
2904     TH1D** hEfficiency;
2905
2906     if (type == kTrackPt) {
2907       hYield = &hYieldPt[0];
2908       hYieldSysError = &hYieldPtSysError[0];
2909       hMCgenPrimYield = &hMCgenPrimYieldPt[0];
2910       hSecSS = &hSecSSPt[0];
2911       hSec = &hSecPt[0];
2912       hEfficiency = &hEfficiencyPt[0];
2913     }
2914     else if (type == kZ) {
2915       hYield = &hYieldZ[0];
2916       hYieldSysError = &hYieldZSysError[0];
2917       hMCgenPrimYield = &hMCgenPrimYieldZ[0];
2918       hSecSS = &hSecSSZ[0];
2919       hSec = &hSecZ[0];
2920       hEfficiency = &hEfficiencyZ[0];
2921     }
2922     else if (type == kXi) {
2923       hYield = &hYieldXi[0];
2924       hYieldSysError = &hYieldXiSysError[0];
2925       hMCgenPrimYield = &hMCgenPrimYieldXi[0];
2926       hSecSS = &hSecSSXi[0];
2927       hSec = &hSecXi[0];
2928       hEfficiency = &hEfficiencyXi[0];
2929     }
2930     else
2931       continue;
2932
2933
2934     // NOTE: The yields have already been normalised to nJets in extractFFs. However, the binning in jetPt might be different now.
2935     // Thus, in order to get the correct result, the old normalisation (just the same binnnig as in the nJet histo) has to be
2936     // undone (was done above) and then it has to be normalised to the number of jets in the new binning again (now).
2937     // Also NOTE: Histos are already normalised to bin width via extractFFs!
2938     const Double_t nJetsGenData = hNjetsGenData ? hNjetsGenData->Integral(lowerCentralityBinLimitData, upperCentralityBinLimitData,
2939                                                                           lowerJetPtBinLimitData, upperJetPtBinLimitData) : 1.;
2940     const Double_t nJetsRecData = hNjetsRecData ? hNjetsRecData->Integral(lowerCentralityBinLimitData, upperCentralityBinLimitData,
2941                                                                           lowerJetPtBinLimitData, upperJetPtBinLimitData) : 1.;
2942
2943     for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
2944       hYield[species]->Scale((nJetsRecData > 0) ? 1. / nJetsRecData : 0.);
2945       hYieldSysError[species]->Scale((nJetsRecData > 0) ? 1. / nJetsRecData : 0.);
2946
2947       hYield[species]->SetStats(kFALSE);
2948       hYieldSysError[species]->SetStats(kFALSE);
2949
2950       if (hMCgenPrimYield[species]) {
2951         hMCgenPrimYield[species]->Scale((nJetsGenData > 0) ? 1. / nJetsGenData : 0.);
2952         hMCgenPrimYield[species]->SetStats(kFALSE);
2953       }
2954
2955       hYield[species]->GetYaxis()->SetTitle(titleYaxis[type].Data());
2956       if (hMCgenPrimYield[species]) {
2957         hMCgenPrimYield[species]->GetYaxis()->SetTitle(titleYaxis[type].Data());
2958         hMCgenPrimYield[species]->SetTitle(Form("%s (MC)", AliPID::ParticleLatexName(species)));
2959         hMCgenPrimYield[species]->SetMarkerStyle(24);
2960         hMCgenPrimYield[species]->SetLineColor(getLineColorAliPID(species));
2961         hMCgenPrimYield[species]->SetMarkerColor(getLineColorAliPID(species));
2962       }
2963
2964       hYield[species]->SetLineColor(getLineColorAliPID(species));
2965       hYield[species]->SetMarkerColor(getLineColorAliPID(species));
2966       hYield[species]->SetMarkerStyle(20);
2967
2968       hYieldSysError[species]->SetLineColor(getLineColorAliPID(species));
2969       hYieldSysError[species]->SetMarkerColor(getLineColorAliPID(species));
2970       hYieldSysError[species]->SetMarkerStyle(20);
2971       hYieldSysError[species]->SetFillStyle(0);
2972
2973       hYieldCorrected[species] = new TH1D(*hYield[species]);
2974       hYieldCorrected[species]->SetName(Form("%s_corrected", hYield[species]->GetName()));
2975       hYieldCorrected[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species)));
2976
2977       hYieldCorrectedSysError[species] = new TH1D(*hYieldSysError[species]);
2978       hYieldCorrectedSysError[species]->SetName(Form("%s_corrected", hYieldSysError[species]->GetName()));
2979       hYieldCorrectedSysError[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species)));
2980
2981       // Correction factor histos can have different binning than data histos (constant factor above some threshold)
2982       // -> Need special functions to multiply and divide such histos
2983       multiplyHistsDifferentBinning(hYieldCorrected[species], scaleStrangeness ? hSecSS[species] : hSec[species], 1., 1.);
2984       divideHistsDifferentBinning(hYieldCorrected[species], hEfficiency[species], 1., 1.);
2985
2986       multiplyHistsDifferentBinning(hYieldCorrectedSysError[species], scaleStrangeness ? hSecSS[species] : hSec[species], 1., 1.,
2987                                     kTRUE);
2988       divideHistsDifferentBinning(hYieldCorrectedSysError[species], hEfficiency[species], 1., 1., kTRUE);
2989
2990       //hYieldCorrected[species]->Multiply(hYieldCorrected[species], scaleStrangeness ? hSecSS[species] : hSec[species], 1., 1.,
2991       //"");
2992       //hYieldCorrected[species]->Divide(hYieldCorrected[species], hEfficiency[species], 1., 1., "");
2993     }
2994
2995     // Calculate the total corrected yield. The original total yield had no error (just the number of detected tracks in a pT bin),
2996     // but due to the correction there is some error for the total yield. Also the error of the fractions introduces uncertainties
2997     // for the yields of individual species
2998     TH1D* hYieldCorrectedTotal = new TH1D(*hYieldCorrected[0]);
2999     hYieldCorrectedTotal->SetLineColor(kBlack);
3000     hYieldCorrectedTotal->SetMarkerColor(kBlack);
3001     hYieldCorrectedTotal->SetMarkerStyle(20);
3002     hYieldCorrectedTotal->SetName(Form("hYieldCorrectedTotal_%s", titles[type].Data()));
3003     hYieldCorrectedTotal->SetTitle("Total");
3004     //hYieldCorrectedTotal->SetTitle("Total, secondary and efficiency x acceptance x pT resolution corrected");
3005
3006     //TODO doesn't it require the covariance matrix here to get the error of the corrected total yield?
3007     // Because before the total count was taken into account just by the fit. Now the fractions and their errors
3008     // influence the MC correction (and its errors), so that there is some error on the corrected total yield.
3009     // However, the errors of the contributing yields are correlated!
3010     for (Int_t i = 1; i < AliPID::kSPECIES; i++)
3011       hYieldCorrectedTotal->Add(hYieldCorrected[i], 1.);
3012
3013     // Calculate the corrected fractions
3014     TH1D* hFractionCorrected[AliPID::kSPECIES];
3015     TH1D* hFractionCorrectedSysError[AliPID::kSPECIES];
3016
3017     for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
3018       hFractionCorrected[species] = new TH1D(*hYield[species]);
3019       TString oldName = hYield[species]->GetName();
3020       TString newName = oldName.ReplaceAll("Yield", "Fraction");
3021       hFractionCorrected[species]->SetName(Form("%s_corrected", newName.Data()));
3022       hFractionCorrected[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species)));
3023       hFractionCorrected[species]->GetYaxis()->SetTitle("Corrected Fraction");
3024       hFractionCorrected[species]->GetXaxis()->SetMoreLogLabels(kTRUE);
3025       hFractionCorrected[species]->GetXaxis()->SetNoExponent(kTRUE);
3026       hFractionCorrected[species]->SetFillStyle(0);
3027
3028
3029       // Binomial error as for efficiencies (numerator and denominator are statistically not independent) for correct error calculation
3030       // (numerator is a "selected" subset of the denominator). It doesn't matter that the error of the histos is not just
3031       // "sqrt(content)" because the error formula also works for weighted histograms (which means that the error can be more
3032       // or less anything)
3033       hFractionCorrected[species]->Divide(hYieldCorrected[species], hYieldCorrectedTotal, 1., 1., "B");
3034
3035
3036
3037       //  The systematic errors just need to be scaled in the same way as the fractions.
3038       // So, just take the ratio to the uncorrected fraction and scale the sys. error accordingly
3039       // or, in this case, just divide by the same total yield as for yield -> fractions
3040       hFractionCorrectedSysError[species] = new TH1D(*hFractionCorrected[species]);
3041       hFractionCorrectedSysError[species]->SetName(Form("%s_sysError", hFractionCorrected[species]->GetName()));
3042       hFractionCorrectedSysError[species]->SetTitle(Form("%s (sys. error)", hFractionCorrected[species]->GetTitle()));
3043
3044       for (Int_t binX = 1; binX <= hFractionCorrectedSysError[species]->GetNbinsX(); binX++) {
3045         const Double_t corrTotalYield = hYieldCorrectedTotal->GetBinContent(binX);
3046         const Double_t scaleFactor = corrTotalYield > 0 ? 1.0 / corrTotalYield : 1.;
3047         hFractionCorrectedSysError[species]->SetBinError(binX,   hYieldCorrectedSysError[species]->GetBinError(binX) * scaleFactor);
3048       }
3049     }
3050
3051     // If MC is available, calculate the generated fractions
3052     TH1D* hMCgenPrimYieldTotal = 0x0;
3053     TH1D* hMCgenPrimFraction[AliPID::kSPECIES];
3054     for (Int_t i = 0; i < AliPID::kSPECIES; i++)
3055       hMCgenPrimFraction[i] = 0x0;
3056
3057     if (numMCgenPrimYieldHistsFound > 0) {
3058       hMCgenPrimYieldTotal = new TH1D(*hMCgenPrimYield[0]);
3059       hMCgenPrimYieldTotal->SetLineColor(kBlack);
3060       hMCgenPrimYieldTotal->SetMarkerColor(kBlack);
3061       hMCgenPrimYieldTotal->SetMarkerStyle(24);
3062       hMCgenPrimYieldTotal->SetName(Form("hMCgenPrimYieldTotal%s", titles[type].Data()));
3063       hMCgenPrimYieldTotal->SetTitle("Total, MC truth");
3064       //hMCgenPrimYieldTotal->SetTitle("Total generated primary yield (MC truth)");
3065
3066       for (Int_t i = 1; i < AliPID::kSPECIES; i++)
3067         hMCgenPrimYieldTotal->Add(hMCgenPrimYield[i], 1.);
3068
3069       // Calculate the MC fractions
3070       for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
3071         hMCgenPrimYield[species]->SetLineColor(getLineColorAliPID(species));
3072         hMCgenPrimYield[species]->SetMarkerColor(getLineColorAliPID(species));
3073         hMCgenPrimYield[species]->SetMarkerStyle(24);
3074
3075         hMCgenPrimYield[species]->SetTitle(Form("%s, MC truth", AliPID::ParticleLatexName(species)));
3076
3077         hMCgenPrimFraction[species] = new TH1D(*hMCgenPrimYield[species]);
3078         TString oldName = hMCgenPrimYield[species]->GetName();
3079         TString newName = oldName.ReplaceAll("Yield", "Fraction");
3080         hMCgenPrimFraction[species]->SetName(newName.Data());
3081
3082         // Binomial error as for efficiencies (numerator and denominator are statistically not independent) for correct error calculation
3083         // (numerator is a "selected" subset of the denominator).
3084         hMCgenPrimFraction[species]->Divide(hMCgenPrimFraction[species], hMCgenPrimYieldTotal, 1., 1., "B");
3085       }
3086     }
3087
3088     TCanvas* cCorrData = new TCanvas(Form("cCorrData_%s", titles[type].Data()), "Corrected data", 0, 300, 900, 900);
3089     cCorrData->Divide(2, 1, 0., 0.01);
3090     cCorrData->GetPad(1)->SetLogx(1);
3091     cCorrData->GetPad(1)->SetLogy(1);
3092     cCorrData->GetPad(2)->SetLogx(1);
3093     cCorrData->GetPad(2)->SetLogy(1);
3094
3095     cCorrData->GetPad(1)->SetRightMargin(0.001);
3096     cCorrData->GetPad(2)->SetRightMargin(0.001);
3097
3098     cCorrData->GetPad(1)->SetLeftMargin(0.2);
3099     cCorrData->GetPad(2)->SetLeftMargin(0.2);
3100
3101     cCorrData->cd(1); // uncorrected
3102     hYield[AliPID::kPion]->GetYaxis()->SetTitleOffset(1.4);
3103     hYield[AliPID::kPion]->Draw("E1");
3104
3105     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
3106       hYieldSysError[i]->Draw("E2 same");
3107
3108       if (i == AliPID::kPion) continue;
3109       hYield[i]->Draw("E1 same");
3110     }
3111
3112
3113     ClearTitleFromHistoInCanvas(cCorrData, 1);
3114
3115     Double_t maximum = hYieldCorrectedTotal->GetBinContent(hYieldCorrectedTotal->GetMaximumBin()) * 10.;
3116     Double_t minimum = maximum;
3117     for (Int_t spec = 0; spec < AliPID::kSPECIES; spec++) {
3118       Double_t temp = hYieldCorrected[spec]->GetBinContent(hYieldCorrected[spec]->FindLastBinAbove(0.)) * 0.1;
3119       if (temp > 0)
3120         minimum = temp;
3121     }
3122
3123
3124     cCorrData->cd(2); // corrected
3125     hYieldCorrectedTotal->GetXaxis()->SetMoreLogLabels(kTRUE);
3126     hYieldCorrectedTotal->GetXaxis()->SetNoExponent(kTRUE);
3127     hYieldCorrectedTotal->GetYaxis()->SetRangeUser(minimum, maximum);
3128
3129     if (hMCgenPrimYieldTotal) {
3130       hMCgenPrimYieldTotal->GetYaxis()->SetTitleOffset(1.4);
3131       hMCgenPrimYieldTotal->GetXaxis()->SetMoreLogLabels(kTRUE);
3132       hMCgenPrimYieldTotal->GetXaxis()->SetNoExponent(kTRUE);
3133       hMCgenPrimYieldTotal->GetXaxis()->SetTitle(hYieldCorrectedTotal->GetXaxis()->GetTitle());
3134       hMCgenPrimYieldTotal->GetYaxis()->SetRangeUser(minimum, maximum);
3135     }
3136
3137     if (hMCgenPrimYieldTotal) {
3138       hMCgenPrimYieldTotal->Draw("E1");
3139       hYieldCorrectedTotal->Draw("E1 same");
3140     }
3141     else
3142       hYieldCorrectedTotal->Draw("E1");
3143
3144     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
3145       if (hMCgenPrimYield[i])
3146         hMCgenPrimYield[i]->Draw("E1 same");
3147       hYieldCorrected[i]->Draw("E1 same");
3148     }
3149
3150     TLegend* legTemp = cCorrData->cd(2)->BuildLegend();//0.25, 0.16, 0.65, 0.51);
3151     legTemp->SetNColumns(2);
3152     legTemp->SetFillColor(kWhite);
3153
3154     // Do not include into legend
3155     for (Int_t i = 0; i < AliPID::kSPECIES; i++)
3156       hYieldCorrectedSysError[i]->Draw("E2 same");
3157
3158     ClearTitleFromHistoInCanvas(cCorrData, 2);
3159
3160     TCanvas* cCorrYieldsRatio = 0x0;
3161
3162     TH1D* hYieldCorrectedRatioToMC[AliPID::kSPECIES];
3163     for (Int_t i = 0; i < AliPID::kSPECIES; i++)
3164       hYieldCorrectedRatioToMC[i] = 0x0;
3165
3166     TH1D* hYieldCorrectedTotalRatioToMC = 0x0;
3167
3168     if (numMCgenPrimYieldHistsFound > 0) {
3169       // Compare with MC truth
3170       for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
3171         hYieldCorrectedRatioToMC[species] = new TH1D(*hYieldCorrected[species]);
3172         hYieldCorrectedRatioToMC[species]->SetName(Form("hYieldCorrectedRatioToMC_%s_%s", titles[type].Data(),
3173                                                         AliPID::ParticleShortName(species)));
3174         hYieldCorrectedRatioToMC[species]->SetTitle(Form("%s", AliPID::ParticleLatexName(species)));
3175         hYieldCorrectedRatioToMC[species]->GetYaxis()->SetTitle("Corrected Yield: Fit / MC Truth");
3176         hYieldCorrectedRatioToMC[species]->Divide(hMCgenPrimYield[species]);
3177       }
3178
3179       hYieldCorrectedTotalRatioToMC = new TH1D(*hYieldCorrectedTotal);
3180       hYieldCorrectedTotalRatioToMC->SetName(Form("hYieldCorrectedTotalRatioToMC_%s", titles[type].Data()));
3181       hYieldCorrectedTotalRatioToMC->SetTitle("Total yield");
3182       hYieldCorrectedTotalRatioToMC->GetYaxis()->SetTitle("Corrected Yield: Fit / MC Truth");
3183       hYieldCorrectedTotalRatioToMC->Divide(hMCgenPrimYieldTotal);
3184
3185       cCorrYieldsRatio = new TCanvas(Form("cCorrYieldsRatio_%s", titles[type].Data()), "Corrected Yields Comparison to MC",
3186                                      0, 300, 900, 900);
3187       cCorrYieldsRatio->SetGridx(1);
3188       cCorrYieldsRatio->SetGridy(1);
3189       cCorrYieldsRatio->SetLogx(1);
3190
3191       hYieldCorrectedTotalRatioToMC->GetYaxis()->SetRangeUser(0.6, 1.6);
3192       hYieldCorrectedTotalRatioToMC->GetYaxis()->SetTitleOffset(0.85);
3193       hYieldCorrectedTotalRatioToMC->Draw("E1");
3194
3195       for (Int_t species = 0; species < AliPID::kSPECIES; species++)
3196         hYieldCorrectedRatioToMC[species]->Draw("E1 same");
3197
3198       cCorrYieldsRatio->BuildLegend()->SetFillColor(kWhite);
3199
3200       ClearTitleFromHistoInCanvas(cCorrYieldsRatio);
3201     }
3202
3203     TCanvas* cCorrFractions = new TCanvas(Form("cCorrFractions_%s", titles[type].Data()), "Corrected particleFractions",
3204                                           0, 300, 900, 900);
3205     cCorrFractions->SetLogx(1);
3206     hFractionCorrected[0]->GetYaxis()->SetRangeUser(0., 1.);
3207     hFractionCorrected[0]->Draw("E1");
3208     if (hMCgenPrimFraction[0])
3209       hMCgenPrimFraction[0]->Draw("E1 same");
3210
3211     for (Int_t i = 1; i < AliPID::kSPECIES; i++) {
3212       hFractionCorrected[i]->Draw("E1 same");
3213       if (hMCgenPrimFraction[i])
3214         hMCgenPrimFraction[i]->Draw("E1 same");
3215     }
3216
3217     cCorrFractions->BuildLegend()->SetFillColor(kWhite);
3218
3219     // Do not include into legend
3220     for (Int_t i = 0; i < AliPID::kSPECIES; i++)
3221       hFractionCorrectedSysError[i]->Draw("E2 same");
3222
3223     ClearTitleFromHistoInCanvas(cCorrFractions);
3224
3225     // Save results
3226     saveFile->cd();
3227
3228     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
3229       if (hYield[i])
3230         hYield[i]->Write();
3231     }
3232
3233     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
3234       if (hYieldCorrected[i])
3235         hYieldCorrected[i]->Write();
3236
3237       if (hYieldCorrectedSysError[i])
3238         hYieldCorrectedSysError[i]->Write();
3239     }
3240
3241     if (hYieldCorrectedTotal)
3242       hYieldCorrectedTotal->Write();
3243
3244     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
3245       if (hMCgenPrimYield[i])
3246         hMCgenPrimYield[i]->Write();
3247     }
3248
3249     if (hMCgenPrimYieldTotal)
3250         hMCgenPrimYieldTotal->Write();
3251
3252     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
3253       if (hYieldCorrectedRatioToMC[i])
3254         hYieldCorrectedRatioToMC[i]->Write();
3255     }
3256
3257     if (hYieldCorrectedTotalRatioToMC)
3258         hYieldCorrectedTotalRatioToMC->Write();
3259
3260     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
3261       if (hFractionCorrected[i])
3262         hFractionCorrected[i]->Write();
3263
3264       if (hFractionCorrectedSysError[i])
3265         hFractionCorrectedSysError[i]->Write();
3266     }
3267
3268     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
3269       if (hMCgenPrimFraction[i])
3270         hMCgenPrimFraction[i]->Write();
3271     }
3272
3273     if (cCorrData)
3274       cCorrData->Write();
3275
3276     if (cCorrYieldsRatio)
3277       cCorrYieldsRatio->Write();
3278
3279     if (cCorrFractions)
3280       cCorrFractions->Write();
3281   }
3282
3283   saveFile->cd();
3284
3285   TNamed* settings = new TNamed(
3286       Form("Settings: Efficiency file \"%s\", numJetsForEfficiency file \"%s\", Data file \"%s\", lowerCentrality %.3f, upperCentrality %.3f, lowerJetPt %.1f, upperJetPt %.1f, constantCorrectionAbovePtThreshold %.3f\n",
3287            pathNameEfficiency.Data(), pathNameDataMC.Data(), pathNameData.Data(), lowerCentrality, upperCentrality, lowerJetPt,
3288            upperJetPt, constantCorrectionAbovePtThreshold), "");
3289   settings->Write();
3290
3291   saveFile->Close();
3292
3293
3294   return 0;
3295 }