PWGJE/macros/PID/extractFFs.C

   1 #include "TH1D.h"
   2 #include "TH2D.h"
   3 #include "TCanvas.h"
   4 #include "TFile.h"
   5 #include "TLegend.h"
   6 #include "TROOT.h"
   7 #include "TSystem.h"
   8
   9 #include "AliPID.h"
  10
  11 #include <iostream>
  12
  13 #include "THnSparseDefinitions.h"
  14
  15 #include "../../UserTasks/AliAnalysisTaskPID.h"
  16 #include "SystematicErrorUtils.h"
  17
  18 enum axisDataProj { kPidPtProj = 0, kPidJetPtProj = 1, kPidZProj = 2, kPidXiProj = 3, kNDimsProj = 4 };
  19
  20
  21
  22 //___________________________________________________________________
  23 void setupHist(TH1* h, TString histName, TString histTitle, TString xAxisTitle, TString yAxisTitle, Int_t color)
  24 {
  25   if (histName != "")
  26     h->SetName(histName.Data());
  27   h->SetTitle(histTitle.Data());
  28
  29   if (xAxisTitle != "")
  30     h->GetXaxis()->SetTitle(xAxisTitle.Data());
  31   if (yAxisTitle != "")
  32     h->GetYaxis()->SetTitle(yAxisTitle.Data());
  33
  34   h->SetMarkerStyle(24);
  35   h->SetLineColor(color);
  36   h->SetMarkerColor(color);
  37
  38   h->SetStats(kFALSE);
  39 }
  40
  41
  42 //____________________________________________________________________________________________________________________
  43 void normaliseYieldHist2D(TH2* hData, TH2* hNumJets, const Int_t lowerCentralityBinLimit, const Int_t upperCentralityBinLimit)
  44 {
  45   // Normalise to 1/numJets and bin width. NOTE: jetPt binning of hData and hNumJets assumed to be the same!
  46
  47   if (!hData)
  48     return;
  49
  50   for (Int_t binJetPt = 0; binJetPt <= hData->GetNbinsY() + 1; binJetPt++) {
  51     // No normalisation to numJets, if no histo is provided
  52     const Double_t numJets = hNumJets ? hNumJets->Integral(lowerCentralityBinLimit, upperCentralityBinLimit, binJetPt, binJetPt) : 1.;
  53     Bool_t noJets = numJets < 1e-13;
  54
  55     for (Int_t binObs = 0; binObs <= hData->GetNbinsX() + 1; binObs++) {
  56       if (noJets) {
  57         if (hData->GetBinContent(binObs, binJetPt) > 0.) {
  58           printf("Error: No jets for jetPt ~ %f, but found content %f at y-coord %f!\n", hData->GetYaxis()->GetBinCenter(binJetPt),
  59                  hData->GetBinContent(binObs, binJetPt),  hData->GetXaxis()->GetBinCenter(binObs));
  60         }
  61         continue;
  62       }
  63       const Double_t dObservable = hData->GetXaxis()->GetBinWidth(binObs);
  64       const Double_t normFactor = 1. / (numJets * dObservable);
  65
  66       hData->SetBinContent(binObs, binJetPt, hData->GetBinContent(binObs, binJetPt) * normFactor);
  67       hData->SetBinError(binObs, binJetPt, hData->GetBinError(binObs, binJetPt) * normFactor);
  68     }
  69   }
  70 }
  71
  72 //____________________________________________________________________________________________________________________
  73 void normaliseYieldHist(TH1* h, Double_t numJets)
  74 {
  75   // Normalise to 1/numJets and bin width
  76
  77   if (numJets <= 0) // Do not normalise
  78     numJets = 1.;
  79
  80   for (Int_t bin = 0; bin <= h->GetNbinsX() + 1; bin++) {
  81     const Double_t dObservable = h->GetBinWidth(bin);
  82     const Double_t normFactor = 1. / (numJets * dObservable);
  83     h->SetBinContent(bin, h->GetBinContent(bin) * normFactor);
  84     h->SetBinError(bin, h->GetBinError(bin) * normFactor);
  85   }
  86 }
  87
  88
  89 //___________________________________________________________________
  90 void GenerateParticleYields(THnSparse* hDataProj, AliAnalysisTaskPID* pidTask, const Double_t centrality,
  91                             TH2D** hIDFFtrackPt, TH2D** hIDFFz, TH2D** hIDFFxi,
  92                             const Bool_t setMean, const Bool_t addErrorsQuadratically, const Bool_t smearByError,
  93                             const Bool_t uniformSystematicError, const Bool_t takeIntoAccountSysError, Int_t nGenerations)
  94 {
  95   if (!hDataProj || !pidTask || !hIDFFtrackPt || !hIDFFz || !hIDFFxi) {
  96     printf("Cannot generate particle fractions - missing input!\n");
  97     return;
  98   }
  99
 100   if (takeIntoAccountSysError && smearByError) {
 101     printf("It doesn't make sense to correlate statistical and systematic errors!\n");
 102     return;
 103   }
 104
 105   const Int_t nDimProj = hDataProj->GetNdimensions();
 106   const Long64_t nBinsTHnSparseProj = hDataProj->GetNbins();
 107   Double_t binContent = 0., binError2 = 0.;
 108   Int_t binCoord[nDimProj];
 109   Double_t binCentre[nDimProj];
 110   Double_t prob[AliPID::kSPECIES];
 111
 112   Bool_t success = kTRUE;
 113
 114   // NOTE: In the following, the bin error is divided according to the fraction. The idea is to divide bin content (and error)
 115   // into independent sub-samples according to the fractions. But then, adding up the errors of the sub-samples quadratically
 116   // should recover the original error (think of 100 tracks, error sqrt(100), divided into 10 samples a 10 tracks => error of
 117   // samples should be sqrt(10) then).
 118   // This means, that one needs to add error^2 * fraction and NOT error^2 * fraction^2!!
 119   // However, the errors are ignored anyway at the moment....
 120
 121   if (!smearByError && !takeIntoAccountSysError) {
 122     nGenerations = 1;
 123     // Just calculate fractions/spectra w/o any smearing
 124     const Int_t smearSpeciesByError = -1;
 125     const Int_t takeIntoAccountSysErrorOfSpecies = -1;
 126
 127     for (Long64_t bin = 0; bin < nBinsTHnSparseProj; bin++) {
 128       binContent = hDataProj->GetBinContent(bin, binCoord);
 129       binError2  = hDataProj->GetBinError2(bin);
 130
 131       // If the bin is empty, do not compute the particle fractions -> This bin contributes nothing anyway
 132       if (binContent < 2. * std::numeric_limits<double>::min())
 133         continue;
 134
 135       for (Int_t dim = 0; dim < nDimProj; dim++)
 136         binCentre[dim] = hDataProj->GetAxis(dim)->GetBinCenter(binCoord[dim]);
 137
 138       // Since there is no smearing, the fraction will stay the same and there is no need to save the result for one iteration in a
 139       // histogram (cf. case smearByError >= 0).
 140       success = pidTask->GetParticleFractions(binCentre[kPidPtProj], binCentre[kPidJetPtProj], centrality, prob,
 141                                               smearSpeciesByError, takeIntoAccountSysErrorOfSpecies, uniformSystematicError);
 142
 143       if (success) {
 144         for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
 145           // Track pT
 146           hIDFFtrackPt[species]->SetBinContent(binCoord[kPidPtProj], binCoord[kPidJetPtProj],
 147                                               hIDFFtrackPt[species]->GetBinContent(binCoord[kPidPtProj], binCoord[kPidJetPtProj])
 148                                               + binContent * prob[species]);
 149           Double_t tempError = hIDFFtrackPt[species]->GetBinError(binCoord[kPidPtProj], binCoord[kPidJetPtProj]);
 150           hIDFFtrackPt[species]->SetBinError(binCoord[kPidPtProj], binCoord[kPidJetPtProj],
 151                                             TMath::Sqrt(tempError * tempError + binError2 * prob[species]));
 152
 153           // z
 154           hIDFFz[species]->SetBinContent(binCoord[kPidZProj], binCoord[kPidJetPtProj],
 155                                         hIDFFz[species]->GetBinContent(binCoord[kPidZProj], binCoord[kPidJetPtProj])
 156                                         + binContent * prob[species]);
 157           tempError = hIDFFz[species]->GetBinError(binCoord[kPidZProj], binCoord[kPidJetPtProj]);
 158           hIDFFz[species]->SetBinError(binCoord[kPidZProj], binCoord[kPidJetPtProj],
 159                                       TMath::Sqrt(tempError * tempError + binError2 * prob[species]));
 160
 161           // xi
 162           hIDFFxi[species]->SetBinContent(binCoord[kPidXiProj], binCoord[kPidJetPtProj],
 163                                           hIDFFxi[species]->GetBinContent(binCoord[kPidXiProj], binCoord[kPidJetPtProj])
 164                                           + binContent * prob[species]);
 165           tempError = hIDFFxi[species]->GetBinError(binCoord[kPidXiProj], binCoord[kPidJetPtProj]);
 166           hIDFFxi[species]->SetBinError(binCoord[kPidXiProj], binCoord[kPidJetPtProj],
 167                                         TMath::Sqrt(tempError * tempError + binError2 * prob[species]));
 168         }
 169       }
 170       else
 171         printf("Problem obtaining fractions for: trackPt %f, jetPt %f, cent %f\n", binCentre[kPidPtProj], binCentre[kPidJetPtProj],
 172               centrality);
 173     }
 174   }
 175   else {
 176     // Calculate fractions/spectra w/ smearing "nGenerations" times for every species and obtain the statistical or systematic error
 177     // from the spread of the results.
 178     // Note: For a given species, only the spread of the variation of exactly that species is considered, i.e. one is not
 179     // looking at the change of the fraction if another species is varied (this would bias the result towards smaller errors
 180     // by construction since the fraction of the considered species changes weighted with its statistical/systematic error)
 181
 182     TH2D* hIDFFtrackPtVaried[AliPID::kSPECIES][nGenerations];
 183     TH2D* hIDFFzVaried[AliPID::kSPECIES][nGenerations];
 184     TH2D* hIDFFxiVaried[AliPID::kSPECIES][nGenerations];
 185
 186     // Create this histo and set all entries to -1 (= not calculated yet) in every loop
 187     TH2D* hPartFractionsSmeared = hDataProj->Projection(kPidJetPtProj, kPidPtProj);
 188     hPartFractionsSmeared->SetName("hPartFractionsSmeared");
 189
 190     // Generate results with varying fractions
 191     for (Int_t smearSpeciesByError = 0; smearSpeciesByError < AliPID::kSPECIES; smearSpeciesByError++) {
 192       for (Int_t i = 0; i < nGenerations; i++) {
 193         hIDFFtrackPtVaried[smearSpeciesByError][i] = new TH2D(*(TH2D*)hIDFFtrackPt[smearSpeciesByError]);
 194         hIDFFtrackPtVaried[smearSpeciesByError][i]->SetName(Form("hIDFFtrackPtVaried_%d_%d", smearSpeciesByError, i));
 195         hIDFFtrackPtVaried[smearSpeciesByError][i]->Reset();
 196
 197         hIDFFzVaried[smearSpeciesByError][i] = new TH2D(*(TH2D*)hIDFFz[smearSpeciesByError]);
 198         hIDFFzVaried[smearSpeciesByError][i]->SetName(Form("hIDFFzVaried_%d_%d", smearSpeciesByError, i));
 199         hIDFFzVaried[smearSpeciesByError][i]->Reset();
 200
 201         hIDFFxiVaried[smearSpeciesByError][i] = new TH2D(*(TH2D*)hIDFFxi[smearSpeciesByError]);
 202         hIDFFxiVaried[smearSpeciesByError][i]->SetName(Form("hIDFFxiVaried_%d_%d", smearSpeciesByError, i));
 203         hIDFFxiVaried[smearSpeciesByError][i]->Reset();
 204
 205         // In each iteration (moving over all bins), one want only ONE fixed particle fraction per bin in the spectra map.
 206         // Otherwise, one would assing different fractions to e.g. the same trackPt bin, if only z and xi is different
 207         // (but jetPt bin is still the same). This would then cause the fluctuations to cancel partially.
 208         // Thus: Calculate the smeared fraction for each bin in the spectra map only once and store it in a histo.
 209         // If it is requested again (i.e. histoEntry >= 0), use the value from the histo.
 210
 211         // Set all entries to -1 (= not calculated yet)
 212         for (Int_t iX = 0; iX <= hPartFractionsSmeared->GetNbinsX(); iX++) {
 213           for (Int_t iY = 0; iY <= hPartFractionsSmeared->GetNbinsY(); iY++) {
 214             hPartFractionsSmeared->SetBinContent(iX, iY, -1);
 215           }
 216         }
 217
 218         for (Long64_t bin = 0; bin < nBinsTHnSparseProj; bin++) {
 219           binContent = hDataProj->GetBinContent(bin, binCoord);
 220           binError2  = hDataProj->GetBinError2(bin);
 221
 222           // If the bin is empty, do not compute the particle fractions -> This bin contributes nothing anyway
 223           if (binContent < 2. * std::numeric_limits<double>::min())
 224             continue;
 225
 226           for (Int_t iS = 0; iS < AliPID::kSPECIES; iS++)
 227             prob[iS] = 0.;
 228
 229           for (Int_t dim = 0; dim < nDimProj; dim++)
 230             binCentre[dim] = hDataProj->GetAxis(dim)->GetBinCenter(binCoord[dim]);
 231
 232           if (hPartFractionsSmeared->GetBinContent(binCoord[kPidPtProj], binCoord[kPidJetPtProj]) >= 0) {
 233             success = kTRUE;
 234             prob[smearSpeciesByError] = hPartFractionsSmeared->GetBinContent(binCoord[kPidPtProj], binCoord[kPidJetPtProj]);
 235           }
 236           else {
 237             const Int_t smearSpeciesByStatisticalError = smearByError ? smearSpeciesByError : -1;
 238             const Int_t takeIntoAccountSysErrorOfSpecies = takeIntoAccountSysError ? smearSpeciesByError : -1;
 239
 240             success = pidTask->GetParticleFractions(binCentre[kPidPtProj], binCentre[kPidJetPtProj], centrality, prob,
 241                                                     smearSpeciesByStatisticalError, takeIntoAccountSysErrorOfSpecies,
 242                                                     uniformSystematicError);
 243
 244             if (success)
 245               hPartFractionsSmeared->SetBinContent(binCoord[kPidPtProj], binCoord[kPidJetPtProj], prob[smearSpeciesByError]);
 246           }
 247
 248           // NOTE: Since only the bin content of the current species is stored in hPartFractionsSmeared,
 249           // only prob[smearSpeciesByError] can be used in the following!!!!
 250
 251           if (success) {
 252             // To make things readable...
 253             TH2D* hIDFFtrackPtVariedCurr = hIDFFtrackPtVaried[smearSpeciesByError][i];
 254             TH2D* hIDFFzVariedCurr       = hIDFFzVaried[smearSpeciesByError][i];
 255             TH2D* hIDFFxiVariedCurr      = hIDFFxiVaried[smearSpeciesByError][i];
 256
 257             // Track pT
 258             hIDFFtrackPtVariedCurr->SetBinContent(binCoord[kPidPtProj], binCoord[kPidJetPtProj],
 259                                                   hIDFFtrackPtVariedCurr->GetBinContent(binCoord[kPidPtProj], binCoord[kPidJetPtProj])
 260                                                   + binContent * prob[smearSpeciesByError]);
 261             Double_t tempError = hIDFFtrackPtVariedCurr->GetBinError(binCoord[kPidPtProj], binCoord[kPidJetPtProj]);
 262             hIDFFtrackPtVariedCurr->SetBinError(binCoord[kPidPtProj], binCoord[kPidJetPtProj],
 263                                                 TMath::Sqrt(tempError * tempError + binError2 * prob[smearSpeciesByError]));
 264
 265             // z
 266             hIDFFzVariedCurr->SetBinContent(binCoord[kPidZProj], binCoord[kPidJetPtProj],
 267                                             hIDFFzVariedCurr->GetBinContent(binCoord[kPidZProj], binCoord[kPidJetPtProj])
 268                                             + binContent * prob[smearSpeciesByError]);
 269             tempError = hIDFFzVariedCurr->GetBinError(binCoord[kPidZProj], binCoord[kPidJetPtProj]);
 270             hIDFFzVariedCurr->SetBinError(binCoord[kPidZProj], binCoord[kPidJetPtProj],
 271                                           TMath::Sqrt(tempError * tempError + binError2 * prob[smearSpeciesByError]));
 272
 273             // xi
 274             hIDFFxiVariedCurr->SetBinContent(binCoord[kPidXiProj], binCoord[kPidJetPtProj],
 275                                               hIDFFxiVariedCurr->GetBinContent(binCoord[kPidXiProj], binCoord[kPidJetPtProj])
 276                                               + binContent * prob[smearSpeciesByError]);
 277             tempError = hIDFFxiVariedCurr->GetBinError(binCoord[kPidXiProj], binCoord[kPidJetPtProj]);
 278             hIDFFxiVariedCurr->SetBinError(binCoord[kPidXiProj], binCoord[kPidJetPtProj],
 279                                             TMath::Sqrt(tempError * tempError + binError2 * prob[smearSpeciesByError]));
 280           }
 281           else
 282             printf("Problem obtaining fractions for: trackPt %f, jetPt %f, cent %f, smearSpeciesByError %d\n",
 283                    binCentre[kPidPtProj], binCentre[kPidJetPtProj], centrality, smearSpeciesByError);
 284         }
 285       }
 286     }
 287
 288     delete hPartFractionsSmeared;
 289
 290     const Int_t nHistos = nGenerations;
 291
 292
 293     /*OLD error for each species by variation of ALL species (will bias towards smaller errors!)
 294     // TODO Still does not store all the values in histogram to avoid different fractions for same map bin in same iteration.
 295     // => If this is build in, one needs a histogram for EACH species!!!
 296
 297     // Calculate fractions/spectra w/ smearing "nGenerations" times for every species and obtain the statistical error from
 298     // the spread of the results
 299
 300     TH2D* hIDFFtrackPtVaried[AliPID::kSPECIES][AliPID::kSPECIES * nGenerations];
 301     TH2D* hIDFFzVaried[AliPID::kSPECIES][AliPID::kSPECIES * nGenerations];
 302     TH2D* hIDFFxiVaried[AliPID::kSPECIES][AliPID::kSPECIES * nGenerations];
 303
 304     // Generate results with varying fractions
 305     for (Int_t smearSpeciesByError = 0; smearSpeciesByError < AliPID::kSPECIES; smearSpeciesByError++) {
 306       for (Int_t i = 0; i < nGenerations; i++) {
 307         for (Int_t consideredSpecies = 0; consideredSpecies < AliPID::kSPECIES; consideredSpecies++) {
 308           hIDFFtrackPtVaried[consideredSpecies][smearSpeciesByError * nGenerations + i] =
 309             new TH2D(*(TH2D*)hIDFFtrackPt[consideredSpecies]);
 310           hIDFFtrackPtVaried[consideredSpecies][smearSpeciesByError * nGenerations + i]->Reset();
 311           hIDFFtrackPtVaried[consideredSpecies][smearSpeciesByError * nGenerations + i]->SetName(Form("hIDFFtrackPtVaried_%d_%d",
 312                                                                                           consideredSpecies,
 313                                                                                           smearSpeciesByError * nGenerations + i));
 314           hIDFFzVaried[consideredSpecies][smearSpeciesByError * nGenerations + i] =
 315             new TH2D(*(TH2D*)hIDFFz[consideredSpecies]);
 316           hIDFFzVaried[consideredSpecies][smearSpeciesByError * nGenerations + i]->Reset();
 317           hIDFFzVaried[consideredSpecies][smearSpeciesByError * nGenerations + i]->SetName(Form("hIDFFzVaried_%d_%d",
 318                                                                                           consideredSpecies,
 319                                                                                           smearSpeciesByError * nGenerations + i));
 320           hIDFFxiVaried[consideredSpecies][smearSpeciesByError * nGenerations + i] =
 321             new TH2D(*(TH2D*)hIDFFxi[consideredSpecies]);
 322           hIDFFxiVaried[consideredSpecies][smearSpeciesByError * nGenerations + i]->Reset();
 323           hIDFFxiVaried[consideredSpecies][smearSpeciesByError * nGenerations + i]->SetName(Form("hIDFFxiVaried_%d_%d",
 324                                                                                           consideredSpecies,
 325                                                                                           smearSpeciesByError * nGenerations + i));
 326         }
 327
 328         for (Long64_t bin = 0; bin < nBinsTHnSparseProj; bin++) {
 329           binContent = hDataProj->GetBinContent(bin, binCoord);
 330           binError2  = hDataProj->GetBinError2(bin);
 331
 332           // If the bin is empty, do not compute the particle fractions -> This bin contributes nothing anyway
 333           if (binContent < 2. * std::numeric_limits<double>::min())
 334             continue;
 335
 336           for (Int_t dim = 0; dim < nDimProj; dim++)
 337             binCentre[dim] = hDataProj->GetAxis(dim)->GetBinCenter(binCoord[dim]);
 338
 339           success = pidTask->GetParticleFractions(binCentre[kPidPtProj], binCentre[kPidJetPtProj], centrality, prob,
 340                                                   smearSpeciesByError);
 341
 342           if (success) {
 343             for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
 344               // To make things readable...
 345               TH2D* hIDFFtrackPtVariedCurr = hIDFFtrackPtVaried[species][smearSpeciesByError * nGenerations + i];
 346               TH2D* hIDFFzVariedCurr       = hIDFFzVaried[species][smearSpeciesByError * nGenerations + i];
 347               TH2D* hIDFFxiVariedCurr      = hIDFFxiVaried[species][smearSpeciesByError * nGenerations + i];
 348
 349               // Track pT
 350               hIDFFtrackPtVariedCurr->SetBinContent(binCoord[kPidPtProj], binCoord[kPidJetPtProj],
 351                                                     hIDFFtrackPtVariedCurr->GetBinContent(binCoord[kPidPtProj], binCoord[kPidJetPtProj])
 352                                                     + binContent * prob[species]);
 353               Double_t tempError = hIDFFtrackPtVariedCurr->GetBinError(binCoord[kPidPtProj], binCoord[kPidJetPtProj]);
 354               hIDFFtrackPtVariedCurr->SetBinError(binCoord[kPidPtProj], binCoord[kPidJetPtProj],
 355                                                   TMath::Sqrt(tempError * tempError + binError2 * prob[species]));
 356
 357               // z
 358               hIDFFzVariedCurr->SetBinContent(binCoord[kPidZProj], binCoord[kPidJetPtProj],
 359                                               hIDFFzVariedCurr->GetBinContent(binCoord[kPidZProj], binCoord[kPidJetPtProj])
 360                                               + binContent * prob[species]);
 361               tempError = hIDFFzVariedCurr->GetBinError(binCoord[kPidZProj], binCoord[kPidJetPtProj]);
 362               hIDFFzVariedCurr->SetBinError(binCoord[kPidZProj], binCoord[kPidJetPtProj],
 363                                             TMath::Sqrt(tempError * tempError + binError2 * prob[species]));
 364
 365               // xi
 366               hIDFFxiVariedCurr->SetBinContent(binCoord[kPidXiProj], binCoord[kPidJetPtProj],
 367                                                hIDFFxiVariedCurr->GetBinContent(binCoord[kPidXiProj], binCoord[kPidJetPtProj])
 368                                                + binContent * prob[species]);
 369               tempError = hIDFFxiVariedCurr->GetBinError(binCoord[kPidXiProj], binCoord[kPidJetPtProj]);
 370               hIDFFxiVariedCurr->SetBinError(binCoord[kPidXiProj], binCoord[kPidJetPtProj],
 371                                              TMath::Sqrt(tempError * tempError + binError2 * prob[species]));
 372             }
 373           }
 374           else
 375             printf("Problem obtaining fractions for: trackPt %f, jetPt %f, cent %f, smearSpeciesByError %d\n",
 376                    binCentre[kPidPtProj], binCentre[kPidJetPtProj], centrality, smearSpeciesByError);
 377         }
 378       }
 379     }
 380
 381     const Int_t nHistos = AliPID::kSPECIES * nGenerations;
 382     */
 383
 384     // Compare results to obtain error
 385     const Bool_t ignoreSigmaErrors = kTRUE;
 386
 387     for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
 388       if (!extractSystematicError(nHistos, hIDFFtrackPtVaried[species], hIDFFtrackPt[species], setMean, addErrorsQuadratically,
 389                                   ignoreSigmaErrors))
 390         printf("Failed to determine systematic error for trackPt, species %d\n", species);
 391
 392       if (!extractSystematicError(nHistos, hIDFFzVaried[species], hIDFFz[species], setMean, addErrorsQuadratically, ignoreSigmaErrors))
 393         printf("Failed to determine systematic error for z, species %d\n", species);
 394
 395       if (!extractSystematicError(nHistos, hIDFFxiVaried[species], hIDFFxi[species], setMean, addErrorsQuadratically,
 396                                   ignoreSigmaErrors))
 397         printf("Failed to determine systematic error for xi, species %d\n", species);
 398
 399       for (Int_t i = 0; i < nHistos; i++) {
 400         delete hIDFFtrackPtVaried[species][i];
 401         hIDFFtrackPtVaried[species][i] = 0x0;
 402
 403         delete hIDFFzVaried[species][i];
 404         hIDFFzVaried[species][i] = 0x0;
 405
 406         delete hIDFFxiVaried[species][i];
 407         hIDFFxiVaried[species][i] = 0x0;
 408       }
 409     }
 410   }
 411 }
 412
 413
 414 //___________________________________________________________________
 415 Int_t extractFFs(TString particleFractionPackagePathName, TString pathNameData, TString listName /*= ""*/,
 416                  Bool_t uniformSystematicError, Int_t chargeMode /*kNegCharge = -1, kAllCharged = 0, kPosCharge = 1*/,
 417                  Double_t lowerCentrality = -2, Double_t upperCentrality = -2, Int_t rebinPt = 1, Int_t rebinZ = 1, Int_t rebinXi = 1)
 418 {
 419   TObjArray* histList = 0x0;
 420
 421   if (listName == "") {
 422     listName = pathNameData;
 423     listName.Replace(0, listName.Last('/') + 1, "");
 424     listName.ReplaceAll(".root", "");
 425   }
 426
 427   TFile* f = TFile::Open(pathNameData.Data());
 428   if (!f)  {
 429     std::cout << std::endl;
 430     std::cout << "Failed to open file \"" << pathNameData.Data() << "\"!" << std::endl;
 431     return -1;
 432   }
 433
 434   histList = (TObjArray*)(f->Get(listName.Data()));
 435   if (!histList) {
 436     std::cout << std::endl;
 437     std::cout << "Failed to load list \"" << listName.Data() << "\"!" << std::endl;
 438     return -1;
 439   }
 440
 441   // Extract the data histogram
 442   THnSparse* hPIDdata = dynamic_cast<THnSparse*>(histList->FindObject("hPIDdataAll"));
 443   if (!hPIDdata) {
 444     std::cout << std::endl;
 445     std::cout << "Failed to load data histo!" << std::endl;
 446     return -1;
 447   }
 448
 449   // Set proper errors, if not yet calculated
 450   if (!hPIDdata->GetCalculateErrors()) {
 451     std::cout << "Re-calculating errors of " << hPIDdata->GetName() << "..." << std::endl;
 452     hPIDdata->Sumw2();
 453     Long64_t nBinsTHnSparse = hPIDdata->GetNbins();
 454     Double_t binContent = 0;
 455
 456     for (Long64_t bin = 0; bin < nBinsTHnSparse; bin++) {
 457       binContent = hPIDdata->GetBinContent(bin);
 458       hPIDdata->SetBinError(bin, TMath::Sqrt(binContent));
 459     }
 460   }
 461
 462   // If desired, restrict centrality axis
 463   Int_t lowerCentralityBinLimit = -1;
 464   Int_t upperCentralityBinLimit = -2; // Integral(lowerCentBinLimit, uppCentBinLimit) will not be restricted if these values are kept
 465   Bool_t restrictCentralityAxis = kFALSE;
 466   Double_t actualLowerCentrality = -1.;
 467   Double_t actualUpperCentrality = -1.;
 468
 469   if (lowerCentrality >= -1 && upperCentrality >= -1) {
 470     // Add subtract a very small number to avoid problems with values right on the border between to bins
 471     lowerCentralityBinLimit = hPIDdata->GetAxis(kPidCentrality)->FindBin(lowerCentrality + 0.001);
 472     upperCentralityBinLimit = hPIDdata->GetAxis(kPidCentrality)->FindBin(upperCentrality - 0.001);
 473
 474     // Check if the values look reasonable
 475     if (lowerCentralityBinLimit <= upperCentralityBinLimit && lowerCentralityBinLimit >= 1
 476         && upperCentralityBinLimit <= hPIDdata->GetAxis(kPidCentrality)->GetNbins()) {
 477       actualLowerCentrality = hPIDdata->GetAxis(kPidCentrality)->GetBinLowEdge(lowerCentralityBinLimit);
 478       actualUpperCentrality = hPIDdata->GetAxis(kPidCentrality)->GetBinUpEdge(upperCentralityBinLimit);
 479
 480       restrictCentralityAxis = kTRUE;
 481     }
 482     else {
 483       std::cout << std::endl;
 484       std::cout << "Requested centrality range out of limits or upper and lower limit are switched!" << std::endl;
 485       return -1;
 486     }
 487   }
 488
 489   std::cout << "centrality: ";
 490   if (restrictCentralityAxis) {
 491     std::cout << actualLowerCentrality << " - " << actualUpperCentrality << std::endl;
 492   }
 493   else {
 494     std::cout << "All" << std::endl;
 495   }
 496
 497   if (restrictCentralityAxis) {
 498     hPIDdata->GetAxis(kPidCentrality)->SetRange(lowerCentralityBinLimit, upperCentralityBinLimit);
 499   }
 500
 501   // If desired, restrict charge axis
 502   std::cout << "Charge selection: ";
 503   if (chargeMode == kAllCharged)
 504     std::cout << "All charged particles" << std::endl;
 505   else if (chargeMode == kNegCharge)
 506     std::cout << "Negative particles only" << std::endl;
 507   else if (chargeMode == kPosCharge)
 508     std::cout << "Positive particles only" << std::endl;
 509   else {
 510     std::cout << "Unknown -> ERROR" << std::endl;
 511     return -1;
 512   }
 513
 514   const Bool_t restrictCharge = (chargeMode != kAllCharged);
 515
 516   const Int_t indexChargeAxisData = GetAxisByTitle(hPIDdata, "Charge (e_{0})");
 517   if (indexChargeAxisData < 0 && restrictCharge) {
 518     std::cout << "Error: Charge axis not found for data histogram!" << std::endl;
 519     return -1;
 520   }
 521   Int_t lowerChargeBinLimitData = -1;
 522   Int_t upperChargeBinLimitData = -2;
 523   Double_t actualLowerChargeData = -999;
 524   Double_t actualUpperChargeData = -999;
 525
 526   if (restrictCharge) {
 527     // Add subtract a very small number to avoid problems with values right on the border between to bins
 528     if (chargeMode == kNegCharge) {
 529       lowerChargeBinLimitData = hPIDdata->GetAxis(indexChargeAxisData)->FindBin(-1. + 0.001);
 530       upperChargeBinLimitData = hPIDdata->GetAxis(indexChargeAxisData)->FindBin(0. - 0.001);
 531     }
 532     else if (chargeMode == kPosCharge) {
 533       lowerChargeBinLimitData = hPIDdata->GetAxis(indexChargeAxisData)->FindBin(0. + 0.001);
 534       upperChargeBinLimitData = hPIDdata->GetAxis(indexChargeAxisData)->FindBin(1. - 0.001);
 535     }
 536
 537     // Check if the values look reasonable
 538     if (lowerChargeBinLimitData <= upperChargeBinLimitData && lowerChargeBinLimitData >= 1
 539         && upperChargeBinLimitData <= hPIDdata->GetAxis(indexChargeAxisData)->GetNbins()) {
 540       actualLowerChargeData = hPIDdata->GetAxis(indexChargeAxisData)->GetBinLowEdge(lowerChargeBinLimitData);
 541       actualUpperChargeData = hPIDdata->GetAxis(indexChargeAxisData)->GetBinUpEdge(upperChargeBinLimitData);
 542
 543       std::cout << "Charge range data: " << actualLowerChargeData << " - " << actualUpperChargeData << std::endl;
 544     }
 545     else {
 546       std::cout << std::endl;
 547       std::cout << "Requested charge range out of limits or upper and lower limit are switched!" << std::endl;
 548       return -1;
 549     }
 550
 551     hPIDdata->GetAxis(indexChargeAxisData)->SetRange(lowerChargeBinLimitData, upperChargeBinLimitData);
 552   }
 553
 554   // Just take one arbitrary selectSpecies to avoid multiple counting
 555   hPIDdata->GetAxis(kPidSelectSpecies)->SetRange(1, 1);
 556
 557   // Get projection on dimensions relevant for FFs
 558   const Int_t nDimProj = kNDimsProj;
 559   Int_t dimProj[nDimProj];
 560   dimProj[kPidPtProj] = kPidPt;
 561   dimProj[kPidJetPtProj] = kPidJetPt;
 562   dimProj[kPidZProj] = kPidZ;
 563   dimProj[kPidXiProj] =kPidXi;
 564
 565   THnSparse* hDataProj = hPIDdata->Projection(nDimProj, dimProj, "e");
 566
 567   // If desired, rebin axes
 568   if (rebinPt != 1 || rebinZ != 1 || rebinXi != 1) {
 569     Int_t group[hDataProj->GetNdimensions()];
 570
 571     for (Int_t i = 0; i < hDataProj->GetNdimensions(); i++) {
 572       if (i == kPidPtProj)
 573         group[i] = rebinPt;
 574       else if (i == kPidZProj)
 575         group[i] = rebinZ;
 576       else if (i == kPidXiProj)
 577         group[i] = rebinXi;
 578       else
 579         group[i] = 1;
 580     }
 581
 582     THnSparse* hTemp = hDataProj->Rebin(group);
 583     hDataProj->SetName("temp");
 584     delete hDataProj;
 585
 586     hDataProj = hTemp;
 587   }
 588
 589   /* OLD - Now normalisation to nJets
 590   Double_t numEvents = -1;
 591   TH1* hNumEvents = dynamic_cast<TH1*>(histList->FindObject("fhEventsProcessed"));
 592   if (!hNumEvents) {
 593     std::cout << std::endl;
 594     std::cout << "Histo with number of processed events not found! Yields will NOT be normalised to this number!" << std::endl
 595               << std::endl;
 596   }
 597   else {
 598     numEvents = hNumEvents->Integral(lowerCentralityBinLimit, upperCentralityBinLimit);
 599
 600     if (numEvents <= 0) {
 601       numEvents = -1;
 602       std::cout << std::endl;
 603       std::cout << "Number of processed events < 1 in selected range! Yields will NOT be normalised to this number!"
 604                 << std::endl << std::endl;
 605     }
 606   }*/
 607
 608
 609   TH2D* hNjetsGen = 0x0;
 610   TH2D* hNjetsRec = 0x0;
 611
 612   TH2D* hMCgenPrimYieldPt[AliPID::kSPECIES];
 613   TH2D* hMCgenPrimYieldZ[AliPID::kSPECIES];
 614   TH2D* hMCgenPrimYieldXi[AliPID::kSPECIES];
 615
 616   hNjetsGen = (TH2D*)histList->FindObject("fh2FFJetPtGen");
 617   hNjetsRec = (TH2D*)histList->FindObject("fh2FFJetPtRec");
 618
 619   if (!hNjetsRec) {
 620     printf("Failed to load number of jets (rec) histo!\n");
 621
 622     // For backward compatibility (TODO REMOVE IN FUTURE): Load info from fixed AnalysisResults file (might be wrong, if other
 623     // period is considered; also: No multiplicity information)
 624     TFile* fBackward = TFile::Open("finalCuts/pp/7TeV/LHC10e.pass2/corrected/finalisedSplines/finalMapsAndTail/Jets/noCutOn_ncl_or_liav/AnalysisResults.root");
 625
 626     TString dirDataInFile = "";
 627     TDirectory* dirData = fBackward ? (TDirectory*)fBackward->Get(fBackward->GetListOfKeys()->At(0)->GetName()) : 0x0;
 628
 629     TList* list = dirData ? (TList*)dirData->Get(dirData->GetListOfKeys()->At(0)->GetName()) : 0x0;
 630
 631     TH1D* hFFJetPtRec = list ? (TH1D*)list->FindObject("fh1FFJetPtRecCutsInc") : 0x0;
 632
 633     if (hFFJetPtRec) {
 634       printf("***WARNING: For backward compatibility, using file \"finalCuts/pp/7TeV/LHC10e.pass2/corrected/finalisedSplines/finalMapsAndTail/Jets/noCutOn_ncl_or_liav/AnalysisResults.root\" to get number of jets. BUT: Might be wrong period and has no mult info!***\n");
 635       printf("ALSO: Using Njets for inclusive jets!!!!\n");
 636
 637       hNjetsRec = new TH2D("fh2FFJetPtRec", "", 1, -1, 1,  hDataProj->GetAxis(kPidJetPtProj)->GetNbins(),
 638                           hDataProj->GetAxis(kPidJetPtProj)->GetXbins()->GetArray());
 639
 640       for (Int_t iJet = 1; iJet <= hNjetsRec->GetNbinsY(); iJet++) {
 641         Int_t lowerBin = hFFJetPtRec->FindBin(hNjetsRec->GetYaxis()->GetBinLowEdge(iJet) + 1e-3);
 642         Int_t upperBin = hFFJetPtRec->FindBin(hNjetsRec->GetYaxis()->GetBinUpEdge(iJet) - 1e-3);
 643         hNjetsRec->SetBinContent(1, iJet, hFFJetPtRec->Integral(lowerBin, upperBin));
 644       }
 645     }
 646
 647     if (!hNjetsRec)
 648       return -1;
 649   }
 650
 651   THnSparse* hMCgeneratedYieldsPrimaries = (THnSparse*)histList->FindObject("fhMCgeneratedYieldsPrimaries");
 652
 653   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
 654     hMCgenPrimYieldPt[i] = 0x0;
 655     hMCgenPrimYieldZ[i] = 0x0;
 656     hMCgenPrimYieldXi[i] = 0x0;
 657   }
 658
 659   if (hMCgeneratedYieldsPrimaries && hMCgeneratedYieldsPrimaries->GetEntries() > 0) {
 660     // Set proper errors, if not yet calculated
 661     if (!hMCgeneratedYieldsPrimaries->GetCalculateErrors()) {
 662       std::cout << "Re-calculating errors of " << hMCgeneratedYieldsPrimaries->GetName() << "..." << std::endl;
 663
 664       hMCgeneratedYieldsPrimaries->Sumw2();
 665
 666       Long64_t nBinsTHnSparseGenYield = hMCgeneratedYieldsPrimaries->GetNbins();
 667       Double_t binContentGenYield = 0;
 668       for (Long64_t bin = 0; bin < nBinsTHnSparseGenYield; bin++) {
 669         binContentGenYield = hMCgeneratedYieldsPrimaries->GetBinContent(bin);
 670         hMCgeneratedYieldsPrimaries->SetBinError(bin, TMath::Sqrt(binContentGenYield));
 671       }
 672     }
 673
 674     // If desired, rebin axes
 675     if (rebinPt != 1 || rebinZ != 1 || rebinXi != 1) {
 676       Int_t group[hMCgeneratedYieldsPrimaries->GetNdimensions()];
 677
 678       for (Int_t k = 0; k < hMCgeneratedYieldsPrimaries->GetNdimensions(); k++) {
 679         if (k == kPidGenYieldPt)
 680           group[k] = rebinPt;
 681         else if (k == kPidGenYieldZ)
 682           group[k] = rebinZ;
 683         else if (k == kPidGenYieldXi)
 684           group[k] = rebinXi;
 685         else
 686           group[k] = 1;
 687       }
 688
 689       THnSparse* hTemp = hMCgeneratedYieldsPrimaries->Rebin(group);
 690       hMCgeneratedYieldsPrimaries->SetName("temp");
 691       delete hMCgeneratedYieldsPrimaries;
 692
 693       hMCgeneratedYieldsPrimaries = hTemp;
 694     }
 695
 696
 697     if (restrictCentralityAxis)
 698       hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldCentrality)->SetRange(lowerCentralityBinLimit, upperCentralityBinLimit);
 699
 700     if (restrictCharge) {
 701       const Int_t indexChargeAxisGenYield = GetAxisByTitle(hMCgeneratedYieldsPrimaries, "Charge (e_{0})");
 702       if (indexChargeAxisGenYield < 0) {
 703         std::cout << "Error: Charge axis not found for gen yield histogram!" << std::endl;
 704         return -1;
 705       }
 706
 707       Int_t lowerChargeBinLimitGenYield = -1;
 708       Int_t upperChargeBinLimitGenYield = -2;
 709       Double_t actualLowerChargeGenYield = -999;
 710       Double_t actualUpperChargeGenYield = -999;
 711
 712       // Add subtract a very small number to avoid problems with values right on the border between to bins
 713       if (chargeMode == kNegCharge) {
 714         lowerChargeBinLimitGenYield = hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->FindBin(-1. + 0.001);
 715         upperChargeBinLimitGenYield = hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->FindBin(0. - 0.001);
 716       }
 717       else if (chargeMode == kPosCharge) {
 718         lowerChargeBinLimitGenYield = hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->FindBin(0. + 0.001);
 719         upperChargeBinLimitGenYield = hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->FindBin(1. - 0.001);
 720       }
 721
 722       // Check if the values look reasonable
 723       if (lowerChargeBinLimitGenYield <= upperChargeBinLimitGenYield && lowerChargeBinLimitGenYield >= 1
 724           && upperChargeBinLimitGenYield <= hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->GetNbins()) {
 725         actualLowerChargeGenYield = hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->GetBinLowEdge(lowerChargeBinLimitGenYield);
 726         actualUpperChargeGenYield = hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->GetBinUpEdge(upperChargeBinLimitGenYield);
 727
 728         if (TMath::Abs(actualLowerChargeGenYield - actualLowerChargeData) > 1e-4 ||
 729             TMath::Abs(actualUpperChargeGenYield - actualUpperChargeData) > 1e-4) {
 730           std::cout << std::endl;
 731           std::cout << "Error: Charge range gen yield: " << actualLowerChargeGenYield << " - " << actualUpperChargeGenYield
 732                     << std::endl << "differs from that of data: " << actualLowerChargeData << " - " << actualUpperChargeData
 733                     << std::endl;
 734           return -1;
 735         }
 736       }
 737       else {
 738         std::cout << std::endl;
 739         std::cout << "Requested charge range (gen yield) out of limits or upper and lower limit are switched!" << std::endl;
 740         return -1;
 741       }
 742
 743       hMCgeneratedYieldsPrimaries->GetAxis(indexChargeAxisGenYield)->SetRange(lowerChargeBinLimitGenYield,
 744                                                                               upperChargeBinLimitGenYield);
 745     }
 746
 747     for (Int_t MCid = 0; MCid < AliPID::kSPECIES; MCid++) {
 748       hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldMCpid)->SetRange(MCid + 1, MCid + 1);
 749
 750       hMCgenPrimYieldPt[MCid] = hMCgeneratedYieldsPrimaries->Projection(kPidGenYieldJetPt, kPidGenYieldPt, "e");
 751       hMCgenPrimYieldPt[MCid]->SetName(Form("hMCgenYieldsPrimPt_%s", AliPID::ParticleShortName(MCid)));
 752       hMCgenPrimYieldPt[MCid]->SetTitle(Form("MC truth generated primary yield, %s", AliPID::ParticleName(MCid)));
 753       hMCgenPrimYieldPt[MCid]->GetZaxis()->SetTitle("1/N_{Jets} dN/dP_{T} (GeV/c)^{-1}");
 754       hMCgenPrimYieldPt[MCid]->SetStats(kFALSE);
 755
 756       hMCgenPrimYieldZ[MCid] = hMCgeneratedYieldsPrimaries->Projection(kPidGenYieldJetPt, kPidGenYieldZ, "e");
 757       hMCgenPrimYieldZ[MCid]->SetName(Form("hMCgenYieldsPrimZ_%s", AliPID::ParticleShortName(MCid)));
 758       hMCgenPrimYieldZ[MCid]->SetTitle(Form("MC truth generated primary yield, %s", AliPID::ParticleName(MCid)));
 759       hMCgenPrimYieldZ[MCid]->GetZaxis()->SetTitle("1/N_{Jets} dN/dz");
 760       hMCgenPrimYieldZ[MCid]->SetStats(kFALSE);
 761
 762       hMCgenPrimYieldXi[MCid] = hMCgeneratedYieldsPrimaries->Projection(kPidGenYieldJetPt, kPidGenYieldXi, "e");
 763       hMCgenPrimYieldXi[MCid]->SetName(Form("hMCgenYieldsPrimXi_%s", AliPID::ParticleShortName(MCid)));
 764       hMCgenPrimYieldXi[MCid]->SetTitle(Form("MC truth generated primary yield, %s", AliPID::ParticleName(MCid)));
 765       hMCgenPrimYieldXi[MCid]->GetZaxis()->SetTitle("1/N_{Jets} dN/d#xi");
 766       hMCgenPrimYieldXi[MCid]->SetStats(kFALSE);
 767
 768       hMCgeneratedYieldsPrimaries->GetAxis(kPidGenYieldMCpid)->SetRange(0, -1);
 769     }
 770   }
 771
 772
 773   AliAnalysisTaskPID *pidTask = new AliAnalysisTaskPID("spectrumExtractorTask");
 774
 775   if (!pidTask->SetParticleFractionHistosFromFile(particleFractionPackagePathName, kFALSE)) {
 776     printf("Failed to load particle fraction package from file \"%s\"!\n", particleFractionPackagePathName.Data());
 777     return -1;
 778   }
 779
 780   if (!pidTask->SetParticleFractionHistosFromFile(particleFractionPackagePathName, kTRUE)) {
 781     printf("Failed to load particle fraction sys error package from file \"%s\"!\n", particleFractionPackagePathName.Data());
 782     return -1;
 783   }
 784
 785   printf("Loaded particle fraction package from file \"%s\"!\n", particleFractionPackagePathName.Data());
 786
 787   TH2D* hIDFFtrackPt[AliPID::kSPECIES];
 788   TH2D* hIDFFz[AliPID::kSPECIES];
 789   TH2D* hIDFFxi[AliPID::kSPECIES];
 790
 791   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
 792     hIDFFtrackPt[i] = hDataProj->Projection(kPidJetPtProj, kPidPtProj);
 793     hIDFFtrackPt[i]->Reset();
 794     if (hIDFFtrackPt[i]->GetSumw2N() <= 0)
 795       hIDFFtrackPt[i]->Sumw2();
 796     setupHist(hIDFFtrackPt[i], Form("hIDFFtrackPt_%s", AliPID::ParticleShortName(i)), Form("%s", AliPID::ParticleShortName(i)),
 797               hDataProj->GetAxis(kPidPtProj)->GetTitle(), hDataProj->GetAxis(kPidJetPtProj)->GetTitle(), kBlack);
 798
 799     hIDFFz[i] = hDataProj->Projection(kPidJetPtProj, kPidZProj);
 800     hIDFFz[i]->Reset();
 801     if (hIDFFz[i]->GetSumw2N() <= 0)
 802       hIDFFz[i]->Sumw2();
 803     setupHist(hIDFFz[i], Form("hIDFFz_%s", AliPID::ParticleShortName(i)), Form("%s", AliPID::ParticleShortName(i)),
 804               hDataProj->GetAxis(kPidZProj)->GetTitle(), hDataProj->GetAxis(kPidJetPtProj)->GetTitle(), kBlack);
 805
 806     hIDFFxi[i] = hDataProj->Projection(kPidJetPtProj, kPidXiProj);
 807     hIDFFxi[i]->Reset();
 808     if (hIDFFxi[i]->GetSumw2N() <= 0)
 809       hIDFFxi[i]->Sumw2();
 810     setupHist(hIDFFxi[i], Form("hIDFFxi_%s", AliPID::ParticleShortName(i)), Form("%s", AliPID::ParticleShortName(i)),
 811               hDataProj->GetAxis(kPidXiProj)->GetTitle(), hDataProj->GetAxis(kPidJetPtProj)->GetTitle(), kBlack);
 812   }
 813
 814   const Double_t centrality = (actualLowerCentrality + actualUpperCentrality) / 2.;
 815
 816   // First iteration: Just take the default fractions to obtain the "mean" of the yields
 817   Bool_t setMean = kTRUE;
 818   Bool_t addErrorsQuadratically = kFALSE;
 819   Bool_t smearByError = kFALSE;
 820   Bool_t takeIntoAccountSysError = kFALSE;
 821   // setMean and all the follwoing parameters are anyway irrelevant, if smearByError = kFALSE and takeIntoAccountSysError = kFALSE
 822   GenerateParticleYields(hDataProj, pidTask, centrality, hIDFFtrackPt, hIDFFz, hIDFFxi, setMean, addErrorsQuadratically,
 823                          smearByError, uniformSystematicError, takeIntoAccountSysError, 1);
 824
 825
 826   // Next iteration: Vary the PID map within statistical errors several times and calculate the resulting statistical error
 827   // of the spectra from this. But since the mean of the fraction is some kind of "best estimate" of the truth, leave the means
 828   // as they have been set during the last step.
 829
 830
 831   // NOTE: Do NOT add the statistical errors from the last step, since the yields/fractions (rel. error is the same) already
 832   // contain the stat. uncertainty of the yield of this species in the corresponding bin! I.e. one just takes some kind of weighted
 833   // mean of the fractions (mean and error) (equivalent of taking the yields with corresponding error).
 834   setMean = kFALSE;
 835   addErrorsQuadratically = kFALSE;
 836   smearByError = kTRUE;
 837   takeIntoAccountSysError = kFALSE;
 838   const Int_t nGenerations = 5000; //TODO set to 5000 in the end, after all testing was successful
 839   GenerateParticleYields(hDataProj, pidTask, centrality, hIDFFtrackPt, hIDFFz, hIDFFxi, setMean, addErrorsQuadratically,
 840                          smearByError, uniformSystematicError, takeIntoAccountSysError, nGenerations);
 841
 842
 843   // Next iteration: Vary the PID map within systematic errors several times and calculate the resulting systematic error
 844   // of the spectra from this. But since the mean of the fraction is some kind of "best estimate" of the truth, leave the means
 845   // as they have been set during the last step.
 846   setMean = kFALSE;
 847   addErrorsQuadratically = kFALSE;
 848   smearByError = kFALSE;
 849   takeIntoAccountSysError = kTRUE;
 850   // Clone histograms with final statistical errors -> Only set systematic errors, but leave mean as it is
 851   TH2D* hIDFFtrackPtSysError[AliPID::kSPECIES];
 852   TH2D* hIDFFzSysError[AliPID::kSPECIES];
 853   TH2D* hIDFFxiSysError[AliPID::kSPECIES];
 854
 855   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
 856     hIDFFtrackPtSysError[i] = new TH2D(*hIDFFtrackPt[i]);
 857     hIDFFtrackPtSysError[i]->SetName(Form("%s_sysError", hIDFFtrackPt[i]->GetName()));
 858     hIDFFtrackPtSysError[i]->SetFillStyle(0);
 859
 860     hIDFFzSysError[i] = new TH2D(*hIDFFz[i]);
 861     hIDFFzSysError[i]->SetName(Form("%s_sysError", hIDFFz[i]->GetName()));
 862     hIDFFzSysError[i]->SetFillStyle(0);
 863
 864     hIDFFxiSysError[i] = new TH2D(*hIDFFxi[i]);
 865     hIDFFxiSysError[i]->SetName(Form("%s_sysError", hIDFFxi[i]->GetName()));
 866     hIDFFxiSysError[i]->SetFillStyle(0);
 867   }
 868
 869   GenerateParticleYields(hDataProj, pidTask, centrality, hIDFFtrackPtSysError, hIDFFzSysError, hIDFFxiSysError, setMean,
 870                          addErrorsQuadratically, smearByError, uniformSystematicError, takeIntoAccountSysError, nGenerations);
 871
 872   delete pidTask;
 873
 874
 875   // Normalise properly
 876   for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
 877     normaliseYieldHist2D(hIDFFtrackPt[species], hNjetsRec, lowerCentralityBinLimit, upperCentralityBinLimit);
 878     normaliseYieldHist2D(hIDFFz[species], hNjetsRec, lowerCentralityBinLimit, upperCentralityBinLimit);
 879     normaliseYieldHist2D(hIDFFxi[species], hNjetsRec, lowerCentralityBinLimit, upperCentralityBinLimit);
 880
 881     normaliseYieldHist2D(hIDFFtrackPtSysError[species], hNjetsRec, lowerCentralityBinLimit, upperCentralityBinLimit);
 882     normaliseYieldHist2D(hIDFFzSysError[species], hNjetsRec, lowerCentralityBinLimit, upperCentralityBinLimit);
 883     normaliseYieldHist2D(hIDFFxiSysError[species], hNjetsRec, lowerCentralityBinLimit, upperCentralityBinLimit);
 884
 885     normaliseYieldHist2D(hMCgenPrimYieldPt[species], hNjetsGen, lowerCentralityBinLimit, upperCentralityBinLimit);
 886     normaliseYieldHist2D(hMCgenPrimYieldZ[species], hNjetsGen, lowerCentralityBinLimit, upperCentralityBinLimit);
 887     normaliseYieldHist2D(hMCgenPrimYieldXi[species], hNjetsGen, lowerCentralityBinLimit, upperCentralityBinLimit);
 888   }
 889
 890
 891
 892
 893   // Save results to file
 894   TString chargeString = "";
 895   if (chargeMode == kPosCharge)
 896     chargeString = "_posCharge";
 897   else if (chargeMode == kNegCharge)
 898     chargeString = "_negCharge";
 899
 900   TString saveFileName = pathNameData;
 901   saveFileName.Replace(0, pathNameData.Last('/') + 1, "");
 902
 903   TString savePath = pathNameData;
 904   savePath.ReplaceAll(Form("/%s", saveFileName.Data()), "");
 905
 906   saveFileName.Prepend("output_extractedFFs_");
 907   saveFileName.ReplaceAll(".root", Form("_centrality_%s%s.root",
 908                                         restrictCentralityAxis ? Form("%.0f_%.0f.root", actualLowerCentrality, actualUpperCentrality)
 909                                                                : "all",
 910                                         chargeString.Data()));
 911
 912   TString saveFilePathName = Form("%s/%s", savePath.Data(), saveFileName.Data());
 913   TFile* saveFile = TFile::Open(saveFilePathName.Data(), "RECREATE");
 914
 915   if (!saveFile) {
 916     printf("Failed to save results to file \"%s\"!\n", saveFilePathName.Data());
 917     return -1;
 918   }
 919
 920   saveFile->cd();
 921
 922
 923   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
 924     if (hIDFFtrackPt[i])
 925       hIDFFtrackPt[i]->Write();
 926
 927     if (hIDFFz[i])
 928       hIDFFz[i]->Write();
 929
 930     if (hIDFFxi[i])
 931       hIDFFxi[i]->Write();
 932
 933
 934     if (hIDFFtrackPtSysError[i])
 935       hIDFFtrackPtSysError[i]->Write();
 936
 937     if (hIDFFzSysError[i])
 938       hIDFFzSysError[i]->Write();
 939
 940     if (hIDFFxiSysError[i])
 941       hIDFFxiSysError[i]->Write();
 942
 943
 944     if (hMCgenPrimYieldPt[i])
 945       hMCgenPrimYieldPt[i]->Write();
 946
 947     if (hMCgenPrimYieldZ[i])
 948       hMCgenPrimYieldZ[i]->Write();
 949
 950     if (hMCgenPrimYieldXi[i])
 951       hMCgenPrimYieldXi[i]->Write();
 952   }
 953
 954   if (hNjetsGen)
 955     hNjetsGen->Write();
 956
 957   if (hNjetsRec)
 958     hNjetsRec->Write();
 959
 960
 961   TNamed* settings = new TNamed(
 962       Form("Settings: Fraction package \"%s\", Data file \"%s\", lowerCentrality %.0f, upperCentrality %.0f, chargeMode %d, uniformSystematicError %d, rebinPt %d, rebinZ %d, rebinXi %d\n",
 963            particleFractionPackagePathName.Data(), pathNameData.Data(), lowerCentrality, upperCentrality, chargeMode,
 964            uniformSystematicError, rebinPt, rebinZ, rebinXi), "");
 965   settings->Write();
 966
 967   saveFile->Close();
 968
 969
 970
 971   /*
 972   Int_t t1 = hIDFFtrackPt[AliPID::kPion]->GetYaxis()->FindBin(5.1);
 973   Int_t t2 = hIDFFtrackPt[AliPID::kPion]->GetYaxis()->FindBin(9.9);
 974
 975   printf("t1 %d, t2 %d\n", t1, t2);
 976   new TCanvas();
 977   TH1D* hTemp = hIDFFtrackPt[AliPID::kPion]->ProjectionX("_pfy1", t1, t2, "e");
 978   hTemp->SetFillStyle(0);
 979   hTemp->Draw("");
 980
 981   hTemp = hIDFFtrackPtSysError[AliPID::kPion]->ProjectionX("sys_pfy1", t1, t2, "e");
 982   hTemp->SetFillStyle(0);
 983   hTemp->Draw("E2same");
 984
 985   if (hMCgenPrimYieldPt[AliPID::kPion]) {
 986     hTemp = hMCgenPrimYieldPt[AliPID::kPion]->ProjectionX("MC_pfy1", t1, t2, "e");
 987     hTemp->SetFillStyle(0);
 988     hTemp->SetMarkerStyle(23);
 989     hTemp->Draw("same");
 990   }
 991
 992
 993   new TCanvas();
 994   hTemp = hIDFFz[AliPID::kPion]->ProjectionX("_pfy2", t1, t2, "e");
 995   hTemp->SetFillStyle(0);
 996   hTemp->Draw("");
 997
 998   hTemp = hIDFFzSysError[AliPID::kPion]->ProjectionX("sys_pfy2", t1, t2, "e");
 999   hTemp->SetFillStyle(0);
1000   hTemp->Draw("E2same");
1001
1002   if (hMCgenPrimYieldZ[AliPID::kPion]) {
1003     hTemp = hMCgenPrimYieldZ[AliPID::kPion]->ProjectionX("MC_pfy2", t1, t2, "e");
1004     hTemp->SetFillStyle(0);
1005     hTemp->SetMarkerStyle(23);
1006     hTemp->Draw("same");
1007   }
1008
1009
1010   new TCanvas();
1011   hTemp = hIDFFxi[AliPID::kPion]->ProjectionX("_pfy3", t1, t2, "e");
1012   hTemp->SetFillStyle(0);
1013   hTemp->Draw("");
1014
1015   hTemp = hIDFFxiSysError[AliPID::kPion]->ProjectionX("sys_pfy3", t1, t2, "e");
1016   hTemp->SetFillStyle(0);
1017   hTemp->Draw("E2same");
1018
1019   if (hMCgenPrimYieldXi[AliPID::kPion]) {
1020     hTemp = hMCgenPrimYieldXi[AliPID::kPion]->ProjectionX("MC_pfy3", t1, t2, "e");
1021     hTemp->SetFillStyle(0);
1022     hTemp->SetMarkerStyle(23);
1023     hTemp->Draw("same");
1024   }
1025
1026
1027
1028
1029   new TCanvas();
1030   hTemp = hIDFFtrackPt[AliPID::kProton]->ProjectionX("_pfy4", t1, t2, "e");
1031   hTemp->SetFillStyle(0);
1032   hTemp->Draw("");
1033
1034   hTemp = hIDFFtrackPtSysError[AliPID::kProton]->ProjectionX("sys_pfy4", t1, t2, "e");
1035   hTemp->SetFillStyle(0);
1036   hTemp->Draw("E2same");
1037
1038   if (hMCgenPrimYieldPt[AliPID::kProton]) {
1039     hTemp = hMCgenPrimYieldPt[AliPID::kProton]->ProjectionX("MC_pfy4", t1, t2, "e");
1040     hTemp->SetFillStyle(0);
1041     hTemp->SetMarkerStyle(23);
1042     hTemp->Draw("same");
1043   }
1044
1045
1046   new TCanvas();
1047   hTemp = hIDFFz[AliPID::kProton]->ProjectionX("_pfy5", t1, t2, "e");
1048   hTemp->SetFillStyle(0);
1049   hTemp->Draw("");
1050
1051   hTemp = hIDFFzSysError[AliPID::kProton]->ProjectionX("sys_pfy5", t1, t2, "e");
1052   hTemp->SetFillStyle(0);
1053   hTemp->Draw("E2same");
1054
1055   if (hMCgenPrimYieldZ[AliPID::kProton]) {
1056     hTemp = hMCgenPrimYieldZ[AliPID::kProton]->ProjectionX("MC_pfy5", t1, t2, "e");
1057     hTemp->SetFillStyle(0);
1058     hTemp->SetMarkerStyle(23);
1059     hTemp->Draw("same");
1060   }
1061
1062
1063   new TCanvas();
1064   hTemp = hIDFFxi[AliPID::kProton]->ProjectionX("_pfy", t1, t2, "e");
1065   hTemp->SetFillStyle(0);
1066   hTemp->Draw("");
1067
1068   hTemp = hIDFFxiSysError[AliPID::kProton]->ProjectionX("sys_pfy", t1, t2, "e");
1069   hTemp->SetFillStyle(0);
1070   hTemp->Draw("E2same");
1071
1072   if (hMCgenPrimYieldXi[AliPID::kProton]) {
1073     hTemp = hMCgenPrimYieldXi[AliPID::kProton]->ProjectionX("MC_pfy", t1, t2, "e");
1074     hTemp->SetFillStyle(0);
1075     hTemp->SetMarkerStyle(23);
1076     hTemp->Draw("same");
1077   }*/
1078
1079   return 0;
1080 }