PWGJE/UserTasks/AliAnalysisTaskPID.cxx

   1 #include "TChain.h"
   2 #include "TFile.h"
   3 #include "TF1.h"
   4 #include "TAxis.h"
   5 #include "TProfile.h"
   6 #include "TRandom3.h"
   7 #include "TFitResultPtr.h"
   8 #include "TFitResult.h"
   9
  10 #include "AliMCParticle.h"
  11
  12 #include "AliAnalysisTask.h"
  13 #include "AliAnalysisManager.h"
  14
  15 #include "AliESDEvent.h"
  16 #include "AliMCEvent.h"
  17 #include "AliESDInputHandler.h"
  18 #include "AliInputEventHandler.h"
  19
  20 #include "AliVVertex.h"
  21 #include "AliAnalysisFilter.h"
  22 #include "AliPID.h"
  23 #include "AliPIDCombined.h"
  24 #include "AliPIDResponse.h"
  25 #include "AliTPCPIDResponse.h"
  26
  27 #include "AliAnalysisTaskPID.h"
  28
  29 /*
  30 This task collects PID output from different detectors.
  31 Only tracks fulfilling some standard quality cuts are taken into account.
  32 At the moment, only data from TPC and TOF is collected. But in future,
  33 data from e.g. HMPID is also foreseen.
  34
  35 Contact: bhess@cern.ch
  36 */
  37
  38 ClassImp(AliAnalysisTaskPID)
  39
  40 const Int_t AliAnalysisTaskPID::fgkNumJetAxes = 3; // Number of additional axes for jets
  41 const Double_t AliAnalysisTaskPID::fgkEpsilon = 1e-8; // Double_t threshold above zero
  42 const Int_t AliAnalysisTaskPID::fgkMaxNumGenEntries = 500; // Maximum number of generated detector responses per track and delta(Prime) and associated species
  43
  44 const Double_t AliAnalysisTaskPID::fgkOneOverSqrt2 = 0.707106781186547462; // = 1. / TMath::Sqrt2();
  45
  46 const Double_t AliAnalysisTaskPID::fgkSigmaReferenceForTransitionPars = 0.05; // Reference sigma chosen to calculate transition
  47
  48 //________________________________________________________________________
  49 AliAnalysisTaskPID::AliAnalysisTaskPID()
  50   : AliAnalysisTaskPIDV0base()
  51   , fRun(-1)
  52   , fPIDcombined(new AliPIDCombined())
  53   , fInputFromOtherTask(kFALSE)
  54   , fDoPID(kTRUE)
  55   , fDoEfficiency(kTRUE)
  56   , fDoPtResolution(kTRUE)
  57   , fStoreCentralityPercentile(kFALSE)
  58   , fStoreAdditionalJetInformation(kFALSE)
  59   , fTakeIntoAccountMuons(kFALSE)
  60   , fUseITS(kFALSE)
  61   , fUseTOF(kFALSE)
  62   , fUsePriors(kFALSE)
  63   , fTPCDefaultPriors(kFALSE)
  64   , fUseMCidForGeneration(kTRUE)
  65   , fUseConvolutedGaus(kFALSE)
  66   , fkConvolutedGausNPar(3)
  67   , fAccuracyNonGaussianTail(1e-8)
  68   , fkDeltaPrimeLowLimit(0.02)
  69   , fkDeltaPrimeUpLimit(40.0)
  70   , fConvolutedGausDeltaPrime(0x0)
  71   , fTOFmode(1)
  72   , fEtaAbsCutLow(0.0)
  73   , fEtaAbsCutUp(0.9)
  74   , fDoAnySystematicStudiesOnTheExpectedSignal(kFALSE)
  75   , fSystematicScalingSplinesThreshold(50.)
  76   , fSystematicScalingSplinesBelowThreshold(1.0)
  77   , fSystematicScalingSplinesAboveThreshold(1.0)
  78   , fSystematicScalingEtaCorrectionMomentumThr(0.35)
  79   , fSystematicScalingEtaCorrectionLowMomenta(1.0)
  80   , fSystematicScalingEtaCorrectionHighMomenta(1.0)
  81   , fSystematicScalingEtaSigmaParaThreshold(250.)
  82   , fSystematicScalingEtaSigmaParaBelowThreshold(1.0)
  83   , fSystematicScalingEtaSigmaParaAboveThreshold(1.0)
  84   , fSystematicScalingMultCorrection(1.0)
  85   , fCentralityEstimator("V0M")
  86   , fhPIDdataAll(0x0)
  87   , fhGenEl(0x0)
  88   , fhGenKa(0x0)
  89   , fhGenPi(0x0)
  90   , fhGenMu(0x0)
  91   , fhGenPr(0x0)
  92   , fGenRespElDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
  93   , fGenRespElDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
  94   , fGenRespElDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
  95   , fGenRespElDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
  96   , fGenRespKaDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
  97   , fGenRespKaDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
  98   , fGenRespKaDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
  99   , fGenRespKaDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
 100   , fGenRespPiDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
 101   , fGenRespPiDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
 102   , fGenRespPiDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
 103   , fGenRespPiDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
 104   , fGenRespMuDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
 105   , fGenRespMuDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
 106   , fGenRespMuDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
 107   , fGenRespMuDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
 108   , fGenRespPrDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
 109   , fGenRespPrDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
 110   , fGenRespPrDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
 111   , fGenRespPrDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
 112   /*
 113   , fGenRespElDeltaEl(new Double_t[fgkMaxNumGenEntries])
 114   , fGenRespElDeltaKa(new Double_t[fgkMaxNumGenEntries])
 115   , fGenRespElDeltaPi(new Double_t[fgkMaxNumGenEntries])
 116   , fGenRespElDeltaPr(new Double_t[fgkMaxNumGenEntries])
 117   , fGenRespKaDeltaEl(new Double_t[fgkMaxNumGenEntries])
 118   , fGenRespKaDeltaKa(new Double_t[fgkMaxNumGenEntries])
 119   , fGenRespKaDeltaPi(new Double_t[fgkMaxNumGenEntries])
 120   , fGenRespKaDeltaPr(new Double_t[fgkMaxNumGenEntries])
 121   , fGenRespPiDeltaEl(new Double_t[fgkMaxNumGenEntries])
 122   , fGenRespPiDeltaKa(new Double_t[fgkMaxNumGenEntries])
 123   , fGenRespPiDeltaPi(new Double_t[fgkMaxNumGenEntries])
 124   , fGenRespPiDeltaPr(new Double_t[fgkMaxNumGenEntries])
 125   , fGenRespMuDeltaEl(new Double_t[fgkMaxNumGenEntries])
 126   , fGenRespMuDeltaKa(new Double_t[fgkMaxNumGenEntries])
 127   , fGenRespMuDeltaPi(new Double_t[fgkMaxNumGenEntries])
 128   , fGenRespMuDeltaPr(new Double_t[fgkMaxNumGenEntries])
 129   , fGenRespPrDeltaEl(new Double_t[fgkMaxNumGenEntries])
 130   , fGenRespPrDeltaKa(new Double_t[fgkMaxNumGenEntries])
 131   , fGenRespPrDeltaPi(new Double_t[fgkMaxNumGenEntries])
 132   , fGenRespPrDeltaPr(new Double_t[fgkMaxNumGenEntries])
 133   */
 134   , fhMaxEtaVariation(0x0)
 135   , fhEventsProcessed(0x0)
 136   , fhEventsTriggerSel(0x0)
 137   , fhEventsTriggerSelVtxCut(0x0)
 138   , fhSkippedTracksForSignalGeneration(0x0)
 139   , fhMCgeneratedYieldsPrimaries(0x0)
 140   , fh2FFJetPtRec(0x0)
 141   , fh2FFJetPtGen(0x0)
 142   , fh1Xsec(0x0)
 143   , fh1Trials(0x0)
 144   , fContainerEff(0x0)
 145   , fOutputContainer(0x0)
 146   , fPtResolutionContainer(0x0)
 147 {
 148   // default Constructor
 149
 150   AliLog::SetClassDebugLevel("AliAnalysisTaskPID", AliLog::kInfo);
 151
 152   fConvolutedGausDeltaPrime = new TF1("convolutedGausDeltaPrime", this, &AliAnalysisTaskPID::ConvolutedGaus,
 153                                       fkDeltaPrimeLowLimit, fkDeltaPrimeUpLimit,
 154                                       fkConvolutedGausNPar, "AliAnalysisTaskPID", "ConvolutedGaus");
 155
 156   // Set some arbitrary parameteres, such that the function call will not crash
 157   // (although it should not be called with these parameters...)
 158   fConvolutedGausDeltaPrime->SetParameter(0, 0);
 159   fConvolutedGausDeltaPrime->SetParameter(1, 1);
 160   fConvolutedGausDeltaPrime->SetParameter(2, 2);
 161
 162
 163   // Initialisation of translation parameters is time consuming.
 164   // Therefore, default values will only be initialised if they are really needed.
 165   // Otherwise, it is left to the user to set the parameter properly.
 166   fConvolutedGaussTransitionPars[0] = -999;
 167   fConvolutedGaussTransitionPars[1] = -999;
 168   fConvolutedGaussTransitionPars[2] = -999;
 169
 170   // Fraction histos
 171   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
 172     fFractionHists[i] = 0x0;
 173     fFractionSysErrorHists[i] = 0x0;
 174
 175     fPtResolution[i] = 0x0;
 176   }
 177 }
 178
 179 //________________________________________________________________________
 180 AliAnalysisTaskPID::AliAnalysisTaskPID(const char *name)
 181   : AliAnalysisTaskPIDV0base(name)
 182   , fRun(-1)
 183   , fPIDcombined(new AliPIDCombined())
 184   , fInputFromOtherTask(kFALSE)
 185   , fDoPID(kTRUE)
 186   , fDoEfficiency(kTRUE)
 187   , fDoPtResolution(kTRUE)
 188   , fStoreCentralityPercentile(kFALSE)
 189   , fStoreAdditionalJetInformation(kFALSE)
 190   , fTakeIntoAccountMuons(kFALSE)
 191   , fUseITS(kFALSE)
 192   , fUseTOF(kFALSE)
 193   , fUsePriors(kFALSE)
 194   , fTPCDefaultPriors(kFALSE)
 195   , fUseMCidForGeneration(kTRUE)
 196   , fUseConvolutedGaus(kFALSE)
 197   , fkConvolutedGausNPar(3)
 198   , fAccuracyNonGaussianTail(1e-8)
 199   , fkDeltaPrimeLowLimit(0.02)
 200   , fkDeltaPrimeUpLimit(40.0)
 201   , fConvolutedGausDeltaPrime(0x0)
 202   , fTOFmode(1)
 203   , fEtaAbsCutLow(0.0)
 204   , fEtaAbsCutUp(0.9)
 205   , fDoAnySystematicStudiesOnTheExpectedSignal(kFALSE)
 206   , fSystematicScalingSplinesThreshold(50.)
 207   , fSystematicScalingSplinesBelowThreshold(1.0)
 208   , fSystematicScalingSplinesAboveThreshold(1.0)
 209   , fSystematicScalingEtaCorrectionMomentumThr(0.35)
 210   , fSystematicScalingEtaCorrectionLowMomenta(1.0)
 211   , fSystematicScalingEtaCorrectionHighMomenta(1.0)
 212   , fSystematicScalingEtaSigmaParaThreshold(250.)
 213   , fSystematicScalingEtaSigmaParaBelowThreshold(1.0)
 214   , fSystematicScalingEtaSigmaParaAboveThreshold(1.0)
 215   , fSystematicScalingMultCorrection(1.0)
 216   , fCentralityEstimator("V0M")
 217   , fhPIDdataAll(0x0)
 218   , fhGenEl(0x0)
 219   , fhGenKa(0x0)
 220   , fhGenPi(0x0)
 221   , fhGenMu(0x0)
 222   , fhGenPr(0x0)
 223   , fGenRespElDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
 224   , fGenRespElDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
 225   , fGenRespElDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
 226   , fGenRespElDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
 227   , fGenRespKaDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
 228   , fGenRespKaDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
 229   , fGenRespKaDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
 230   , fGenRespKaDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
 231   , fGenRespPiDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
 232   , fGenRespPiDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
 233   , fGenRespPiDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
 234   , fGenRespPiDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
 235   , fGenRespMuDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
 236   , fGenRespMuDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
 237   , fGenRespMuDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
 238   , fGenRespMuDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
 239   , fGenRespPrDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
 240   , fGenRespPrDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
 241   , fGenRespPrDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
 242   , fGenRespPrDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
 243   /*
 244   , fGenRespElDeltaEl(new Double_t[fgkMaxNumGenEntries])
 245   , fGenRespElDeltaKa(new Double_t[fgkMaxNumGenEntries])
 246   , fGenRespElDeltaPi(new Double_t[fgkMaxNumGenEntries])
 247   , fGenRespElDeltaPr(new Double_t[fgkMaxNumGenEntries])
 248   , fGenRespKaDeltaEl(new Double_t[fgkMaxNumGenEntries])
 249   , fGenRespKaDeltaKa(new Double_t[fgkMaxNumGenEntries])
 250   , fGenRespKaDeltaPi(new Double_t[fgkMaxNumGenEntries])
 251   , fGenRespKaDeltaPr(new Double_t[fgkMaxNumGenEntries])
 252   , fGenRespPiDeltaEl(new Double_t[fgkMaxNumGenEntries])
 253   , fGenRespPiDeltaKa(new Double_t[fgkMaxNumGenEntries])
 254   , fGenRespPiDeltaPi(new Double_t[fgkMaxNumGenEntries])
 255   , fGenRespPiDeltaPr(new Double_t[fgkMaxNumGenEntries])
 256   , fGenRespMuDeltaEl(new Double_t[fgkMaxNumGenEntries])
 257   , fGenRespMuDeltaKa(new Double_t[fgkMaxNumGenEntries])
 258   , fGenRespMuDeltaPi(new Double_t[fgkMaxNumGenEntries])
 259   , fGenRespMuDeltaPr(new Double_t[fgkMaxNumGenEntries])
 260   , fGenRespPrDeltaEl(new Double_t[fgkMaxNumGenEntries])
 261   , fGenRespPrDeltaKa(new Double_t[fgkMaxNumGenEntries])
 262   , fGenRespPrDeltaPi(new Double_t[fgkMaxNumGenEntries])
 263   , fGenRespPrDeltaPr(new Double_t[fgkMaxNumGenEntries])
 264   */
 265   , fhMaxEtaVariation(0x0)
 266   , fhEventsProcessed(0x0)
 267   , fhEventsTriggerSel(0x0)
 268   , fhEventsTriggerSelVtxCut(0x0)
 269   , fhSkippedTracksForSignalGeneration(0x0)
 270   , fhMCgeneratedYieldsPrimaries(0x0)
 271   , fh2FFJetPtRec(0x0)
 272   , fh2FFJetPtGen(0x0)
 273   , fh1Xsec(0x0)
 274   , fh1Trials(0x0)
 275   , fContainerEff(0x0)
 276   , fOutputContainer(0x0)
 277   , fPtResolutionContainer(0x0)
 278 {
 279   // Constructor
 280
 281   AliLog::SetClassDebugLevel("AliAnalysisTaskPID", AliLog::kInfo);
 282
 283   fConvolutedGausDeltaPrime = new TF1("convolutedGausDeltaPrime", this, &AliAnalysisTaskPID::ConvolutedGaus,
 284                                       fkDeltaPrimeLowLimit, fkDeltaPrimeUpLimit,
 285                                       fkConvolutedGausNPar, "AliAnalysisTaskPID", "ConvolutedGaus");
 286
 287   // Set some arbitrary parameteres, such that the function call will not crash
 288   // (although it should not be called with these parameters...)
 289   fConvolutedGausDeltaPrime->SetParameter(0, 0);
 290   fConvolutedGausDeltaPrime->SetParameter(1, 1);
 291   fConvolutedGausDeltaPrime->SetParameter(2, 2);
 292
 293
 294   // Initialisation of translation parameters is time consuming.
 295   // Therefore, default values will only be initialised if they are really needed.
 296   // Otherwise, it is left to the user to set the parameter properly.
 297   fConvolutedGaussTransitionPars[0] = -999;
 298   fConvolutedGaussTransitionPars[1] = -999;
 299   fConvolutedGaussTransitionPars[2] = -999;
 300
 301   // Fraction histos
 302   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
 303     fFractionHists[i] = 0x0;
 304     fFractionSysErrorHists[i] = 0x0;
 305
 306     fPtResolution[i] = 0x0;
 307   }
 308
 309   // Define input and output slots here
 310   // Input slot #0 works with a TChain
 311   DefineInput(0, TChain::Class());
 312
 313   DefineOutput(1, TObjArray::Class());
 314
 315   DefineOutput(2, AliCFContainer::Class());
 316
 317   DefineOutput(3, TObjArray::Class());
 318 }
 319
 320
 321 //________________________________________________________________________
 322 AliAnalysisTaskPID::~AliAnalysisTaskPID()
 323 {
 324   // dtor
 325
 326   CleanupParticleFractionHistos();
 327
 328   delete fOutputContainer;
 329   fOutputContainer = 0x0;
 330
 331   delete fPtResolutionContainer;
 332   fPtResolutionContainer = 0x0;
 333
 334   delete fConvolutedGausDeltaPrime;
 335   fConvolutedGausDeltaPrime = 0x0;
 336
 337   delete [] fGenRespElDeltaPrimeEl;
 338   delete [] fGenRespElDeltaPrimeKa;
 339   delete [] fGenRespElDeltaPrimePi;
 340   delete [] fGenRespElDeltaPrimePr;
 341
 342   fGenRespElDeltaPrimeEl = 0x0;
 343   fGenRespElDeltaPrimeKa = 0x0;
 344   fGenRespElDeltaPrimePi = 0x0;
 345   fGenRespElDeltaPrimePr = 0x0;
 346
 347   delete [] fGenRespKaDeltaPrimeEl;
 348   delete [] fGenRespKaDeltaPrimeKa;
 349   delete [] fGenRespKaDeltaPrimePi;
 350   delete [] fGenRespKaDeltaPrimePr;
 351
 352   fGenRespKaDeltaPrimeEl = 0x0;
 353   fGenRespKaDeltaPrimeKa = 0x0;
 354   fGenRespKaDeltaPrimePi = 0x0;
 355   fGenRespKaDeltaPrimePr = 0x0;
 356
 357   delete [] fGenRespPiDeltaPrimeEl;
 358   delete [] fGenRespPiDeltaPrimeKa;
 359   delete [] fGenRespPiDeltaPrimePi;
 360   delete [] fGenRespPiDeltaPrimePr;
 361
 362   fGenRespPiDeltaPrimeEl = 0x0;
 363   fGenRespPiDeltaPrimeKa = 0x0;
 364   fGenRespPiDeltaPrimePi = 0x0;
 365   fGenRespPiDeltaPrimePr = 0x0;
 366
 367   delete [] fGenRespMuDeltaPrimeEl;
 368   delete [] fGenRespMuDeltaPrimeKa;
 369   delete [] fGenRespMuDeltaPrimePi;
 370   delete [] fGenRespMuDeltaPrimePr;
 371
 372   fGenRespMuDeltaPrimeEl = 0x0;
 373   fGenRespMuDeltaPrimeKa = 0x0;
 374   fGenRespMuDeltaPrimePi = 0x0;
 375   fGenRespMuDeltaPrimePr = 0x0;
 376
 377   delete [] fGenRespPrDeltaPrimeEl;
 378   delete [] fGenRespPrDeltaPrimeKa;
 379   delete [] fGenRespPrDeltaPrimePi;
 380   delete [] fGenRespPrDeltaPrimePr;
 381
 382   fGenRespPrDeltaPrimeEl = 0x0;
 383   fGenRespPrDeltaPrimeKa = 0x0;
 384   fGenRespPrDeltaPrimePi = 0x0;
 385   fGenRespPrDeltaPrimePr = 0x0;
 386
 387   delete fhMaxEtaVariation;
 388   fhMaxEtaVariation = 0x0;
 389
 390   /*OLD with deltaSpecies
 391   delete [] fGenRespElDeltaEl;
 392   delete [] fGenRespElDeltaKa;
 393   delete [] fGenRespElDeltaPi;
 394   delete [] fGenRespElDeltaPr;
 395
 396   fGenRespElDeltaEl = 0x0;
 397   fGenRespElDeltaKa = 0x0;
 398   fGenRespElDeltaPi = 0x0;
 399   fGenRespElDeltaPr = 0x0;
 400
 401   delete [] fGenRespKaDeltaEl;
 402   delete [] fGenRespKaDeltaKa;
 403   delete [] fGenRespKaDeltaPi;
 404   delete [] fGenRespKaDeltaPr;
 405
 406   fGenRespKaDeltaEl = 0x0;
 407   fGenRespKaDeltaKa = 0x0;
 408   fGenRespKaDeltaPi = 0x0;
 409   fGenRespKaDeltaPr = 0x0;
 410
 411   delete [] fGenRespPiDeltaEl;
 412   delete [] fGenRespPiDeltaKa;
 413   delete [] fGenRespPiDeltaPi;
 414   delete [] fGenRespPiDeltaPr;
 415
 416   fGenRespPiDeltaEl = 0x0;
 417   fGenRespPiDeltaKa = 0x0;
 418   fGenRespPiDeltaPi = 0x0;
 419   fGenRespPiDeltaPr = 0x0;
 420
 421   delete [] fGenRespMuDeltaEl;
 422   delete [] fGenRespMuDeltaKa;
 423   delete [] fGenRespMuDeltaPi;
 424   delete [] fGenRespMuDeltaPr;
 425
 426   fGenRespMuDeltaEl = 0x0;
 427   fGenRespMuDeltaKa = 0x0;
 428   fGenRespMuDeltaPi = 0x0;
 429   fGenRespMuDeltaPr = 0x0;
 430
 431   delete [] fGenRespPrDeltaEl;
 432   delete [] fGenRespPrDeltaKa;
 433   delete [] fGenRespPrDeltaPi;
 434   delete [] fGenRespPrDeltaPr;
 435
 436   fGenRespPrDeltaEl = 0x0;
 437   fGenRespPrDeltaKa = 0x0;
 438   fGenRespPrDeltaPi = 0x0;
 439   fGenRespPrDeltaPr = 0x0;
 440   */
 441 }
 442
 443
 444 //________________________________________________________________________
 445 void AliAnalysisTaskPID::SetUpPIDcombined()
 446 {
 447   // Initialise the PIDcombined object
 448
 449   if (!fDoPID)
 450     return;
 451
 452   if(fDebug > 1)
 453     printf("File: %s, Line: %d: SetUpPIDcombined\n", (char*)__FILE__, __LINE__);
 454
 455   if (!fPIDcombined) {
 456     AliFatal("No PIDcombined object!\n");
 457     return;
 458   }
 459
 460   fPIDcombined->SetSelectedSpecies(AliPID::kSPECIESC);
 461   fPIDcombined->SetEnablePriors(fUsePriors);
 462
 463   if (fTPCDefaultPriors)
 464     fPIDcombined->SetDefaultTPCPriors();
 465
 466   //TODO use individual priors...
 467
 468   // Change detector mask (TPC,TOF,ITS)
 469   Int_t detectorMask = AliPIDResponse::kDetTPC;
 470
 471   // Reject mismatch mask - mismatch only relevant for TOF at the moment - other detectors do not use it
 472   Int_t rejectMismatchMask = AliPIDResponse::kDetTPC;
 473
 474
 475   if (fUseITS) {
 476     detectorMask = detectorMask | AliPIDResponse::kDetITS;
 477     rejectMismatchMask = rejectMismatchMask | AliPIDResponse::kDetITS;
 478   }
 479   if (fUseTOF) {
 480     detectorMask = detectorMask | AliPIDResponse::kDetTOF;
 481     rejectMismatchMask = rejectMismatchMask | AliPIDResponse::kDetTOF;
 482   }
 483
 484   fPIDcombined->SetDetectorMask(detectorMask);
 485   fPIDcombined->SetRejectMismatchMask(rejectMismatchMask);
 486
 487   if(fDebug > 1)
 488     printf("File: %s, Line: %d: SetUpPIDcombined done\n", (char*)__FILE__, __LINE__);
 489 }
 490
 491
 492 //________________________________________________________________________
 493 Bool_t AliAnalysisTaskPID::CalculateMaxEtaVariationMapFromPIDResponse()
 494 {
 495   // Calculate the maximum deviation from unity of the eta correction factors for each row in 1/dEdx(splines)
 496   // from the eta correction map of the TPCPIDResponse. The result is stored in fhMaxEtaVariation.
 497
 498   if (!fPIDResponse) {
 499     AliError("No PID response!");
 500     return kFALSE;
 501   }
 502
 503   delete fhMaxEtaVariation;
 504
 505   const TH2D* hEta = fPIDResponse->GetTPCResponse().GetEtaCorrMap();
 506   if (!hEta) {
 507     AliError("No eta correction map!");
 508     return kFALSE;
 509   }
 510
 511   // Take binning from hEta in Y for fhMaxEtaVariation
 512   fhMaxEtaVariation = hEta->ProjectionY("hMaxEtaVariation");
 513   fhMaxEtaVariation->SetDirectory(0);
 514   fhMaxEtaVariation->Reset();
 515
 516   // For each bin in 1/dEdx, loop of all tanTheta bins and find the maximum deviation from unity.
 517   // Store the result in fhMaxEtaVariation
 518
 519   for (Int_t binY = 1; binY <= fhMaxEtaVariation->GetNbinsX(); binY++) {
 520     Double_t maxAbs = -1;
 521     for (Int_t binX = 1; binX <= hEta->GetNbinsX(); binX++) {
 522       Double_t curr = TMath::Abs(hEta->GetBinContent(binX, binY) - 1.);
 523       if (curr > maxAbs)
 524         maxAbs = curr;
 525     }
 526
 527     if (maxAbs < 1e-12) {
 528       AliError(Form("Maximum deviation from unity is zero for 1/dEdx = %f (bin %d)", hEta->GetYaxis()->GetBinCenter(binY), binY));
 529       delete fhMaxEtaVariation;
 530       return kFALSE;
 531     }
 532
 533     fhMaxEtaVariation->SetBinContent(binY, maxAbs);
 534   }
 535
 536   printf("AliAnalysisTaskPID: Calculated max eta variation.\n");
 537
 538   return kTRUE;
 539 }
 540
 541
 542 //________________________________________________________________________
 543 void AliAnalysisTaskPID::UserCreateOutputObjects()
 544 {
 545   // Create histograms
 546   // Called once
 547
 548   if(fDebug > 1)
 549     printf("File: %s, Line: %d: UserCreateOutputObjects\n", (char*)__FILE__, __LINE__);
 550
 551   SetUpPIDcombined();
 552
 553   // Input handler
 554   AliAnalysisManager* man = AliAnalysisManager::GetAnalysisManager();
 555   AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
 556
 557   if (!inputHandler)
 558     AliFatal("Input handler needed");
 559   else {
 560     // PID response object
 561     fPIDResponse = inputHandler->GetPIDResponse();
 562     if (!fPIDResponse)
 563       AliFatal("PIDResponse object was not created");
 564   }
 565
 566   if(fDebug > 2)
 567     printf("File: %s, Line: %d: UserCreateOutputObjects -> Retrieved PIDresponse object\n", (char*)__FILE__, __LINE__);
 568
 569   OpenFile(1);
 570
 571   if(fDebug > 2)
 572     printf("File: %s, Line: %d: UserCreateOutputObjects -> OpenFile(1) successful\n", (char*)__FILE__, __LINE__);
 573
 574   fOutputContainer = new TObjArray(1);
 575   fOutputContainer->SetName(GetName());
 576   fOutputContainer->SetOwner(kTRUE);
 577
 578   const Int_t nPtBins = 68;
 579   Double_t binsPt[nPtBins+1] = {0. ,  0.05, 0.1,  0.15, 0.2,  0.25, 0.3,  0.35, 0.4,  0.45,
 580            0.5,  0.55, 0.6,  0.65, 0.7,  0.75, 0.8,  0.85, 0.9,  0.95,
 581            1.0,  1.1 , 1.2,  1.3 , 1.4,  1.5 , 1.6,  1.7 , 1.8,  1.9 ,
 582            2.0,  2.2 , 2.4,  2.6 , 2.8,  3.0 , 3.2,  3.4 , 3.6,  3.8 ,
 583            4.0,  4.5 , 5.0,  5.5 , 6.0,  6.5 , 7.0,  8.0 , 9.0,  10.0,
 584            11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0,
 585            26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 };
 586
 587   const Int_t nCentBins = 12;
 588   //-1 for pp; 90-100 has huge electro-magnetic impurities
 589   Double_t binsCent[nCentBins+1] = {-1, 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 };
 590
 591   const Int_t nJetPtBins = 11;
 592   Double_t binsJetPt[nJetPtBins+1] = {0, 2, 5, 10, 15, 20, 30, 40, 60, 80, 120, 200};
 593
 594   const Int_t nChargeBins = 2;
 595   const Double_t binsCharge[nChargeBins+1] = { -1.0 - 1e-4, 0.0, 1.0 + 1e-4 };
 596
 597   const Int_t nBinsJets = kDataNumAxes;
 598   const Int_t nBinsNoJets = nBinsJets - fgkNumJetAxes;
 599
 600   const Int_t nBins = fStoreAdditionalJetInformation ? nBinsJets : nBinsNoJets;
 601
 602   // deltaPrimeSpecies binning
 603   const Int_t deltaPrimeNBins = 600;
 604   Double_t deltaPrimeBins[deltaPrimeNBins + 1];
 605
 606   const Double_t fromLow = fkDeltaPrimeLowLimit;
 607   const Double_t toHigh = fkDeltaPrimeUpLimit;
 608   const Double_t factor = TMath::Power(toHigh/fromLow, 1./deltaPrimeNBins);
 609
 610   // Log binning for whole deltaPrime range
 611   deltaPrimeBins[0] = fromLow;
 612   for (Int_t i = 0 + 1; i <= deltaPrimeNBins; i++) {
 613     deltaPrimeBins[i] = factor * deltaPrimeBins[i - 1];
 614   }
 615
 616   const Int_t nMCPIDbins = 5;
 617   const Double_t mcPIDmin = 0.;
 618   const Double_t mcPIDmax = 5.;
 619
 620   const Int_t nSelSpeciesBins = 4;
 621   const Double_t selSpeciesMin = 0.;
 622   const Double_t selSpeciesMax = 4.;
 623
 624   const Int_t nZBins = 20;
 625   const Double_t zMin = 0.;
 626   const Double_t zMax = 1.;
 627
 628   const Int_t nXiBins = 70;
 629   const Double_t xiMin = 0.;
 630   const Double_t xiMax = 7.;
 631
 632   const Int_t nTOFpidInfoBins = kNumTOFpidInfoBins;
 633   const Double_t tofPIDinfoMin = kNoTOFinfo;
 634   const Double_t tofPIDinfoMax = kNoTOFinfo + kNumTOFpidInfoBins;
 635
 636   // MC PID, SelectSpecies, pT, deltaPrimeSpecies, centrality percentile, jet pT, z = track_pT/jet_pT, xi = log(1/z)
 637   Int_t binsNoJets[nBinsNoJets] =    { nMCPIDbins,
 638                                        nSelSpeciesBins,
 639                                        nPtBins,
 640                                        deltaPrimeNBins,
 641                                        nCentBins,
 642                                        nChargeBins,
 643                                        nTOFpidInfoBins };
 644
 645   Int_t binsJets[nBinsJets]     =    { nMCPIDbins,
 646                                        nSelSpeciesBins,
 647                                        nPtBins,
 648                                        deltaPrimeNBins,
 649                                        nCentBins,
 650                                        nJetPtBins,
 651                                        nZBins,
 652                                        nXiBins,
 653                                        nChargeBins,
 654                                        nTOFpidInfoBins };
 655
 656   Int_t *bins = fStoreAdditionalJetInformation ? &binsJets[0] : &binsNoJets[0];
 657
 658   Double_t xminNoJets[nBinsNoJets] = { mcPIDmin,
 659                                        selSpeciesMin,
 660                                        binsPt[0],
 661                                        deltaPrimeBins[0],
 662                                        binsCent[0],
 663                                        binsCharge[0],
 664                                        tofPIDinfoMin };
 665
 666   Double_t xminJets[nBinsJets] =     { mcPIDmin,
 667                                        selSpeciesMin,
 668                                        binsPt[0],
 669                                        deltaPrimeBins[0],
 670                                        binsCent[0],
 671                                        binsJetPt[0],
 672                                        zMin,
 673                                        xiMin,
 674                                        binsCharge[0],
 675                                        tofPIDinfoMin };
 676
 677   Double_t *xmin = fStoreAdditionalJetInformation? &xminJets[0] : &xminNoJets[0];
 678
 679   Double_t xmaxNoJets[nBinsNoJets] = { mcPIDmax,
 680                                        selSpeciesMax,
 681                                        binsPt[nPtBins],
 682                                        deltaPrimeBins[deltaPrimeNBins],
 683                                        binsCent[nCentBins],
 684                                        binsCharge[nChargeBins],
 685                                        tofPIDinfoMax };
 686
 687   Double_t xmaxJets[nBinsJets] =     { mcPIDmax,
 688                                        selSpeciesMax,
 689                                        binsPt[nPtBins],
 690                                        deltaPrimeBins[deltaPrimeNBins],
 691                                        binsCent[nCentBins],
 692                                        binsJetPt[nJetPtBins],
 693                                        zMax,
 694                                        xiMax,
 695                                        binsCharge[nChargeBins],
 696                                        tofPIDinfoMax };
 697
 698   Double_t *xmax = fStoreAdditionalJetInformation? &xmaxJets[0] : &xmaxNoJets[0];
 699
 700   fConvolutedGausDeltaPrime->SetNpx(deltaPrimeNBins);
 701
 702   if (fDoPID) {
 703     fhPIDdataAll = new THnSparseD("hPIDdataAll","", nBins, bins, xmin, xmax);
 704     SetUpHist(fhPIDdataAll, binsPt, deltaPrimeBins, binsCent, binsJetPt);
 705     fOutputContainer->Add(fhPIDdataAll);
 706   }
 707
 708   // Generated histograms (so far, bins are the same as for primary THnSparse)
 709   const Int_t nGenBins = fStoreAdditionalJetInformation ? nBinsJets : nBinsNoJets;
 710   // MC PID, SelectSpecies, Pt, deltaPrimeSpecies, jet pT, z = track_pT/jet_pT, xi = log(1/z)
 711
 712   Int_t *genBins = fStoreAdditionalJetInformation ? &binsJets[0] : &binsNoJets[0];
 713   Double_t *genXmin = fStoreAdditionalJetInformation? &xminJets[0] : &xminNoJets[0];
 714   Double_t *genXmax = fStoreAdditionalJetInformation? &xmaxJets[0] : &xmaxNoJets[0];
 715
 716   if (fDoPID) {
 717     fhGenEl = new THnSparseD("hGenEl", "", nGenBins, genBins, genXmin, genXmax);
 718     SetUpGenHist(fhGenEl, binsPt, deltaPrimeBins, binsCent, binsJetPt);
 719     fOutputContainer->Add(fhGenEl);
 720
 721     fhGenKa = new THnSparseD("hGenKa", "", nGenBins, genBins, genXmin, genXmax);
 722     SetUpGenHist(fhGenKa, binsPt, deltaPrimeBins, binsCent, binsJetPt);
 723     fOutputContainer->Add(fhGenKa);
 724
 725     fhGenPi = new THnSparseD("hGenPi", "", nGenBins, genBins, genXmin, genXmax);
 726     SetUpGenHist(fhGenPi, binsPt, deltaPrimeBins, binsCent, binsJetPt);
 727     fOutputContainer->Add(fhGenPi);
 728
 729     if (fTakeIntoAccountMuons) {
 730       fhGenMu = new THnSparseD("hGenMu", "", nGenBins, genBins, genXmin, genXmax);
 731       SetUpGenHist(fhGenMu, binsPt, deltaPrimeBins, binsCent, binsJetPt);
 732       fOutputContainer->Add(fhGenMu);
 733     }
 734
 735     fhGenPr = new THnSparseD("hGenPr", "", nGenBins, genBins, genXmin, genXmax);
 736     SetUpGenHist(fhGenPr, binsPt, deltaPrimeBins, binsCent, binsJetPt);
 737     fOutputContainer->Add(fhGenPr);
 738
 739     fhSkippedTracksForSignalGeneration = new TH2D("fhSkippedTracksForSignalGeneration",
 740                                                   "Number of tracks skipped for the signal generation;p_{T}^{gen} (GeV/c);TPC signal N",
 741                                                   nPtBins, binsPt, 161, -0.5, 160.5);
 742     fhSkippedTracksForSignalGeneration->Sumw2();
 743     fOutputContainer->Add(fhSkippedTracksForSignalGeneration);
 744   }
 745
 746
 747   fhEventsProcessed = new TH1D("fhEventsProcessed",
 748                                "Number of events passing trigger selection, vtx and zvtx cuts;Centrality percentile",
 749                                nCentBins, binsCent);
 750   fhEventsProcessed->Sumw2();
 751   fOutputContainer->Add(fhEventsProcessed);
 752
 753   fhEventsTriggerSelVtxCut = new TH1D("fhEventsTriggerSelVtxCut",
 754                                       "Number of events passing trigger selection and vtx cut;Centrality percentile",
 755                                       nCentBins, binsCent);
 756   fhEventsTriggerSelVtxCut->Sumw2();
 757   fOutputContainer->Add(fhEventsTriggerSelVtxCut);
 758
 759   fhEventsTriggerSel = new TH1D("fhEventsTriggerSel",
 760                                 "Number of events passing trigger selection;Centrality percentile",
 761                                 nCentBins, binsCent);
 762   fOutputContainer->Add(fhEventsTriggerSel);
 763   fhEventsTriggerSel->Sumw2();
 764
 765
 766   // Generated yields within acceptance
 767   const Int_t nBinsGenYields = fStoreAdditionalJetInformation ? kGenYieldNumAxes : kGenYieldNumAxes - 3;
 768   Int_t genYieldsBins[kGenYieldNumAxes]    = { nMCPIDbins,         nPtBins,           nCentBins,            nJetPtBins, nZBins, nXiBins,
 769                                                nChargeBins };
 770   genYieldsBins[GetIndexOfChargeAxisGenYield()] = nChargeBins;
 771   Double_t genYieldsXmin[kGenYieldNumAxes] = {   mcPIDmin,       binsPt[0],         binsCent[0],          binsJetPt[0],   zMin,   xiMin,
 772                                                binsCharge[0] };
 773   genYieldsXmin[GetIndexOfChargeAxisGenYield()] = binsCharge[0];
 774   Double_t genYieldsXmax[kGenYieldNumAxes] = {   mcPIDmax, binsPt[nPtBins], binsCent[nCentBins], binsJetPt[nJetPtBins],   zMax,   xiMax,
 775                                                binsCharge[nChargeBins] };
 776   genYieldsXmax[GetIndexOfChargeAxisGenYield()] = binsCharge[nChargeBins];
 777
 778   if (fDoPID) {
 779     fhMCgeneratedYieldsPrimaries = new THnSparseD("fhMCgeneratedYieldsPrimaries",
 780                                                   "Generated yields w/o reco and cuts inside acceptance (physical primaries)",
 781                                                   nBinsGenYields, genYieldsBins, genYieldsXmin, genYieldsXmax);
 782     SetUpGenYieldHist(fhMCgeneratedYieldsPrimaries, binsPt, binsCent, binsJetPt);
 783     fOutputContainer->Add(fhMCgeneratedYieldsPrimaries);
 784   }
 785
 786   // Container with several process steps (generated and reconstructed level with some variations)
 787   if (fDoEfficiency) {
 788     OpenFile(2);
 789
 790     if(fDebug > 2)
 791     printf("File: %s, Line: %d: UserCreateOutputObjects -> OpenFile(2) successful\n", (char*)__FILE__, __LINE__);
 792
 793     // Array for the number of bins in each dimension
 794     // Dimensions: MC-ID, trackPt, trackEta, trackCharge, cenrality percentile, jetPt, z, xi TODO phi???
 795     const Int_t nEffDims = fStoreAdditionalJetInformation ? kEffNumAxes : kEffNumAxes - 3; // Number of dimensions for the efficiency
 796
 797     const Int_t nMCIDbins = AliPID::kSPECIES;
 798     Double_t binsMCID[nMCIDbins + 1];
 799
 800     for(Int_t i = 0; i <= nMCIDbins; i++) {
 801       binsMCID[i]= i;
 802     }
 803
 804     const Int_t nEtaBins = 18;
 805     const Double_t binsEta[nEtaBins+1] = {-0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1,
 806                                           0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 };
 807
 808     const Int_t nEffBins[kEffNumAxes] = { nMCIDbins, nPtBins, nEtaBins, nChargeBins, nCentBins, nJetPtBins, nZBins, nXiBins };
 809
 810     fContainerEff = new AliCFContainer("containerEff", "Reconstruction Efficiency x Acceptance x Resolution and Secondary Correction",
 811                                       kNumSteps, nEffDims, nEffBins);
 812
 813     // Setting the bin limits
 814     fContainerEff->SetBinLimits(kEffMCID, binsMCID);
 815     fContainerEff->SetBinLimits(kEffTrackPt, binsPt);
 816     fContainerEff->SetBinLimits(kEffTrackEta, binsEta);
 817     fContainerEff->SetBinLimits(kEffTrackCharge, binsCharge);
 818     fContainerEff->SetBinLimits(kEffCentrality, binsCent);
 819     if (fStoreAdditionalJetInformation) {
 820       fContainerEff->SetBinLimits(kEffJetPt, binsJetPt);
 821       fContainerEff->SetBinLimits(kEffZ, zMin, zMax);
 822       fContainerEff->SetBinLimits(kEffXi, xiMin, xiMax);
 823     }
 824
 825     fContainerEff->SetVarTitle(kEffMCID,"MC ID");
 826     fContainerEff->SetVarTitle(kEffTrackPt,"p_{T} (GeV/c)");
 827     fContainerEff->SetVarTitle(kEffTrackEta,"#eta");
 828     fContainerEff->SetVarTitle(kEffTrackCharge,"Charge (e_{0})");
 829     fContainerEff->SetVarTitle(kEffCentrality, "Centrality Percentile");
 830     if (fStoreAdditionalJetInformation) {
 831       fContainerEff->SetVarTitle(kEffJetPt, "p_{T}^{jet} (GeV/c)");
 832       fContainerEff->SetVarTitle(kEffZ, "z = p_{T}^{track} / p_{T}^{jet}");
 833       fContainerEff->SetVarTitle(kEffXi, "#xi = ln(p_{T}^{jet} / p_{T}^{track})");
 834     }
 835
 836     // Define clean MC sample
 837     fContainerEff->SetStepTitle(kStepGenWithGenCuts, "Particle level, cuts on particle level");
 838     // For Acceptance x Efficiency correction of primaries
 839     fContainerEff->SetStepTitle(kStepRecWithGenCuts, "Detector level (rec) with cuts on particle level");
 840     // For (pT) resolution correction
 841     fContainerEff->SetStepTitle(kStepRecWithGenCutsMeasuredObs,
 842                                 "Detector level (rec) with cuts on particle level with measured observables");
 843     // For secondary correction
 844     fContainerEff->SetStepTitle(kStepRecWithRecCutsMeasuredObs,
 845                                 "Detector level, all cuts on detector level with measured observables");
 846     fContainerEff->SetStepTitle(kStepRecWithRecCutsPrimaries,
 847                                 "Detector level, all cuts on detector level, only MC primaries");
 848     fContainerEff->SetStepTitle(kStepRecWithRecCutsMeasuredObsPrimaries,
 849                                 "Detector level, all cuts on detector level with measured observables, only MC primaries");
 850     fContainerEff->SetStepTitle(kStepRecWithRecCutsMeasuredObsStrangenessScaled,
 851                                 "Detector level (strangeness scaled), all cuts on detector level with measured observables");
 852   }
 853
 854   if (fDoPID || fDoEfficiency) {
 855     // Generated jets
 856     fh2FFJetPtRec = new TH2D("fh2FFJetPtRec", "Number of reconstructed jets;Centrality Percentile;p_{T}^{jet} (GeV/c)",
 857                             nCentBins, binsCent, nJetPtBins, binsJetPt);
 858     fh2FFJetPtRec->Sumw2();
 859     fOutputContainer->Add(fh2FFJetPtRec);
 860     fh2FFJetPtGen = new TH2D("fh2FFJetPtGen", "Number of generated jets;Centrality Percentile;p_{T}^{jet} (GeV/c)",
 861                             nCentBins, binsCent, nJetPtBins, binsJetPt);
 862     fh2FFJetPtGen->Sumw2();
 863     fOutputContainer->Add(fh2FFJetPtGen);
 864   }
 865
 866   // Pythia information
 867   fh1Xsec = new TProfile("fh1Xsec", "xsec from pyxsec.root", 1, 0, 1);
 868   fh1Xsec->Sumw2();
 869   fh1Xsec->GetXaxis()->SetBinLabel(1, "<#sigma>");
 870   fh1Trials = new TH1D("fh1Trials", "trials from pyxsec.root", 1, 0, 1);
 871   fh1Trials->Sumw2();
 872   fh1Trials->GetXaxis()->SetBinLabel(1, "#sum{ntrials}");
 873
 874   fOutputContainer->Add(fh1Xsec);
 875   fOutputContainer->Add(fh1Trials);
 876
 877   if (fDoPtResolution) {
 878     OpenFile(3);
 879
 880     if(fDebug > 2)
 881     printf("File: %s, Line: %d: UserCreateOutputObjects -> OpenFile(3) successful\n", (char*)__FILE__, __LINE__);
 882
 883     fPtResolutionContainer = new TObjArray(1);
 884     fPtResolutionContainer->SetName(Form("%s_PtResolution", GetName()));
 885     fPtResolutionContainer->SetOwner(kTRUE);
 886
 887     const Int_t nPtBinsRes = 100;
 888     Double_t pTbinsRes[nPtBinsRes + 1];
 889
 890     const Double_t fromLowPtRes = 0.15;
 891     const Double_t toHighPtRes = 50.;
 892     const Double_t factorPtRes = TMath::Power(toHighPtRes/fromLowPtRes, 1./nPtBinsRes);
 893     // Log binning for whole pT range
 894     pTbinsRes[0] = fromLowPtRes;
 895     for (Int_t i = 0 + 1; i <= nPtBinsRes; i++) {
 896       pTbinsRes[i] = factorPtRes * pTbinsRes[i - 1];
 897     }
 898
 899     const Int_t nBinsPtResolution = kPtResNumAxes;
 900     Int_t ptResolutionBins[kPtResNumAxes]    = { nJetPtBins,                       nPtBinsRes,            nPtBinsRes,
 901                                                  nChargeBins,             nCentBins };
 902     Double_t ptResolutionXmin[kPtResNumAxes] = { binsJetPt[0],                   pTbinsRes[0],          pTbinsRes[0],
 903                                                  binsCharge[0],           binsCent[0] };
 904     Double_t ptResolutionXmax[kPtResNumAxes] = { binsJetPt[nJetPtBins], pTbinsRes[nPtBinsRes], pTbinsRes[nPtBinsRes],
 905                                                  binsCharge[nChargeBins], binsCent[nCentBins] };
 906
 907     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
 908       fPtResolution[i] = new THnSparseD(Form("fPtResolution_%s", AliPID::ParticleShortName(i)),
 909                                         Form("Pt resolution for primaries, %s", AliPID::ParticleLatexName(i)),
 910                                         nBinsPtResolution, ptResolutionBins, ptResolutionXmin, ptResolutionXmax);
 911       SetUpPtResHist(fPtResolution[i], pTbinsRes, binsJetPt, binsCent);
 912       fPtResolutionContainer->Add(fPtResolution[i]);
 913     }
 914   }
 915
 916   if(fDebug > 2)
 917     printf("File: %s, Line: %d: UserCreateOutputObjects -> Posting output data\n", (char*)__FILE__, __LINE__);
 918
 919   PostData(1, fOutputContainer);
 920   PostData(2, fContainerEff);
 921   PostData(3, fPtResolutionContainer);
 922
 923   if(fDebug > 2)
 924     printf("File: %s, Line: %d: UserCreateOutputObjects -> Done\n", (char*)__FILE__, __LINE__);
 925 }
 926
 927
 928 //________________________________________________________________________
 929 void AliAnalysisTaskPID::UserExec(Option_t *)
 930 {
 931   // Main loop
 932   // Called for each event
 933
 934   if(fDebug > 1)
 935     printf("File: %s, Line: %d: UserExec\n", (char*)__FILE__, __LINE__);
 936
 937   Int_t run = InputEvent()->GetRunNumber();
 938
 939   if (run != fRun){
 940     // If systematics on eta is investigated, need to calculate the maxEtaVariationMap
 941     if ((TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon) ||
 942         (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon)) {
 943       if (!CalculateMaxEtaVariationMapFromPIDResponse())
 944         AliFatal("Systematics on eta correction requested, but failed to calculate max eta varation map!");
 945     }
 946   }
 947
 948   fRun = run;
 949
 950   // No processing of event, if input is fed in directly from another task
 951   if (fInputFromOtherTask)
 952     return;
 953
 954   if(fDebug > 1)
 955     printf("File: %s, Line: %d: UserExec -> Processing started\n", (char*)__FILE__, __LINE__);
 956
 957   fEvent = dynamic_cast<AliVEvent*>(InputEvent());
 958   if (!fEvent) {
 959     Printf("ERROR: fEvent not available");
 960     return;
 961   }
 962
 963   fMC = dynamic_cast<AliMCEvent*>(MCEvent());
 964
 965   if (!fPIDResponse || !fPIDcombined)
 966     return;
 967
 968   Double_t centralityPercentile = -1;
 969   if (fStoreCentralityPercentile)
 970     centralityPercentile = fEvent->GetCentrality()->GetCentralityPercentile(fCentralityEstimator.Data());
 971
 972   IncrementEventCounter(centralityPercentile, kTriggerSel);
 973
 974   // Check if vertex is ok, but don't apply cut on z position
 975   if (!GetVertexIsOk(fEvent, kFALSE))
 976     return;
 977
 978   fESD = dynamic_cast<AliESDEvent*>(fEvent);
 979   const AliVVertex* primaryVertex = fESD ? fESD->GetPrimaryVertexTracks() : fEvent->GetPrimaryVertex();
 980   if (!primaryVertex)
 981     return;
 982
 983   if(primaryVertex->GetNContributors() <= 0)
 984     return;
 985
 986   IncrementEventCounter(centralityPercentile, kTriggerSelAndVtxCut);
 987
 988   // Now check again, but also require z position to be in desired range
 989   if (!GetVertexIsOk(fEvent, kTRUE))
 990     return;
 991
 992   IncrementEventCounter(centralityPercentile, kTriggerSelAndVtxCutAndZvtxCut);
 993
 994   Double_t magField = fEvent->GetMagneticField();
 995
 996   if (fMC) {
 997     if (fDoPID || fDoEfficiency) {
 998       for (Int_t iPart = 0; iPart < fMC->GetNumberOfTracks(); iPart++) {
 999         AliMCParticle *mcPart  = dynamic_cast<AliMCParticle*>(fMC->GetTrack(iPart));
1000
1001         if (!mcPart)
1002             continue;
1003
1004         // Define clean MC sample with corresponding particle level track cuts:
1005         // - MC-track must be in desired eta range
1006         // - MC-track must be physical primary
1007         // - Species must be one of those in question (everything else goes to the overflow bin of mcID)
1008
1009         // Geometrie should be the same as on the reconstructed level -> By definition analysis within this eta interval
1010         if (!IsInAcceptedEtaRange(TMath::Abs(mcPart->Eta()))) continue;
1011
1012         Int_t mcID = PDGtoMCID(mcPart->PdgCode());
1013
1014         // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
1015         Double_t chargeMC = mcPart->Charge() / 3.;
1016
1017         if (TMath::Abs(chargeMC) < 0.01)
1018           continue; // Reject neutral particles (only relevant, if mcID is not used)
1019
1020         if (!fMC->IsPhysicalPrimary(iPart))
1021             continue;
1022
1023         if (fDoPID) {
1024           Double_t valuesGenYield[kGenYieldNumAxes] = { mcID, mcPart->Pt(), centralityPercentile, -1, -1, -1, -1 };
1025           valuesGenYield[GetIndexOfChargeAxisGenYield()] = chargeMC;
1026
1027           fhMCgeneratedYieldsPrimaries->Fill(valuesGenYield);
1028         }
1029
1030
1031         if (fDoEfficiency) {
1032           Double_t valueEff[kEffNumAxes] = { mcID, mcPart->Pt(), mcPart->Eta(), chargeMC, centralityPercentile,
1033                                             -1, -1, -1 };
1034
1035           fContainerEff->Fill(valueEff, kStepGenWithGenCuts);
1036         }
1037       }
1038     }
1039   }
1040
1041   // Track loop to fill a Train spectrum
1042   for (Int_t iTracks = 0; iTracks < fEvent->GetNumberOfTracks(); iTracks++) {
1043     AliVTrack* track = dynamic_cast<AliVTrack*>(fEvent->GetTrack(iTracks));
1044     if (!track) {
1045       Printf("ERROR: Could not retrieve track %d", iTracks);
1046       continue;
1047     }
1048
1049
1050     // Apply detector level track cuts
1051     const Bool_t tuneOnDataTPC = fPIDResponse->IsTunedOnData() &&
1052                                  ((fPIDResponse->GetTunedOnDataMask() & AliPIDResponse::kDetTPC) == AliPIDResponse::kDetTPC);
1053     Double_t dEdxTPC = tuneOnDataTPC ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
1054     if (dEdxTPC <= 0)
1055       continue;
1056
1057     if(fTrackFilter && !fTrackFilter->IsSelected(track))
1058       continue;
1059
1060     if (GetUseTPCCutMIGeo()) {
1061       if (!TPCCutMIGeo(track, fEvent))
1062         continue;
1063     }
1064     else if (GetUseTPCnclCut()) {
1065       if (!TPCnclCut(track))
1066         continue;
1067     }
1068
1069     if(fUsePhiCut) {
1070       if (!PhiPrimeCut(track, magField))
1071         continue; // reject track
1072     }
1073
1074     Int_t pdg =  0; // = 0 indicates data for the moment
1075     AliMCParticle* mcTrack = 0x0;
1076     Int_t mcID = AliPID::kUnknown;
1077     Int_t label = 0;
1078
1079     if (fMC) {
1080       label = track->GetLabel();
1081
1082       //if (label < 0)
1083       //  continue;
1084
1085       mcTrack = dynamic_cast<AliMCParticle*>(fMC->GetTrack(TMath::Abs(label)));
1086       if (!mcTrack) {
1087         Printf("ERROR: Could not retrieve mcTrack with label %d for track %d", label, iTracks);
1088         continue;
1089       }
1090
1091       pdg = mcTrack->PdgCode();
1092       mcID = PDGtoMCID(mcTrack->PdgCode());
1093
1094       if (fDoEfficiency) {
1095         // For efficiency: Reconstructed track has survived all cuts on the detector level (excluding acceptance)
1096         // and has an associated MC track which is a physical primary and was generated inside the acceptance
1097         if (fMC->IsPhysicalPrimary(TMath::Abs(label)) &&
1098             IsInAcceptedEtaRange(TMath::Abs(mcTrack->Eta()))) {
1099
1100           // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
1101           Double_t value[kEffNumAxes] = { mcID, mcTrack->Pt(), mcTrack->Eta(), mcTrack->Charge() / 3., centralityPercentile,
1102                                           -1, -1, -1 };
1103           fContainerEff->Fill(value, kStepRecWithGenCuts);
1104
1105           Double_t valueMeas[kEffNumAxes] = { mcID, track->Pt(), track->Eta(), track->Charge(), centralityPercentile,
1106                                               -1, -1, -1 };
1107           fContainerEff->Fill(valueMeas, kStepRecWithGenCutsMeasuredObs);
1108         }
1109       }
1110     }
1111
1112     // Only process tracks inside the desired eta window
1113     if (!IsInAcceptedEtaRange(TMath::Abs(track->Eta()))) continue;
1114
1115     if (fDoPID)
1116       ProcessTrack(track, pdg, centralityPercentile, -1); // No jet information in this case -> Set jet pT to -1
1117
1118     if (fDoPtResolution) {
1119       if (mcTrack && fMC->IsPhysicalPrimary(TMath::Abs(label))) {
1120         // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
1121         Double_t valuePtRes[kPtResNumAxes] = { -1, mcTrack->Pt(), track->Pt(), mcTrack->Charge() / 3., centralityPercentile };
1122         FillPtResolution(mcID, valuePtRes);
1123       }
1124     }
1125
1126     if (fDoEfficiency) {
1127       if (mcTrack) {
1128         Double_t valueRecAllCuts[kEffNumAxes] = { mcID, track->Pt(), track->Eta(), track->Charge(), centralityPercentile,
1129                                                   -1, -1, -1 };
1130         fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObs);
1131
1132         Double_t weight = IsSecondaryWithStrangeMotherMC(fMC, TMath::Abs(label)) ?
1133                                                                 GetMCStrangenessFactorCMS(fMC, mcTrack) : 1.0;
1134         fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObsStrangenessScaled, weight);
1135
1136         // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
1137         Double_t valueGenAllCuts[kEffNumAxes] = { mcID, mcTrack->Pt(), mcTrack->Eta(), mcTrack->Charge() / 3.,
1138                                                   centralityPercentile, -1, -1, -1 };
1139         if (fMC->IsPhysicalPrimary(TMath::Abs(label))) {
1140           fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObsPrimaries);
1141           fContainerEff->Fill(valueGenAllCuts, kStepRecWithRecCutsPrimaries);
1142         }
1143       }
1144     }
1145   } //track loop
1146
1147   if(fDebug > 2)
1148     printf("File: %s, Line: %d: UserExec -> Processing done\n", (char*)__FILE__, __LINE__);
1149
1150   PostOutputData();
1151
1152   if(fDebug > 2)
1153     printf("File: %s, Line: %d: UserExec -> Done\n", (char*)__FILE__, __LINE__);
1154 }
1155
1156 //________________________________________________________________________
1157 void AliAnalysisTaskPID::Terminate(const Option_t *)
1158 {
1159   // Draw result to the screen
1160   // Called once at the end of the query
1161 }
1162
1163
1164 //_____________________________________________________________________________
1165 void AliAnalysisTaskPID::CheckDoAnyStematicStudiesOnTheExpectedSignal()
1166 {
1167   // Check whether at least one scale factor indicates the ussage of systematic studies
1168   // and set internal flag accordingly.
1169
1170   fDoAnySystematicStudiesOnTheExpectedSignal = kFALSE;
1171
1172
1173   if ((TMath::Abs(fSystematicScalingSplinesBelowThreshold - 1.0) > fgkEpsilon) ||
1174       (TMath::Abs(fSystematicScalingSplinesAboveThreshold - 1.0) > fgkEpsilon)) {
1175     fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1176     return;
1177   }
1178
1179   if ((TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon) ||
1180       (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon)) {
1181     fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1182     return;
1183   }
1184
1185   if ((TMath::Abs(fSystematicScalingEtaSigmaParaBelowThreshold - 1.0) > fgkEpsilon) ||
1186       (TMath::Abs(fSystematicScalingEtaSigmaParaAboveThreshold - 1.0) > fgkEpsilon)) {
1187     fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1188     return;
1189   }
1190
1191   if (TMath::Abs(fSystematicScalingMultCorrection - 1.0) > fgkEpsilon) {
1192     fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1193     return;
1194   }
1195 }
1196
1197
1198 //_____________________________________________________________________________
1199 Int_t AliAnalysisTaskPID::PDGtoMCID(Int_t pdg)
1200 {
1201   // Returns the corresponding AliPID index to the given pdg code.
1202   // Returns AliPID::kUnkown if pdg belongs to a not considered species.
1203
1204   Int_t absPDGcode = TMath::Abs(pdg);
1205   if (absPDGcode == 211) {//Pion
1206     return AliPID::kPion;
1207   }
1208   else if (absPDGcode == 321) {//Kaon
1209     return AliPID::kKaon;
1210   }
1211   else if (absPDGcode == 2212) {//Proton
1212     return AliPID::kProton;
1213   }
1214   else if (absPDGcode == 11) {//Electron
1215     return AliPID::kElectron;
1216   }
1217   else if (absPDGcode == 13) {//Muon
1218     return AliPID::kMuon;
1219   }
1220
1221   return AliPID::kUnknown;
1222 }
1223
1224
1225 //_____________________________________________________________________________
1226 void AliAnalysisTaskPID::GetJetTrackObservables(Double_t trackPt, Double_t jetPt, Double_t& z, Double_t& xi)
1227 {
1228   // Uses trackPt and jetPt to obtain z and xi.
1229
1230   z = (jetPt > 0 && trackPt >= 0) ? (trackPt / jetPt) : -1;
1231   xi = (z > 0) ? TMath::Log(1. / z) : -1;
1232
1233   if(trackPt > (1. - 1e-06) * jetPt && trackPt < (1. + 1e-06) * jetPt) { // case z=1 : move entry to last histo bin <1
1234     z  = 1. - 1e-06;
1235     xi = 1e-06;
1236   }
1237 }
1238
1239
1240 //_____________________________________________________________________________
1241 void AliAnalysisTaskPID::CleanupParticleFractionHistos()
1242 {
1243   // Delete histos with particle fractions
1244
1245   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1246     delete fFractionHists[i];
1247     fFractionHists[i] = 0x0;
1248
1249     delete fFractionSysErrorHists[i];
1250     fFractionSysErrorHists[i] = 0x0;
1251   }
1252 }
1253
1254
1255 //_____________________________________________________________________________
1256 Double_t AliAnalysisTaskPID::ConvolutedGaus(const Double_t* xx, const Double_t* par) const
1257 {
1258   // Convolutes gauss with an exponential tail which describes dEdx-response better than pure gaussian
1259
1260   const Double_t mean = par[0];
1261   const Double_t sigma = par[1];
1262   const Double_t lambda = par[2];
1263
1264   if(fDebug > 5)
1265     printf("File: %s, Line: %d: ConvolutedGaus: mean %e, sigma %e, lambda %e\n", (char*)__FILE__, __LINE__, mean, sigma, lambda);
1266
1267   return lambda/sigma*TMath::Exp(-lambda/sigma*(xx[0]-mean)+lambda*lambda*0.5)*0.5*TMath::Erfc((-xx[0]+mean+sigma*lambda)/sigma*fgkOneOverSqrt2);
1268 }
1269
1270
1271 //_____________________________________________________________________________
1272 inline Double_t AliAnalysisTaskPID::FastGaus(Double_t x, Double_t mean, Double_t sigma) const
1273 {
1274   // Calculate an unnormalised gaussian function with mean and sigma.
1275
1276   if (sigma < fgkEpsilon)
1277     return 1.e30;
1278
1279   const Double_t arg = (x - mean) / sigma;
1280   return exp(-0.5 * arg * arg);
1281 }
1282
1283
1284 //_____________________________________________________________________________
1285 inline Double_t AliAnalysisTaskPID::FastNormalisedGaus(Double_t x, Double_t mean, Double_t sigma) const
1286 {
1287   // Calculate a normalised (divided by sqrt(2*Pi)*sigma) gaussian function with mean and sigma.
1288
1289   if (sigma < fgkEpsilon)
1290     return 1.e30;
1291
1292   const Double_t arg = (x - mean) / sigma;
1293   const Double_t res = exp(-0.5 * arg * arg);
1294   return res / (2.50662827463100024 * sigma); //sqrt(2*Pi)=2.50662827463100024
1295 }
1296
1297
1298 //_____________________________________________________________________________
1299 Int_t AliAnalysisTaskPID::FindBinWithinRange(TAxis* axis, Double_t value) const
1300 {
1301   // Find the corresponding bin of the axis. Values outside the range (also under and overflow) will be set to the first/last
1302   // available bin
1303
1304   Int_t bin = axis->FindFixBin(value);
1305
1306   if (bin <= 0)
1307     bin = 1;
1308   if (bin > axis->GetNbins())
1309     bin = axis->GetNbins();
1310
1311   return bin;
1312 }
1313
1314
1315 //_____________________________________________________________________________
1316 Int_t AliAnalysisTaskPID::FindFirstBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yBin,
1317                                                       Int_t zBin) const
1318 {
1319   // Kind of projects a TH3 to 1 bin combination in y and z
1320   // and looks for the first x bin above a threshold for this projection.
1321   // If no such bin is found, -1 is returned.
1322
1323   if (!hist)
1324     return -1;
1325
1326   Int_t nBinsX = hist->GetNbinsX();
1327   for (Int_t xBin = 1; xBin <= nBinsX; xBin++) {
1328     if (hist->GetBinContent(xBin, yBin, zBin) > threshold)
1329       return xBin;
1330   }
1331
1332   return -1;
1333 }
1334
1335
1336 //_____________________________________________________________________________
1337 Int_t AliAnalysisTaskPID::FindLastBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yBin,
1338                                                      Int_t zBin) const
1339 {
1340   // Kind of projects a TH3 to 1 bin combination in y and z
1341   // and looks for the last x bin above a threshold for this projection.
1342   // If no such bin is found, -1 is returned.
1343
1344   if (!hist)
1345     return -1;
1346
1347   Int_t nBinsX = hist->GetNbinsX();
1348   for (Int_t xBin = nBinsX; xBin >= 1; xBin--) {
1349     if (hist->GetBinContent(xBin, yBin, zBin) > threshold)
1350       return xBin;
1351   }
1352
1353   return -1;
1354 }
1355
1356
1357 //_____________________________________________________________________________
1358 Bool_t AliAnalysisTaskPID::GetParticleFraction(Double_t trackPt, Double_t jetPt, Double_t centralityPercentile,
1359                                                AliPID::EParticleType species,
1360                                                Double_t& fraction, Double_t& fractionErrorStat, Double_t& fractionErrorSys) const
1361 {
1362   // Computes the particle fraction for the corresponding species for the given trackPt, jetPt and centrality.
1363   // Use jetPt = -1 for inclusive spectra and centralityPercentile = -1 for pp.
1364   // On success (return value kTRUE), fraction contains the particle fraction, fractionErrorStat(Sys) the sigma of its
1365   // statistical (systematic) error
1366
1367   fraction = -999.;
1368   fractionErrorStat = 999.;
1369   fractionErrorSys = 999.;
1370
1371   if (species > AliPID::kProton || species < AliPID::kElectron) {
1372     AliError(Form("Only fractions for species index %d to %d availabe, but not for the requested one: %d", 0, AliPID::kProton, species));
1373     return kFALSE;
1374   }
1375
1376   if (!fFractionHists[species]) {
1377     AliError(Form("Histo with particle fractions for species %d not loaded!", species));
1378
1379     return kFALSE;
1380   }
1381
1382   Int_t jetPtBin = FindBinWithinRange(fFractionHists[species]->GetYaxis(), jetPt);
1383   Int_t centBin  = FindBinWithinRange(fFractionHists[species]->GetZaxis(), centralityPercentile);
1384
1385   // The following interpolation takes the bin content of the first/last available bin,
1386   // if requested point lies beyond bin center of first/last bin.
1387   // The interpolation is only done for the x-axis (track pT), i.e. jetPtBin and centBin are fix,
1388   // because the analysis will anyhow be run in bins of jetPt and centrality and
1389   // it is not desired to correlate different jetPt bins via interpolation.
1390
1391   // The same procedure is used for the error of the fraction
1392   TAxis* xAxis = fFractionHists[species]->GetXaxis();
1393
1394   // No interpolation to values beyond the centers of the first/last bins (we don't know exactly where the spectra start or stop,
1395   // thus, search for the first and last bin above 0.0 to constrain the range
1396   Int_t firstBin = TMath::Max(1, FindFirstBinAboveIn3dSubset(fFractionHists[species], 0.0, jetPtBin, centBin));
1397   Int_t lastBin  = TMath::Min(fFractionHists[species]->GetNbinsX(),
1398                               FindLastBinAboveIn3dSubset(fFractionHists[species], 0.0, jetPtBin, centBin));
1399
1400   if (trackPt <= xAxis->GetBinCenter(firstBin)) {
1401     fraction = fFractionHists[species]->GetBinContent(firstBin, jetPtBin, centBin);
1402     fractionErrorStat = fFractionHists[species]->GetBinError(firstBin, jetPtBin, centBin);
1403     fractionErrorSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(firstBin, jetPtBin, centBin) : 0.;
1404   }
1405   else if (trackPt >= xAxis->GetBinCenter(lastBin)) {
1406     fraction = fFractionHists[species]->GetBinContent(lastBin, jetPtBin, centBin);
1407     fractionErrorStat = fFractionHists[species]->GetBinError(lastBin, jetPtBin, centBin);
1408     fractionErrorSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(lastBin, jetPtBin, centBin) : 0.;
1409   }
1410   else {
1411     Double_t x0 = 0., x1 = 0., y0 = 0., y1 = 0.;
1412     Double_t y0errStat = 0., y1errStat = 0., y0errSys = 0., y1errSys = 0.;
1413     Int_t trackPtBin = xAxis->FindBin(trackPt);
1414
1415     // Linear interpolation between nearest neighbours in trackPt
1416     if (trackPt <= xAxis->GetBinCenter(trackPtBin)) {
1417         y0 = fFractionHists[species]->GetBinContent(trackPtBin - 1, jetPtBin, centBin);
1418         y0errStat = fFractionHists[species]->GetBinError(trackPtBin - 1, jetPtBin, centBin);
1419         y0errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin - 1, jetPtBin, centBin)
1420                                                    : 0.;
1421         x0 = xAxis->GetBinCenter(trackPtBin - 1);
1422         y1 = fFractionHists[species]->GetBinContent(trackPtBin, jetPtBin, centBin);
1423         y1errStat = fFractionHists[species]->GetBinError(trackPtBin, jetPtBin, centBin);
1424         y1errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin, jetPtBin, centBin)
1425                                                    : 0.;
1426         x1 = xAxis->GetBinCenter(trackPtBin);
1427     }
1428     else {
1429         y0 = fFractionHists[species]->GetBinContent(trackPtBin, jetPtBin, centBin);
1430         y0errStat = fFractionHists[species]->GetBinError(trackPtBin, jetPtBin, centBin);
1431         y0errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin, jetPtBin, centBin)
1432                                                    : 0.;
1433         x0 = xAxis->GetBinCenter(trackPtBin);
1434         y1 = fFractionHists[species]->GetBinContent(trackPtBin + 1, jetPtBin, centBin);
1435         y1errStat = fFractionHists[species]->GetBinError(trackPtBin + 1, jetPtBin, centBin);
1436         y1errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin + 1, jetPtBin, centBin)
1437                                                    : 0.;
1438         x1 = xAxis->GetBinCenter(trackPtBin + 1);
1439     }
1440
1441     // Per construction: x0 < trackPt < x1
1442     fraction = y0 + (trackPt - x0) * ((y1 - y0) / (x1 - x0));
1443     fractionErrorStat = y0errStat + (trackPt - x0) * ((y1errStat - y0errStat) / (x1 - x0));
1444     fractionErrorSys = fFractionSysErrorHists[species] ? (y0errSys + (trackPt - x0) * ((y1errSys - y0errSys) / (x1 - x0))) : 0.;
1445   }
1446
1447   return kTRUE;
1448 }
1449
1450
1451 //_____________________________________________________________________________
1452 Bool_t AliAnalysisTaskPID::GetParticleFractions(Double_t trackPt, Double_t jetPt, Double_t centralityPercentile,
1453                                                 Double_t* prob, Int_t smearSpeciesByError,
1454                                                 Int_t takeIntoAccountSpeciesSysError, Bool_t uniformSystematicError) const
1455 {
1456   // Fills the particle fractions for the given trackPt, jetPt and centrality into "prob".
1457   // Use jetPt = -1 for inclusive spectra and centralityPercentile = -1 for pp.
1458   // If smearSpeciesByError is >= 0 && < AliPID::kSPECIES, the returned fractions will be a random number distributed
1459   // with a gauss with mean being the corresponding particle fraction and sigma it's error for the considered species
1460   // "smearSpeciesByError".
1461   // Note that in this case the fractions for all species will NOT sum up to 1!
1462   // Thus, all other species fractions will be re-scaled weighted with their corresponding statistical error.
1463   // A similar procedure is used for "takeIntoAccountSpeciesSysError":  The systematic error of the corresponding species
1464   // is used to generate a random number with uniform distribution in [mean - sysError, mean + sysError] for the new mean
1465   // (in cace of uniformSystematicError = kTRUE, otherwise it will be a gaus(mean, sysError)),
1466   // then the other species will be re-scaled according to their systematic errors.
1467   // First, the systematic error uncertainty procedure will be performed (that is including re-scaling), then the statistical
1468   // uncertainty procedure.
1469   // On success, kTRUE is returned.
1470
1471   if (!prob || smearSpeciesByError >= AliPID::kSPECIES || takeIntoAccountSpeciesSysError >= AliPID::kSPECIES)
1472     return kFALSE;
1473
1474   Double_t probTemp[AliPID::kSPECIES];
1475   Double_t probErrorStat[AliPID::kSPECIES];
1476   Double_t probErrorSys[AliPID::kSPECIES];
1477
1478   Bool_t success = kTRUE;
1479   success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kElectron,
1480                                            probTemp[AliPID::kElectron], probErrorStat[AliPID::kElectron],
1481                                            probErrorSys[AliPID::kElectron]);
1482   success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kMuon,
1483                                            probTemp[AliPID::kMuon], probErrorStat[AliPID::kMuon], probErrorSys[AliPID::kMuon]);
1484   success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kPion,
1485                                            probTemp[AliPID::kPion], probErrorStat[AliPID::kPion], probErrorSys[AliPID::kPion]);
1486   success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kKaon,
1487                                            probTemp[AliPID::kKaon], probErrorStat[AliPID::kKaon], probErrorSys[AliPID::kKaon]);
1488   success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kProton,
1489                                            probTemp[AliPID::kProton], probErrorStat[AliPID::kProton], probErrorSys[AliPID::kProton]);
1490
1491   if (!success)
1492     return kFALSE;
1493
1494   // If desired, take into account the systematic error of the corresponding species and re-generate probTemp accordingly
1495   if (takeIntoAccountSpeciesSysError >= 0) {
1496     // Generate random fraction of the considered species "smearSpeciesByError" according to mean and sigma
1497     Double_t generatedFraction = uniformSystematicError
1498                                    ? fRandom->Rndm() * 2. * probErrorSys[takeIntoAccountSpeciesSysError]
1499                                      - probErrorSys[takeIntoAccountSpeciesSysError]
1500                                      + probTemp[takeIntoAccountSpeciesSysError]
1501                                    : fRandom->Gaus(probTemp[takeIntoAccountSpeciesSysError],
1502                                                    probErrorSys[takeIntoAccountSpeciesSysError]);
1503
1504     // Catch cases with invalid fraction (can happen for large errors), i.e. fraction < 0 or > 1
1505     if (generatedFraction < 0.)
1506       generatedFraction = 0.;
1507     else if (generatedFraction > 1.)
1508       generatedFraction = 1.;
1509
1510     // Calculate difference from original fraction (original fractions sum up to 1!)
1511     Double_t deltaFraction = generatedFraction - probTemp[takeIntoAccountSpeciesSysError];
1512
1513     // Fractions must (including errors) lie inside [0,1] -> Adapt weights accordingly by setting the errors
1514     if (deltaFraction > 0) {
1515       // Some part will be SUBTRACTED from the other fractions
1516       for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1517         if (probTemp[i] - probErrorSys[i] < 0)
1518           probErrorSys[i] = probTemp[i];
1519       }
1520     }
1521     else {
1522       // Some part will be ADDED to the other fractions
1523       for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1524         if (probTemp[i] + probErrorSys[i] > 1)
1525           probErrorSys[i] = 1. - probTemp[i];
1526       }
1527     }
1528
1529     // Compute summed weight of all fractions except for the considered one
1530     Double_t summedWeight = 0.;
1531     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1532       if (i != takeIntoAccountSpeciesSysError)
1533         summedWeight += probErrorSys[i];
1534     }
1535
1536     // Compute the weight for the other species
1537     /*
1538     if (summedWeight <= 1e-13) {
1539       // If this happens for some reason (it should not!), just assume flat weight
1540       printf("Error: summedWeight (sys error) ~ 0 for trackPt %f, jetPt %f, centralityPercentile %f. Setting flat weight!\n",
1541              trackPt, jetPt, centralityPercentile);
1542     }*/
1543
1544     Double_t weight[AliPID::kSPECIES];
1545     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1546       if (i != takeIntoAccountSpeciesSysError) {
1547         if (summedWeight > 1e-13)
1548           weight[i] = probErrorSys[i] / summedWeight;
1549         else
1550           weight[i] = probErrorSys[i] / (AliPID::kSPECIES - 1);
1551       }
1552     }
1553
1554     // For the final generated fractions, set the generated value for the considered species
1555     // and the generated value minus delta times statistical weight
1556     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1557       if (i != takeIntoAccountSpeciesSysError)
1558         probTemp[i] = probTemp[i] - weight[i] * deltaFraction;
1559       else
1560         probTemp[i] = generatedFraction;
1561     }
1562   }
1563
1564   // Using the values of probTemp (either the original ones or those after taking into account the systematic error),
1565   // calculate the final fractions - if the statistical error is to be taken into account, smear the corresponding
1566   // fraction. If not, just write probTemp to the final result array.
1567   if (smearSpeciesByError >= 0) {
1568     // Generate random fraction of the considered species "smearSpeciesByError" according to mean and sigma
1569     Double_t generatedFraction = fRandom->Gaus(probTemp[smearSpeciesByError], probErrorStat[smearSpeciesByError]);
1570
1571     // Catch cases with invalid fraction (can happen for large errors), i.e. fraction < 0 or > 1
1572     if (generatedFraction < 0.)
1573       generatedFraction = 0.;
1574     else if (generatedFraction > 1.)
1575       generatedFraction = 1.;
1576
1577     // Calculate difference from original fraction (original fractions sum up to 1!)
1578     Double_t deltaFraction = generatedFraction - probTemp[smearSpeciesByError];
1579
1580     // Fractions must (including errors) lie inside [0,1] -> Adapt weights accordingly by setting the errors
1581     if (deltaFraction > 0) {
1582       // Some part will be SUBTRACTED from the other fractions
1583       for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1584         if (probTemp[i] - probErrorStat[i] < 0)
1585           probErrorStat[i] = probTemp[i];
1586       }
1587     }
1588     else {
1589       // Some part will be ADDED to the other fractions
1590       for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1591         if (probTemp[i] + probErrorStat[i] > 1)
1592           probErrorStat[i] = 1. - probTemp[i];
1593       }
1594     }
1595
1596     // Compute summed weight of all fractions except for the considered one
1597     Double_t summedWeight = 0.;
1598     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1599       if (i != smearSpeciesByError)
1600         summedWeight += probErrorStat[i];
1601     }
1602
1603     // Compute the weight for the other species
1604     /*
1605     if (summedWeight <= 1e-13) {
1606       // If this happens for some reason (it should not!), just assume flat weight
1607       printf("Error: summedWeight (stat error) ~ 0 for trackPt %f, jetPt %f, centralityPercentile %f. Setting flat weight!\n",
1608              trackPt, jetPt, centralityPercentile);
1609     }*/
1610
1611     Double_t weight[AliPID::kSPECIES];
1612     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1613       if (i != smearSpeciesByError) {
1614         if (summedWeight > 1e-13)
1615           weight[i] = probErrorStat[i] / summedWeight;
1616         else
1617           weight[i] = probErrorStat[i] / (AliPID::kSPECIES - 1);
1618       }
1619     }
1620
1621     // For the final generated fractions, set the generated value for the considered species
1622     // and the generated value minus delta times statistical weight
1623     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1624       if (i != smearSpeciesByError)
1625         prob[i] = probTemp[i] - weight[i] * deltaFraction;
1626       else
1627         prob[i] = generatedFraction;
1628     }
1629   }
1630   else {
1631     // Just take the generated values
1632     for (Int_t i = 0; i < AliPID::kSPECIES; i++)
1633       prob[i] = probTemp[i];
1634   }
1635
1636
1637   // Should already be normalised, but make sure that it really is:
1638   Double_t probSum = 0.;
1639   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1640     probSum += prob[i];
1641   }
1642
1643   if (probSum <= 0)
1644     return kFALSE;
1645
1646   if (TMath::Abs(probSum - 1.0) > 1e-4) {
1647     printf("Warning: Re-normalising sum of fractions: Sum is %e\n", probSum);
1648     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1649       prob[i] /= probSum;
1650     }
1651   }
1652
1653   return kTRUE;
1654 }
1655
1656
1657 //_____________________________________________________________________________
1658 const TH3D* AliAnalysisTaskPID::GetParticleFractionHisto(Int_t species, Bool_t sysError) const
1659 {
1660   if (species < AliPID::kElectron || species > AliPID::kProton)
1661     return 0x0;
1662
1663   return sysError ? fFractionSysErrorHists[species] : fFractionHists[species];
1664 }
1665
1666
1667 //_____________________________________________________________________________
1668 Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(Int_t motherPDG, Double_t motherGenPt)
1669 {
1670   // Strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
1671   // -> Based on function in PWGJE/AliAnalysisTaskFragmentationFunction, which uses
1672   // the following data from CMS pp @ 7 TeV inclusive (JHEP 05 (2011) 064)
1673
1674   Double_t fac = 1;
1675
1676   const Int_t absMotherPDG = TMath::Abs(motherPDG);
1677
1678   if (absMotherPDG == 310 || absMotherPDG == 321) { // K0s / K+ / K-
1679     if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.768049;
1680     else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.732933;
1681     else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.650298;
1682     else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.571332;
1683     else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.518734;
1684     else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.492543;
1685     else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.482704;
1686     else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.488056;
1687     else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.488861;
1688     else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.492862;
1689     else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.504332;
1690     else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.501858;
1691     else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.512970;
1692     else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.524131;
1693     else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.539130;
1694     else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.554101;
1695     else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.560348;
1696     else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.568869;
1697     else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.583310;
1698     else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.604818;
1699     else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.632630;
1700     else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.710070;
1701     else if(6.00 <= motherGenPt && motherGenPt < 8.00) fac = 0.736365;
1702     else if(8.00 <= motherGenPt && motherGenPt < 10.00) fac = 0.835865;
1703   }
1704
1705   if (absMotherPDG == 3122) { // Lambda
1706     if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.645162;
1707     else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.627431;
1708     else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.457136;
1709     else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.384369;
1710     else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.330597;
1711     else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.309571;
1712     else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.293620;
1713     else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.283709;
1714     else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.282047;
1715     else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.277261;
1716     else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.275772;
1717     else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.280726;
1718     else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.288540;
1719     else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.288315;
1720     else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.296619;
1721     else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.302993;
1722     else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.338121;
1723     else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.349800;
1724     else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.356802;
1725     else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.391202;
1726     else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.422573;
1727     else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.573815;
1728     else if(6.00 <= motherGenPt && motherGenPt < 8.00) fac = 0.786984;
1729     else if(8.00 <= motherGenPt && motherGenPt < 10.00) fac = 1.020021;
1730   }
1731
1732   if (absMotherPDG == 3312 || absMotherPDG == 3322) { // xi
1733     if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.666620;
1734     else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.575908;
1735     else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.433198;
1736     else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.340901;
1737     else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.290896;
1738     else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.236074;
1739     else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.218681;
1740     else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.207763;
1741     else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.222848;
1742     else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.208806;
1743     else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.197275;
1744     else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.183645;
1745     else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.188788;
1746     else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.188282;
1747     else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.207442;
1748     else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.240388;
1749     else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.241916;
1750     else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.208276;
1751     else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.234550;
1752     else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.251689;
1753     else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.310204;
1754     else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.343492;
1755   }
1756
1757   const Double_t weight = 1. / fac;
1758
1759   return weight;
1760 }
1761
1762
1763 //_____________________________________________________________________________
1764 Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(AliMCEvent* mcEvent, AliMCParticle* daughter)
1765 {
1766   // Strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
1767   // -> Based on function in PWGJE/AliAnalysisTaskFragmentationFunction
1768
1769   if (!mcEvent)
1770     return 1.;
1771
1772   AliMCParticle* currentMother   = daughter;
1773   AliMCParticle* currentDaughter = daughter;
1774
1775
1776   // find first primary mother K0s, Lambda or Xi
1777   while(1) {
1778     Int_t daughterPDG = currentDaughter->PdgCode();
1779
1780     Int_t motherLabel = currentDaughter->GetMother();
1781     if(motherLabel >= mcEvent->GetNumberOfTracks()){ // protection
1782       currentMother = currentDaughter;
1783       break;
1784     }
1785
1786     currentMother = (AliMCParticle*)mcEvent->GetTrack(motherLabel);
1787
1788     if (!currentMother) {
1789       currentMother = currentDaughter;
1790       break;
1791     }
1792
1793     Int_t motherPDG = currentMother->PdgCode();
1794
1795     // phys. primary found ?
1796     if (mcEvent->IsPhysicalPrimary(motherLabel))
1797       break;
1798
1799     if (TMath::Abs(daughterPDG) == 321) {
1800       // K+/K- e.g. from phi (ref data not feeddown corrected)
1801       currentMother = currentDaughter;
1802       break;
1803     }
1804     if (TMath::Abs(motherPDG) == 310) {
1805       // K0s e.g. from phi (ref data not feeddown corrected)
1806       break;
1807     }
1808     if (TMath::Abs(motherPDG) == 3212 && TMath::Abs(daughterPDG) == 3122) {
1809       // Mother Sigma0, daughter Lambda (this case not included in feeddown corr.)
1810       currentMother = currentDaughter;
1811       break;
1812     }
1813
1814     currentDaughter = currentMother;
1815   }
1816
1817
1818   Int_t motherPDG = currentMother->PdgCode();
1819   Double_t motherGenPt = currentMother->Pt();
1820
1821   return GetMCStrangenessFactorCMS(motherPDG, motherGenPt);
1822 }
1823
1824
1825 // _________________________________________________________________________________
1826 AliAnalysisTaskPID::TOFpidInfo AliAnalysisTaskPID::GetTOFType(const AliVTrack* track, Int_t tofMode) const
1827 {
1828   // Get the (locally defined) particle type judged by TOF
1829
1830   if (!fPIDResponse) {
1831     Printf("ERROR: fEvent not available -> Cannot determine TOF type!");
1832     return kNoTOFinfo;
1833   }
1834
1835   // Check kTOFout, kTIME, mismatch
1836   const AliPIDResponse::EDetPidStatus tofStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, track);
1837   if (tofStatus != AliPIDResponse::kDetPidOk)
1838     return kNoTOFinfo;
1839
1840   Double_t nsigma[kNumTOFspecies + 1] = { -999., -999., -999., -999. };
1841   nsigma[kTOFpion]   = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion);
1842   nsigma[kTOFkaon]   = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon);
1843   nsigma[kTOFproton] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton);
1844
1845   Double_t inclusion = -999;
1846   Double_t exclusion = -999;
1847
1848   if (tofMode == 0) {
1849     inclusion = 3.;
1850     exclusion = 2.5;
1851   }
1852   else if (tofMode == 1) { // default
1853     inclusion = 3.;
1854     exclusion = 3.;
1855   }
1856   else if (tofMode == 2) {
1857     inclusion = 3.;
1858     exclusion = 3.5;
1859   }
1860   else {
1861     Printf("ERROR: Bad TOF mode: %d!", tofMode);
1862     return kNoTOFinfo;
1863   }
1864
1865   // Do not cut on nSigma electron because this would also remove pions in a large pT range.
1866   // The electron contamination of the pion sample can be treated by TPC dEdx fits afterwards.
1867   if (TMath::Abs(nsigma[kTOFpion]) < inclusion && TMath::Abs(nsigma[kTOFkaon]) > exclusion && TMath::Abs(nsigma[kTOFproton]) > exclusion)
1868     return kTOFpion;
1869   if (TMath::Abs(nsigma[kTOFpion]) > exclusion && TMath::Abs(nsigma[kTOFkaon]) < inclusion && TMath::Abs(nsigma[kTOFproton]) > exclusion)
1870     return kTOFkaon;
1871   if (TMath::Abs(nsigma[kTOFpion]) > exclusion && TMath::Abs(nsigma[kTOFkaon]) > exclusion && TMath::Abs(nsigma[kTOFproton]) < inclusion)
1872     return kTOFproton;
1873
1874   // There are no TOF electrons selected because the purity is rather bad, even for small momenta
1875   // (also a small mismatch probability significantly affects electrons because their fraction
1876   // is small). There is no need for TOF electrons anyway, since the dEdx distribution of kaons and
1877   // protons in a given pT bin (also at the dEdx crossings) is very different from that of electrons.
1878   // This is due to the steeply falling dEdx of p and K with momentum, whereas the electron dEdx stays
1879   // rather constant.
1880   // As a result, the TPC fit yields a more accurate electron fraction than the TOF selection can do.
1881
1882   return kNoTOFpid;
1883 }
1884
1885
1886 // _________________________________________________________________________________
1887 Bool_t AliAnalysisTaskPID::IsSecondaryWithStrangeMotherMC(AliMCEvent* mcEvent, Int_t partLabel)
1888 {
1889   // Check whether particle is a secondary with strange mother, i.e. returns kTRUE if a strange mother is found
1890   // and the particle is NOT a physical primary. In all other cases kFALSE is returned
1891
1892   if (!mcEvent || partLabel < 0)
1893     return kFALSE;
1894
1895   AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(partLabel);
1896
1897   if (!part)
1898     return kFALSE;
1899
1900   if (mcEvent->IsPhysicalPrimary(partLabel))
1901     return kFALSE;
1902
1903   Int_t iMother = part->GetMother();
1904   if (iMother < 0)
1905     return kFALSE;
1906
1907
1908   AliMCParticle* partM = (AliMCParticle*)mcEvent->GetTrack(iMother);
1909   if (!partM)
1910    return kFALSE;
1911
1912   Int_t codeM = TMath::Abs(partM->PdgCode());
1913   Int_t mfl = Int_t(codeM / TMath::Power(10, Int_t(TMath::Log10(codeM))));
1914   if (mfl == 3 && codeM != 3) // codeM = 3 is for s quark
1915     return kTRUE;
1916
1917   return kFALSE;
1918 }
1919
1920
1921 //_____________________________________________________________________________
1922 Bool_t AliAnalysisTaskPID::SetParticleFractionHisto(const TH3D* hist, Int_t species, Bool_t sysError)
1923 {
1924   // Store a clone of hist (containing the particle fractions of the corresponding species with statistical error (sysError = kFALSE)
1925   // or systematic error (sysError = kTRUE), respectively), internally
1926
1927   if (species < AliPID::kElectron || species > AliPID::kProton) {
1928     AliError(Form("Only fractions for species index %d to %d can be set, but not for the requested one: %d", 0,
1929                   AliPID::kProton, species));
1930     return kFALSE;
1931   }
1932
1933   if (sysError) {
1934     delete fFractionSysErrorHists[species];
1935
1936     fFractionSysErrorHists[species] = new TH3D(*hist);
1937   }
1938   else {
1939     delete fFractionHists[species];
1940
1941     fFractionHists[species] = new TH3D(*hist);
1942   }
1943
1944   return kTRUE;
1945 }
1946
1947
1948 //_____________________________________________________________________________
1949 Bool_t AliAnalysisTaskPID::SetParticleFractionHistosFromFile(const TString filePathName, Bool_t sysError)
1950 {
1951   // Loads particle fractions for all species from the desired file and returns kTRUE on success.
1952   // The maps are assumed to be of Type TH3D, to sit in the main directory and to have names
1953   // Form("hFraction_%e", AliPID::ParticleName(i)) for sysError = kFALSE and
1954   // Form("hFractionSysError_%e", AliPID::ParticleName(i)) for sysError = kTRUE.
1955
1956   TFile* f = TFile::Open(filePathName.Data());
1957   if (!f)  {
1958     std::cout << "Failed to open file with particle fractions \"" << filePathName.Data() << "\"!" << std::endl;
1959     return kFALSE;
1960   }
1961
1962   TH3D* hist = 0x0;
1963   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1964     TString histName = Form("hFraction%s_%s", sysError ? "SysError" : "", AliPID::ParticleName(i));
1965     hist = dynamic_cast<TH3D*>(f->Get(histName.Data()));
1966     if (!hist) {
1967       std::cout << "Failed to load particle fractions for " << histName.Data() << "!";
1968       std::cout << std::endl << "Cleaning up particle fraction histos!" << std::endl;
1969       CleanupParticleFractionHistos();
1970       return kFALSE;
1971     }
1972
1973     if (!SetParticleFractionHisto(hist, i, sysError)) {
1974       std::cout << "Failed to load particle fractions for " << histName.Data() << "!";
1975       std::cout << std::endl << "Cleaning up particle fraction histos!" << std::endl;
1976       CleanupParticleFractionHistos();
1977       return kFALSE;
1978     }
1979   }
1980
1981   delete hist;
1982
1983   return kTRUE;
1984
1985 }
1986
1987
1988 //_____________________________________________________________________________
1989 Double_t AliAnalysisTaskPID::GetMaxEtaVariation(Double_t dEdxSplines)
1990 {
1991   // Returns the maximum eta variation (i.e. deviation of eta correction factor from unity) for the
1992   // given (spline) dEdx
1993
1994   if (dEdxSplines < 1. || !fhMaxEtaVariation) {
1995     Printf("Error GetMaxEtaVariation: No map or invalid dEdxSplines (%f)!", dEdxSplines);
1996     return 999.;
1997   }
1998
1999   Int_t bin = fhMaxEtaVariation->GetXaxis()->FindFixBin(1. / dEdxSplines);
2000
2001   if (bin == 0)
2002     bin = 1;
2003   if (bin > fhMaxEtaVariation->GetXaxis()->GetNbins())
2004     bin = fhMaxEtaVariation->GetXaxis()->GetNbins();
2005
2006   return fhMaxEtaVariation->GetBinContent(bin);
2007 }
2008
2009
2010 //_____________________________________________________________________________
2011 Int_t AliAnalysisTaskPID::GetRandomParticleTypeAccordingToParticleFractions(Double_t trackPt, Double_t jetPt,
2012                                                                             Double_t centralityPercentile,
2013                                                                             Bool_t smearByError,
2014                                                                             Bool_t takeIntoAccountSysError) const
2015 {
2016   // Uses the stored histograms with the particle fractions to generate a random particle type according to these fractions.
2017   // In case of problems (e.g. histo missing), AliPID::kUnknown is returned.
2018   // If smearByError is kTRUE, the used fractions will be random numbers distributed with a gauss with mean
2019   // being the corresponding particle fraction and sigma it's error.
2020   // Note that in this case only the fraction of a random species is varied in this way. The other fractions
2021   // will be re-normalised according their statistical errors.
2022   // The same holds for the systematic error of species "takeIntoAccountSpeciesSysError", but the random number will be
2023   // uniformly distributed within [mean - sys, mean + sys] and the re-normalisation will be weighted with the systematic errors.
2024   // Note that the fractions will be calculated first with only the systematic error taken into account (if desired), including
2025   // re-normalisation. Then, the resulting fractions will be used to calculate the final fractions - either with statistical error
2026   // or without. The species, for which the error will be used for smearing, is the same for sys and stat error.
2027
2028   Double_t prob[AliPID::kSPECIES];
2029   Int_t randomSpecies = (smearByError || takeIntoAccountSysError) ? (Int_t)(fRandom->Rndm() * AliPID::kSPECIES) : -1;
2030   Bool_t success = GetParticleFractions(trackPt, jetPt, centralityPercentile, prob, randomSpecies, randomSpecies);
2031
2032   if (!success)
2033     return AliPID::kUnknown;
2034
2035   Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
2036
2037   if (rnd <= prob[AliPID::kPion])
2038     return AliPID::kPion;
2039   else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon])
2040     return AliPID::kKaon;
2041   else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon] + prob[AliPID::kProton])
2042     return AliPID::kProton;
2043   else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon] + prob[AliPID::kProton] + prob[AliPID::kElectron])
2044     return AliPID::kElectron;
2045
2046   return AliPID::kMuon;  //else it must be a muon (only species left)
2047 }
2048
2049
2050 //_____________________________________________________________________________
2051 AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::GenerateDetectorResponse(AliAnalysisTaskPID::ErrorCode errCode,
2052                                                                            Double_t mean, Double_t sigma,
2053                                                                            Double_t* responses, Int_t nResponses,
2054                                                                            Bool_t usePureGaus)
2055 {
2056   // Generate detector response. If a previous generation was not successful or there is something wrong with this signal generation,
2057   // the function will return kFALSE
2058   if (!responses)
2059     return kError;
2060
2061   // Reset response array
2062   for (Int_t i = 0; i < nResponses; i++)
2063     responses[i] = -999;
2064
2065   if (errCode == kError)
2066     return kError;
2067
2068   ErrorCode ownErrCode = kNoErrors;
2069
2070   if (fUseConvolutedGaus && !usePureGaus) {
2071     // In case of convoluted gauss, calculate the probability density only once to save a lot of time!
2072
2073     TH1* hProbDensity = 0x0;
2074     ownErrCode = SetParamsForConvolutedGaus(mean, sigma);
2075     if (ownErrCode == kError)
2076       return kError;
2077
2078     hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
2079
2080     for (Int_t i = 0; i < nResponses; i++) {
2081       responses[i] = hProbDensity->GetRandom();
2082       //responses[i] fConvolutedGausDeltaPrime->GetRandom(); // MUCH slower than using the binned version via the histogram
2083     }
2084   }
2085   else {
2086     for (Int_t i = 0; i < nResponses; i++) {
2087       responses[i] = fRandom->Gaus(mean, sigma);
2088     }
2089   }
2090
2091   // If forwarded error code was a warning (error case has been handled before), return a warning
2092   if (errCode == kWarning)
2093     return kWarning;
2094
2095   return ownErrCode; // Forward success/warning
2096 }
2097
2098
2099 //_____________________________________________________________________________
2100 void AliAnalysisTaskPID::PrintSettings(Bool_t printSystematicsSettings) const
2101 {
2102   // Print current settings.
2103
2104   printf("\n\nSettings for task %s:\n", GetName());
2105   printf("Is pPb/Pbp: %d -> %s\n", GetIsPbpOrpPb(), GetIsPbpOrpPb() ? "Adapting vertex cuts" : "Using standard vertex cuts");
2106   printf("Track cuts: %s\n", fTrackFilter ? fTrackFilter->GetTitle() : "-");
2107   printf("Eta cut: %.2f <= |eta| <= %.2f\n", GetEtaAbsCutLow(), GetEtaAbsCutUp());
2108   printf("Phi' cut: %d\n", GetUsePhiCut());
2109   printf("TPCCutMIGeo: %d\n", GetUseTPCCutMIGeo());
2110   if (GetUseTPCCutMIGeo()) {
2111     printf("GetCutGeo: %f\n", GetCutGeo());
2112     printf("GetCutNcr: %f\n", GetCutNcr());
2113     printf("GetCutNcl: %f\n", GetCutNcl());
2114   }
2115   printf("TPCnclCut: %d\n", GetUseTPCnclCut());
2116   if (GetUseTPCnclCut()) {
2117     printf("GetCutPureNcl: %d\n", GetCutPureNcl());
2118   }
2119
2120   printf("\n");
2121
2122   printf("Centrality estimator: %s\n", GetCentralityEstimator().Data());
2123
2124   printf("\n");
2125
2126   printf("Use MC-ID for signal generation: %d\n", GetUseMCidForGeneration());
2127   printf("Use ITS: %d\n", GetUseITS());
2128   printf("Use TOF: %d\n", GetUseTOF());
2129   printf("Use priors: %d\n", GetUsePriors());
2130   printf("Use TPC default priors: %d\n", GetUseTPCDefaultPriors());
2131   printf("Use convoluted Gauss: %d\n", GetUseConvolutedGaus());
2132   printf("Accuracy of non-Gaussian tail: %e\n", GetAccuracyNonGaussianTail());
2133   printf("Take into account muons: %d\n", GetTakeIntoAccountMuons());
2134   printf("TOF mode: %d\n", GetTOFmode());
2135   printf("\nParams for transition from gauss to asymmetric shape:\n");
2136   printf("[0]: %e\n", GetConvolutedGaussTransitionPar(0));
2137   printf("[1]: %e\n", GetConvolutedGaussTransitionPar(1));
2138   printf("[2]: %e\n", GetConvolutedGaussTransitionPar(2));
2139
2140   printf("\n");
2141
2142   printf("Do PID: %d\n", fDoPID);
2143   printf("Do Efficiency: %d\n", fDoEfficiency);
2144   printf("Do PtResolution: %d\n", fDoPtResolution);
2145
2146   printf("\n");
2147
2148   printf("Input from other task: %d\n", GetInputFromOtherTask());
2149   printf("Store additional jet information: %d\n", GetStoreAdditionalJetInformation());
2150   printf("Store centrality percentile: %d", GetStoreCentralityPercentile());
2151
2152   if (printSystematicsSettings)
2153     PrintSystematicsSettings();
2154   else
2155     printf("\n\n\n");
2156 }
2157
2158
2159 //_____________________________________________________________________________
2160 void AliAnalysisTaskPID::PrintSystematicsSettings() const
2161 {
2162   // Print current settings for systematic studies.
2163
2164   printf("\n\nSettings for systematics for task %s:\n", GetName());
2165   printf("SplinesThr:\t%f\n", GetSystematicScalingSplinesThreshold());
2166   printf("SplinesBelowThr:\t%f\n", GetSystematicScalingSplinesBelowThreshold());
2167   printf("SplinesAboveThr:\t%f\n", GetSystematicScalingSplinesAboveThreshold());
2168   printf("EtaCorrMomThr:\t%f\n", GetSystematicScalingEtaCorrectionMomentumThr());
2169   printf("EtaCorrLowP:\t%f\n", GetSystematicScalingEtaCorrectionLowMomenta());
2170   printf("EtaCorrHighP:\t%f\n", GetSystematicScalingEtaCorrectionHighMomenta());
2171   printf("SigmaParaThr:\t%f\n", GetSystematicScalingEtaSigmaParaThreshold());
2172   printf("SigmaParaBelowThr:\t%f\n", GetSystematicScalingEtaSigmaParaBelowThreshold());
2173   printf("SigmaParaAboveThr:\t%f\n", GetSystematicScalingEtaSigmaParaAboveThreshold());
2174   printf("MultCorr:\t%f\n", GetSystematicScalingMultCorrection());
2175   printf("TOF mode: %d\n", GetTOFmode());
2176
2177   printf("\n\n");
2178 }
2179
2180
2181 //_____________________________________________________________________________
2182 Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePDGcode, Double_t centralityPercentile,
2183                                         Double_t jetPt)
2184 {
2185   // Process the track (generate expected response, fill histos, etc.).
2186   // particlePDGcode == 0 means data. Otherwise, the corresponding MC ID will be assumed.
2187
2188   //Printf("Debug: Task %s is starting to process track: dEdx %f, pTPC %f, eta %f, ncl %d\n", GetName(), track->GetTPCsignal(), track->GetTPCmomentum(),
2189   //       track->Eta(), track->GetTPCsignalN());
2190
2191   if(fDebug > 1)
2192     printf("File: %s, Line: %d: ProcessTrack\n", (char*)__FILE__, __LINE__);
2193
2194   if (!fDoPID)
2195     return kFALSE;
2196
2197   if(fDebug > 2)
2198     printf("File: %s, Line: %d: ProcessTrack -> Processing started\n", (char*)__FILE__, __LINE__);
2199
2200   const Bool_t isMC = (particlePDGcode == 0) ? kFALSE : kTRUE;
2201
2202   Int_t binMC = -1;
2203
2204   if (isMC) {
2205     if (TMath::Abs(particlePDGcode) == 211) {//Pion
2206       binMC = 3;
2207     }
2208     else if (TMath::Abs(particlePDGcode) == 321) {//Kaon
2209       binMC = 1;
2210     }
2211     else if (TMath::Abs(particlePDGcode) == 2212) {//Proton
2212       binMC = 4;
2213     }
2214     else if (TMath::Abs(particlePDGcode) == 11) {//Electron
2215       binMC = 0;
2216     }
2217     else if (TMath::Abs(particlePDGcode) == 13) {//Muon
2218       binMC = 2;
2219     }
2220     else // In MC-ID case, set to underflow bin such that the response from this track is only used for unidentified signal generation
2221          // or signal generation with PID response and the track is still there (as in data) - e.g. relevant w.r.t. deuterons.
2222          // This is important to be as much as possible consistent with data. And the tracks can still be removed by disabling the
2223          // underflow bin for the projections
2224       binMC = -1;
2225   }
2226
2227   // Momenta
2228   //Double_t p = track->GetP();
2229   const Double_t pTPC = track->GetTPCmomentum();
2230   Double_t pT = track->Pt();
2231
2232   Double_t z = -1, xi = -1;
2233   GetJetTrackObservables(pT, jetPt, z, xi);
2234
2235
2236   Double_t trackCharge = track->Charge();
2237
2238   // TPC signal
2239   const Bool_t tuneOnDataTPC = fPIDResponse->IsTunedOnData() &&
2240                                ((fPIDResponse->GetTunedOnDataMask() & AliPIDResponse::kDetTPC) == AliPIDResponse::kDetTPC);
2241   Double_t dEdxTPC = tuneOnDataTPC ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
2242
2243   if (dEdxTPC <= 0) {
2244     Printf("Skipping track with strange dEdx value: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), pTPC,
2245            track->Eta(), track->GetTPCsignalN());
2246     return kFALSE;
2247   }
2248
2249   // Completely remove tracks from light nuclei defined by el and pr <dEdx> as function of pT.
2250   // A very loose cut to be sure not to throw away electrons and/or protons
2251   Double_t nSigmaPr = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton);
2252   Double_t nSigmaEl = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
2253
2254   if ((pTPC >= 0.3 && (nSigmaPr > 10 && nSigmaEl > 10)) ||
2255       (pTPC <  0.3 && (nSigmaPr > 15 && nSigmaEl > 15))) {
2256     Printf("Skipping track which seems to be a light nuclei: dEdx %f, pTPC %f, pT %f, eta %f, ncl %d, nSigmaPr %f, nSigmaEl %f\n",
2257            track->GetTPCsignal(), pTPC, pT, track->Eta(), track->GetTPCsignalN(), nSigmaPr, nSigmaEl);
2258     return kFALSE;
2259   }
2260
2261
2262   Double_t dEdxEl, dEdxKa, dEdxPi, dEdxMu, dEdxPr;
2263   Double_t sigmaEl, sigmaKa, sigmaPi, sigmaMu, sigmaPr;
2264
2265   if (fDoAnySystematicStudiesOnTheExpectedSignal) {
2266     // Get the uncorrected signal first and the corresponding correction factors.
2267     // Then modify the correction factors and properly recalculate the corrected dEdx
2268
2269     // Get pure spline values for dEdx_expected, without any correction
2270     dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2271     dEdxKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2272     dEdxPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2273     dEdxMu = !fTakeIntoAccountMuons ? -1 :
2274                             fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2275     dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2276
2277     // Scale splines, if desired
2278     if ((TMath::Abs(fSystematicScalingSplinesBelowThreshold - 1.0) > fgkEpsilon) ||
2279         (TMath::Abs(fSystematicScalingSplinesAboveThreshold - 1.0) > fgkEpsilon)) {
2280
2281       // Tune turn-on of correction for pions (not so relevant for the other species, since at very large momenta)
2282       Double_t bg = 0;
2283       Double_t scaleFactor = 1.;
2284       Double_t usedSystematicScalingSplines = 1.;
2285
2286       bg = pTPC / AliPID::ParticleMass(AliPID::kElectron);
2287       scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2288       usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2289                                             + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2290       dEdxEl *= usedSystematicScalingSplines;
2291
2292       bg = pTPC / AliPID::ParticleMass(AliPID::kKaon);
2293       scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2294       usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2295                                             + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2296       dEdxKa *= usedSystematicScalingSplines;
2297
2298       bg = pTPC / AliPID::ParticleMass(AliPID::kPion);
2299       scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2300       usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2301                                             + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2302       dEdxPi *= usedSystematicScalingSplines;
2303
2304       if (fTakeIntoAccountMuons) {
2305         bg = pTPC / AliPID::ParticleMass(AliPID::kMuon);
2306         scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2307         usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2308                                               + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2309         dEdxMu *= usedSystematicScalingSplines;
2310       }
2311
2312
2313       bg = pTPC / AliPID::ParticleMass(AliPID::kProton);
2314       scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2315       usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2316                                             + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2317       dEdxPr *= usedSystematicScalingSplines;
2318     }
2319
2320     // Get the eta correction factors for the (modified) expected dEdx
2321     Double_t etaCorrEl = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxEl) : 1.;
2322     Double_t etaCorrKa = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxKa) : 1.;
2323     Double_t etaCorrPi = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPi) : 1.;
2324     Double_t etaCorrMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCEtaCorrection() ?
2325                             fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxMu) : 1.;
2326     Double_t etaCorrPr = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPr) : 1.;
2327
2328     // Scale eta correction factors, if desired (and eta correction maps are to be used, otherwise it is not possible!)
2329     if (fPIDResponse->UseTPCEtaCorrection() &&
2330         (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon ||
2331          TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon)) {
2332       // Since we do not want to scale the splines with this, but only the eta variation, only scale the deviation of the correction factor!
2333       // E.g. if we would have a flat eta depence fixed at 1.0, we would shift the whole thing equal to shifting the splines by the same factor!
2334
2335
2336       // Due to additional azimuthal effects, there is an additional eta dependence for low momenta which is not corrected successfully so far.
2337       // One can assign a different (higher) systematic scale factor for this low-p region and a threshold which separates low- and high-p.
2338       // An ERF will be used to get (as a function of p_TPC) from one correction factor to the other within roughly 0.2 GeV/c
2339       Double_t usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionHighMomenta;
2340
2341       if (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - fSystematicScalingEtaCorrectionLowMomenta) > fgkEpsilon) {
2342         const Double_t fractionHighMomentumScaleFactor = 0.5 * (1. + TMath::Erf((pTPC - fSystematicScalingEtaCorrectionMomentumThr) / 0.1));
2343         usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionLowMomenta * (1 - fractionHighMomentumScaleFactor)
2344                                              + fSystematicScalingEtaCorrectionHighMomenta * fractionHighMomentumScaleFactor;
2345       }
2346
2347       Double_t maxEtaVariationEl = GetMaxEtaVariation(dEdxEl);
2348       etaCorrEl = etaCorrEl * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrEl - 1.0) / maxEtaVariationEl);
2349
2350       Double_t maxEtaVariationKa = GetMaxEtaVariation(dEdxKa);
2351       etaCorrKa = etaCorrKa * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrKa - 1.0) / maxEtaVariationKa);
2352
2353       Double_t maxEtaVariationPi = GetMaxEtaVariation(dEdxPi);
2354       etaCorrPi = etaCorrPi * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrPi - 1.0) / maxEtaVariationPi);
2355
2356       if (fTakeIntoAccountMuons) {
2357         Double_t maxEtaVariationMu = GetMaxEtaVariation(dEdxMu);
2358         etaCorrMu = etaCorrMu * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrMu - 1.0) / maxEtaVariationMu);
2359       }
2360       else
2361         etaCorrMu = 1.0;
2362
2363       Double_t maxEtaVariationPr = GetMaxEtaVariation(dEdxPr);
2364       etaCorrPr = etaCorrPr * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrPr - 1.0) / maxEtaVariationPr);
2365
2366
2367       /*OLD
2368       etaCorrEl = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrEl - 1.0);
2369       etaCorrKa = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrKa - 1.0);
2370       etaCorrPi = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPi - 1.0);
2371       etaCorrMu = fTakeIntoAccountMuons ? (1.0 + usedSystematicScalingEtaCorrection * (etaCorrMu - 1.0)) : 1.0;
2372       etaCorrPr = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPr - 1.0);
2373       */
2374     }
2375
2376     // Get the multiplicity correction factors for the (modified) expected dEdx
2377     const Int_t currEvtMultiplicity = fPIDResponse->GetTPCResponse().GetCurrentEventMultiplicity();
2378
2379     Double_t multiplicityCorrEl = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2380                                     dEdxEl * etaCorrEl, currEvtMultiplicity) : 1.;
2381     Double_t multiplicityCorrKa = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2382                                     dEdxKa * etaCorrKa, currEvtMultiplicity) : 1.;
2383     Double_t multiplicityCorrPi = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2384                                     dEdxPi * etaCorrPi, currEvtMultiplicity) : 1.;
2385     Double_t multiplicityCorrMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCMultiplicityCorrection() ?
2386                                     fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track, dEdxMu * etaCorrMu, currEvtMultiplicity) : 1.;
2387     Double_t multiplicityCorrPr = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2388                                     dEdxPr * etaCorrPr, currEvtMultiplicity) : 1.;
2389
2390     Double_t multiplicityCorrSigmaEl = fPIDResponse->UseTPCMultiplicityCorrection() ?
2391                                         fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxEl * etaCorrEl, currEvtMultiplicity) : 1.;
2392     Double_t multiplicityCorrSigmaKa = fPIDResponse->UseTPCMultiplicityCorrection() ?
2393                                         fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxKa * etaCorrKa, currEvtMultiplicity) : 1.;
2394     Double_t multiplicityCorrSigmaPi = fPIDResponse->UseTPCMultiplicityCorrection() ?
2395                                         fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxPi * etaCorrPi, currEvtMultiplicity) : 1.;
2396     Double_t multiplicityCorrSigmaMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCMultiplicityCorrection() ?
2397                                         fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxMu * etaCorrMu, currEvtMultiplicity) : 1.;
2398     Double_t multiplicityCorrSigmaPr = fPIDResponse->UseTPCMultiplicityCorrection() ?
2399                                         fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxPr * etaCorrPr, currEvtMultiplicity) : 1.;
2400
2401     // Scale multiplicity correction factors, if desired (and multiplicity correction functions are to be used, otherwise it is not possible!)
2402     if (fPIDResponse->UseTPCMultiplicityCorrection() && TMath::Abs(fSystematicScalingMultCorrection - 1.0) > fgkEpsilon) {
2403       // Since we do not want to scale the splines with this, but only the multiplicity variation, only scale the deviation of the correction factor!
2404       // E.g. if we would have a flat mult depence fix at 1.0, we would shift the whole thing equal to shifting the splines by the same factor!
2405
2406       multiplicityCorrEl = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrEl - 1.0);
2407       multiplicityCorrKa = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrKa - 1.0);
2408       multiplicityCorrPi = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrPi - 1.0);
2409       multiplicityCorrMu = fTakeIntoAccountMuons ? (1.0 + fSystematicScalingMultCorrection * (multiplicityCorrMu - 1.0)) : 1.0;
2410       multiplicityCorrPr = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrPr - 1.0);
2411     }
2412
2413     // eta correction must be enabled in order to use the new sigma parametrisation maps. Since this is the absolute sigma
2414     // for a track calculated with the unscaled paramaters, we have to devide out dEdxExpectedEtaCorrected and then need
2415     // to scale with the multiplicitySigmaCorrFactor * fSystematicScalingEtaSigmaPara. In the end, one has to scale with the
2416     // (modified) dEdx to get the absolute sigma
2417     // This means there is no extra parameter for the multiplicitySigmaCorrFactor, but only for the sigma map itself.
2418     // This is valid, since it appears only as a product. One has to assume a larger systematic shift in case of additional
2419     // multiplicity dependence....
2420
2421     Bool_t doSigmaSystematics = (TMath::Abs(fSystematicScalingEtaSigmaParaBelowThreshold - 1.0) > fgkEpsilon) ||
2422                                 (TMath::Abs(fSystematicScalingEtaSigmaParaAboveThreshold - 1.0) > fgkEpsilon);
2423
2424
2425     Double_t dEdxElExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2426                                                                                kTRUE, kFALSE);
2427     Double_t systematicScalingEtaSigmaParaEl = 1.;
2428     if (doSigmaSystematics) {
2429       Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxElExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2430       systematicScalingEtaSigmaParaEl = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2431                                         + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2432     }
2433     Double_t sigmaRelEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2434                                                                           kTRUE, kFALSE)
2435                           / dEdxElExpected * systematicScalingEtaSigmaParaEl * multiplicityCorrSigmaEl;
2436
2437
2438     Double_t dEdxKaExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2439                                                                                kTRUE, kFALSE);
2440     Double_t systematicScalingEtaSigmaParaKa = 1.;
2441     if (doSigmaSystematics) {
2442       Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxKaExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2443       systematicScalingEtaSigmaParaKa = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2444                                         + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2445     }
2446     Double_t sigmaRelKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2447                                                                           kTRUE, kFALSE)
2448                           / dEdxKaExpected * systematicScalingEtaSigmaParaKa * multiplicityCorrSigmaKa;
2449
2450
2451     Double_t dEdxPiExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2452                                                                                kTRUE, kFALSE);
2453     Double_t systematicScalingEtaSigmaParaPi = 1.;
2454     if (doSigmaSystematics) {
2455       Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxPiExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2456       systematicScalingEtaSigmaParaPi = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2457                                         + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2458     }
2459     Double_t sigmaRelPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2460                                                                           kTRUE, kFALSE)
2461                           / dEdxPiExpected * systematicScalingEtaSigmaParaPi * multiplicityCorrSigmaPi;
2462
2463
2464     Double_t sigmaRelMu = 999.;
2465     if (fTakeIntoAccountMuons) {
2466       Double_t dEdxMuExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2467                                                                                  kTRUE, kFALSE);
2468       Double_t systematicScalingEtaSigmaParaMu = 1.;
2469       if (doSigmaSystematics) {
2470         Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxMuExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2471         systematicScalingEtaSigmaParaMu = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2472                                           + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2473       }
2474       sigmaRelMu = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2475                    / dEdxMuExpected * systematicScalingEtaSigmaParaMu * multiplicityCorrSigmaMu;
2476     }
2477
2478
2479     Double_t dEdxPrExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2480                                                                                kTRUE, kFALSE);
2481     Double_t systematicScalingEtaSigmaParaPr = 1.;
2482     if (doSigmaSystematics) {
2483       Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxPrExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2484       systematicScalingEtaSigmaParaPr = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2485                                         + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2486     }
2487     Double_t sigmaRelPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2488                                                                           kTRUE, kFALSE)
2489                           / dEdxPrExpected * systematicScalingEtaSigmaParaPr * multiplicityCorrSigmaPr;
2490
2491     // Now scale the (possibly modified) spline values with the (possibly modified) correction factors
2492     dEdxEl *= etaCorrEl * multiplicityCorrEl;
2493     dEdxKa *= etaCorrKa * multiplicityCorrKa;
2494     dEdxPi *= etaCorrPi * multiplicityCorrPi;
2495     dEdxMu *= etaCorrMu * multiplicityCorrMu;
2496     dEdxPr *= etaCorrPr * multiplicityCorrPr;
2497
2498     // Finally, get the absolute sigma
2499     sigmaEl = sigmaRelEl * dEdxEl;
2500     sigmaKa = sigmaRelKa * dEdxKa;
2501     sigmaPi = sigmaRelPi * dEdxPi;
2502     sigmaMu = sigmaRelMu * dEdxMu;
2503     sigmaPr = sigmaRelPr * dEdxPr;
2504
2505   }
2506   else {
2507     // No systematic studies on expected signal - just take it as it comve from the TPCPIDResponse
2508     dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2509                                                               fPIDResponse->UseTPCEtaCorrection(),
2510                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2511     dEdxKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2512                                                               fPIDResponse->UseTPCEtaCorrection(),
2513                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2514     dEdxPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2515                                                               fPIDResponse->UseTPCEtaCorrection(),
2516                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2517     dEdxMu = !fTakeIntoAccountMuons ? -1 :
2518                             fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2519                                                                              fPIDResponse->UseTPCEtaCorrection(),
2520                                                                              fPIDResponse->UseTPCMultiplicityCorrection());
2521     dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2522                                                               fPIDResponse->UseTPCEtaCorrection(),
2523                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2524
2525     sigmaEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2526                                                               fPIDResponse->UseTPCEtaCorrection(),
2527                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2528     sigmaKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2529                                                               fPIDResponse->UseTPCEtaCorrection(),
2530                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2531     sigmaPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2532                                                               fPIDResponse->UseTPCEtaCorrection(),
2533                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2534     sigmaMu = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2535                                                               fPIDResponse->UseTPCEtaCorrection(),
2536                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2537     sigmaPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2538                                                               fPIDResponse->UseTPCEtaCorrection(),
2539                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2540   }
2541
2542   Double_t deltaPrimeElectron = (dEdxEl > 0) ? dEdxTPC / dEdxEl : -1;
2543   if (dEdxEl <= 0)  {
2544     Printf("Error: Expected TPC signal <= 0 for electron hypothesis");
2545     return kFALSE;
2546   }
2547
2548   Double_t deltaPrimeKaon = (dEdxKa > 0) ? dEdxTPC / dEdxKa : -1;
2549   if (dEdxKa <= 0)  {
2550     Printf("Error: Expected TPC signal <= 0 for kaon hypothesis");
2551     return kFALSE;
2552   }
2553
2554   Double_t deltaPrimePion = (dEdxPi > 0) ? dEdxTPC / dEdxPi : -1;
2555   if (dEdxPi <= 0)  {
2556     Printf("Error: Expected TPC signal <= 0 for pion hypothesis");
2557     return kFALSE;
2558   }
2559
2560   Double_t deltaPrimeProton = (dEdxPr > 0) ? dEdxTPC / dEdxPr : -1;
2561   if (dEdxPr <= 0)  {
2562     Printf("Error: Expected TPC signal <= 0 for proton hypothesis");
2563     return kFALSE;
2564   }
2565
2566   if(fDebug > 2)
2567     printf("File: %s, Line: %d: ProcessTrack -> Compute probabilities\n", (char*)__FILE__, __LINE__);
2568
2569
2570   // Use probabilities to weigh the response generation for the different species.
2571   // Also determine the most probable particle type.
2572   Double_t prob[AliPID::kSPECIESC];
2573   for (Int_t i = 0; i < AliPID::kSPECIESC; i++)
2574     prob[i] = 0;
2575
2576   fPIDcombined->ComputeProbabilities(track, fPIDResponse, prob);
2577
2578   // If muons are not to be taken into account, just set their probability to zero and normalise the remaining probabilities later
2579   if (!fTakeIntoAccountMuons)
2580     prob[AliPID::kMuon] = 0;
2581
2582   // Priors might cause problems in case of mismatches, i.e. mismatch particles will be identified as pions.
2583   // Can be easily caught for low p_TPC by just setting pion (and muon) probs to zero, if dEdx is larger than
2584   // expected dEdx of electrons and kaons
2585   if (pTPC < 1. && dEdxTPC > dEdxEl && dEdxTPC > dEdxKa) {
2586     prob[AliPID::kMuon] = 0;
2587     prob[AliPID::kPion] = 0;
2588   }
2589
2590
2591   // Bug: One can set the number of species for PIDcombined, but PIDcombined will call PIDresponse, which writes without testing
2592   // the probs for kSPECIESC (including light nuclei) into the array.
2593   // In this case, when only kSPECIES are considered, the probabilities have to be rescaled!
2594
2595   // Exceptions:
2596   // 1. If the most probable PID is a light nucleus and above expected dEdx + 5 sigma of el to be sure not to mix up things in the
2597   //    high-pT region (also contribution there completely negligable!)
2598   // 2. If dEdx is larger than the expected dEdx + 5 sigma of el and pr, it is most likely a light nucleus
2599   //    (if dEdx is more than 5 sigma away from expectation, flat TPC probs are set... and the light nuclei splines are not
2600   //     too accurate)
2601   // In these cases:
2602   // The highest expected dEdx is either for pr or for el. Assign 100% probability to the species with highest expected dEdx then.
2603   // In all other cases: Throw away light nuclei probs and rescale other probs to 100%.
2604
2605   // Find most probable species for the ID
2606   Int_t maxIndex = TMath::LocMax(AliPID::kSPECIESC, prob);
2607
2608   if ((prob[maxIndex] > 0 && maxIndex >= AliPID::kSPECIES && dEdxTPC > dEdxEl + 5. * sigmaEl) ||
2609       (dEdxTPC > dEdxEl + 5. * sigmaEl && dEdxTPC > dEdxPr + 5. * sigmaPr)) {
2610     for (Int_t i = 0; i < AliPID::kSPECIESC; i++)
2611       prob[i] = 0;
2612
2613     if (dEdxEl > dEdxPr) {
2614       prob[AliPID::kElectron] = 1.;
2615       maxIndex = AliPID::kElectron;
2616     }
2617     else {
2618       prob[AliPID::kProton] = 1.;
2619       maxIndex = AliPID::kProton;
2620     }
2621   }
2622   else {
2623     for (Int_t i = AliPID::kSPECIES; i < AliPID::kSPECIESC; i++)
2624       prob[i] = 0;
2625
2626     Double_t probSum = 0.;
2627     for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2628       probSum += prob[i];
2629
2630     if (probSum > 0) {
2631       for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2632         prob[i] /= probSum;
2633     }
2634
2635     // If all probabilities are equal (should only happen if no priors are used), set most probable PID to pion
2636     // because LocMax returns just the first maximal value (i.e. AliPID::kElectron)
2637     if (TMath::Abs(prob[maxIndex] - 1. / AliPID::kSPECIES) < 1e-5)
2638       maxIndex = AliPID::kPion;
2639
2640     // In all other cases, the maximum remains untouched from the scaling (and is < AliPID::kSPECIES by construction)
2641   }
2642
2643   if (!isMC) {
2644     // Clean up the most probable PID for low momenta (not an issue for the probabilities, since fractions small,
2645     // but to get clean (and nice looking) templates (if most probable PID is used instead of generated responses), it should be done).
2646     // Problem: If more than 5 sigma away (and since priors also included), most probable PID for dEdx around protons could be pions.
2647     // Idea: Only clean selection for K and p possible. So, just take for the most probable PID the probs from TPC only calculated
2648     // by hand without restriction to 5 sigma and without priors. This will not affect the template generation, but afterwards one
2649     // can replace the K and p templates with the most probable PID distributions, so all other species templates remain the same.
2650     // I.e. it doesn't hurt if the most probable PID is distributed incorrectly between el, pi, mu.
2651     Bool_t maxIndexSetManually = kFALSE;
2652
2653     if (pTPC < 1.) {
2654       Double_t probManualTPC[AliPID::kSPECIES];
2655        for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2656         probManualTPC[i] = 0;
2657
2658       probManualTPC[AliPID::kElectron] = TMath::Exp(-0.5*(dEdxTPC-dEdxEl)*(dEdxTPC-dEdxEl)/(sigmaEl*sigmaEl))/sigmaEl;
2659       // Muons are not important here, just ignore and save processing time
2660       probManualTPC[AliPID::kMuon] = 0.;//TMath::Exp(-0.5*(dEdxTPC-dEdxMu)*(dEdxTPC-dEdxMu)/(sigmaMu*sigmaMu))/sigmaMu;
2661       probManualTPC[AliPID::kPion] = TMath::Exp(-0.5*(dEdxTPC-dEdxPi)*(dEdxTPC-dEdxPi)/(sigmaPi*sigmaPi))/sigmaPi;
2662       probManualTPC[AliPID::kKaon] = TMath::Exp(-0.5*(dEdxTPC-dEdxKa)*(dEdxTPC-dEdxKa)/(sigmaKa*sigmaKa))/sigmaKa;
2663       probManualTPC[AliPID::kProton] = TMath::Exp(-0.5*(dEdxTPC-dEdxPr)*(dEdxTPC-dEdxPr)/(sigmaPr*sigmaPr))/sigmaPr;
2664
2665       const Int_t maxIndexManualTPC = TMath::LocMax(AliPID::kSPECIES, probManualTPC);
2666       // Should always be > 0, but in case the dEdx is far off such that the machine accuracy sets every prob to zero,
2667       // better take the "old" result
2668       if (probManualTPC[maxIndexManualTPC] > 0.)
2669         maxIndex = maxIndexManualTPC;
2670
2671       maxIndexSetManually = kTRUE;
2672     }
2673
2674
2675     // Translate from AliPID numbering to numbering of this class
2676     if (prob[maxIndex] > 0 || maxIndexSetManually) {
2677       if (maxIndex == AliPID::kElectron)
2678         binMC = 0;
2679       else if (maxIndex == AliPID::kKaon)
2680         binMC = 1;
2681       else if (maxIndex == AliPID::kMuon)
2682         binMC = 2;
2683       else if (maxIndex == AliPID::kPion)
2684         binMC = 3;
2685       else if (maxIndex == AliPID::kProton)
2686         binMC = 4;
2687       else
2688         binMC = -1;
2689     }
2690     else {
2691       // Only take track into account for expectation values, if valid pid response is available.. Otherwise: Set to underflow bin.
2692       binMC = -1;
2693     }
2694   }
2695
2696   /*
2697   //For testing: Swap p<->pT to analyse pure p-dependence => Needs to be removed later
2698   Double_t temp = pT;
2699   pT = pTPC;
2700   pTPC = temp;
2701   */
2702
2703   TOFpidInfo tofPIDinfo = GetTOFType(track, fTOFmode);
2704
2705   Double_t entry[fStoreAdditionalJetInformation ? kDataNumAxes : kDataNumAxes - fgkNumJetAxes];
2706   entry[kDataMCID] = binMC;
2707   entry[kDataSelectSpecies] = 0;
2708   entry[kDataPt] = pT;
2709   entry[kDataDeltaPrimeSpecies] = deltaPrimeElectron;
2710   entry[kDataCentrality] = centralityPercentile;
2711
2712   if (fStoreAdditionalJetInformation) {
2713     entry[kDataJetPt] = jetPt;
2714     entry[kDataZ] = z;
2715     entry[kDataXi] = xi;
2716   }
2717
2718   entry[GetIndexOfChargeAxisData()] = trackCharge;
2719   entry[GetIndexOfTOFpidInfoAxisData()] = tofPIDinfo;
2720
2721   fhPIDdataAll->Fill(entry);
2722
2723   entry[kDataSelectSpecies] = 1;
2724   entry[kDataDeltaPrimeSpecies] = deltaPrimeKaon;
2725   fhPIDdataAll->Fill(entry);
2726
2727   entry[kDataSelectSpecies] = 2;
2728   entry[kDataDeltaPrimeSpecies] = deltaPrimePion;
2729   fhPIDdataAll->Fill(entry);
2730
2731   entry[kDataSelectSpecies] = 3;
2732   entry[kDataDeltaPrimeSpecies] = deltaPrimeProton;
2733   fhPIDdataAll->Fill(entry);
2734
2735
2736   // Construct the expected shape for the transition p -> pT
2737
2738   Double_t genEntry[fStoreAdditionalJetInformation ? kGenNumAxes : kGenNumAxes - fgkNumJetAxes];
2739   genEntry[kGenMCID] = binMC;
2740   genEntry[kGenSelectSpecies] = 0;
2741   genEntry[kGenPt] = pT;
2742   genEntry[kGenDeltaPrimeSpecies] = -999;
2743   genEntry[kGenCentrality] = centralityPercentile;
2744
2745   if (fStoreAdditionalJetInformation) {
2746     genEntry[kGenJetPt] = jetPt;
2747     genEntry[kGenZ] = z;
2748     genEntry[kGenXi] = xi;
2749   }
2750
2751   genEntry[GetIndexOfChargeAxisGen()] = trackCharge;
2752   genEntry[GetIndexOfTOFpidInfoAxisGen()] = tofPIDinfo;
2753
2754   // Generate numGenEntries "responses" with fluctuations around the expected values.
2755   // fgkMaxNumGenEntries = 500 turned out to give reasonable templates even for highest track pT in 15-20 GeV/c jet pT bin.
2756   Int_t numGenEntries = fgkMaxNumGenEntries; // fgkMaxNumGenEntries = 500
2757
2758   /*OLD: Different number of responses depending on pT and jet pT for fgkMaxNumGenEntries = 1000
2759    * => Problem: If threshold to higher number of responses inside a bin (or after rebinning), then the template
2760    * is biased to the higher pT.
2761   // Generate numGenEntries "responses" with fluctuations around the expected values.
2762   // The higher the (transverse) momentum, the more "responses" will be generated in order not to run out of statistics too fast.
2763   Int_t numGenEntries = 10;
2764
2765   // Jets have even worse statistics, therefore, scale numGenEntries further
2766   if (jetPt >= 40)
2767     numGenEntries *= 20;
2768   else if (jetPt >= 20)
2769     numGenEntries *= 10;
2770   else if (jetPt >= 10)
2771     numGenEntries *= 2;
2772
2773
2774
2775   // Do not generate more entries than available in memory!
2776   if (numGenEntries > fgkMaxNumGenEntries)// fgkMaxNumGenEntries = 1000
2777     numGenEntries = fgkMaxNumGenEntries;
2778   */
2779
2780
2781   ErrorCode errCode = kNoErrors;
2782
2783   if(fDebug > 2)
2784     printf("File: %s, Line: %d: ProcessTrack -> Generate Responses\n", (char*)__FILE__, __LINE__);
2785
2786   // Electrons
2787   errCode = GenerateDetectorResponse(errCode, 1.,              sigmaEl / dEdxEl, fGenRespElDeltaPrimeEl, numGenEntries);
2788   errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxKa, sigmaEl / dEdxKa, fGenRespElDeltaPrimeKa, numGenEntries);
2789   errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxPi, sigmaEl / dEdxPi, fGenRespElDeltaPrimePi, numGenEntries);
2790   errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxPr, sigmaEl / dEdxPr, fGenRespElDeltaPrimePr, numGenEntries);
2791
2792   // Kaons
2793   errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxEl, sigmaKa / dEdxEl, fGenRespKaDeltaPrimeEl, numGenEntries);
2794   errCode = GenerateDetectorResponse(errCode, 1.,              sigmaKa / dEdxKa, fGenRespKaDeltaPrimeKa, numGenEntries);
2795   errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxPi, sigmaKa / dEdxPi, fGenRespKaDeltaPrimePi, numGenEntries);
2796   errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxPr, sigmaKa / dEdxPr, fGenRespKaDeltaPrimePr, numGenEntries);
2797
2798   // Pions
2799   errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxEl, sigmaPi / dEdxEl, fGenRespPiDeltaPrimeEl, numGenEntries);
2800   errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxKa, sigmaPi / dEdxKa, fGenRespPiDeltaPrimeKa, numGenEntries);
2801   errCode = GenerateDetectorResponse(errCode, 1.,              sigmaPi / dEdxPi, fGenRespPiDeltaPrimePi, numGenEntries);
2802   errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxPr, sigmaPi / dEdxPr, fGenRespPiDeltaPrimePr, numGenEntries);
2803
2804   // Muons, if desired
2805   if (fTakeIntoAccountMuons) {
2806     errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxEl, sigmaMu / dEdxEl, fGenRespMuDeltaPrimeEl, numGenEntries);
2807     errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxKa, sigmaMu / dEdxKa, fGenRespMuDeltaPrimeKa, numGenEntries);
2808     errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxPi, sigmaMu / dEdxPi, fGenRespMuDeltaPrimePi, numGenEntries);
2809     errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxPr, sigmaMu / dEdxPr, fGenRespMuDeltaPrimePr, numGenEntries);
2810   }
2811
2812   // Protons
2813   errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxEl, sigmaPr / dEdxEl, fGenRespPrDeltaPrimeEl, numGenEntries);
2814   errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxKa, sigmaPr / dEdxKa, fGenRespPrDeltaPrimeKa, numGenEntries);
2815   errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxPi, sigmaPr / dEdxPi, fGenRespPrDeltaPrimePi, numGenEntries);
2816   errCode = GenerateDetectorResponse(errCode, 1.,              sigmaPr / dEdxPr, fGenRespPrDeltaPrimePr, numGenEntries);
2817
2818   if (errCode != kNoErrors) {
2819     if (errCode == kWarning && fDebug > 1) {
2820       Printf("Warning: Questionable detector response for track, most likely due to very low number of PID clusters! Debug output (dEdx_expected, sigma_expected):");
2821     }
2822     else
2823       Printf("Error: Failed to generate detector response for track - skipped! Debug output (dEdx_expected, sigma_expected):");
2824
2825     if (fDebug > 1) {
2826       Printf("Pr: %e, %e", dEdxPr, sigmaPr);
2827       Printf("Pi: %e, %e", dEdxPi, sigmaPi);
2828       Printf("El: %e, %e", dEdxEl, sigmaEl);
2829       Printf("Mu: %e, %e", dEdxMu, sigmaMu);
2830       Printf("Ka: %e, %e", dEdxKa, sigmaKa);
2831       Printf("track: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), track->GetTPCmomentum(), track->Eta(),
2832             track->GetTPCsignalN());
2833     }
2834
2835     if (errCode != kWarning) {
2836       fhSkippedTracksForSignalGeneration->Fill(track->GetTPCmomentum(), track->GetTPCsignalN());
2837       return kFALSE;// Skip generated response in case of error
2838     }
2839   }
2840
2841   for (Int_t n = 0; n < numGenEntries; n++)  {
2842     if (!isMC || !fUseMCidForGeneration) {
2843       // If no MC info is available or shall not be used, use weighting with priors to generate entries for the different species
2844       Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
2845
2846       // Consider generated response as originating from...
2847       if (rnd <= prob[AliPID::kElectron])
2848         genEntry[kGenMCID] = 0; // ... an electron
2849       else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon])
2850         genEntry[kGenMCID] = 1;  // ... a kaon
2851       else if (rnd <=  prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon])
2852         genEntry[kGenMCID] = 2;  // ... a muon -> NOTE: prob[AliPID::kMuon] = 0 in case of fTakeIntoAccountMuons = kFALSE
2853       else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon] + prob[AliPID::kPion])
2854         genEntry[kGenMCID] = 3;  // ... a pion
2855       else
2856         genEntry[kGenMCID] = 4;  // ... a proton
2857     }
2858
2859     // Electrons
2860     genEntry[kGenSelectSpecies] = 0;
2861     genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeEl[n];
2862     fhGenEl->Fill(genEntry);
2863
2864     genEntry[kGenSelectSpecies] = 1;
2865     genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeKa[n];
2866     fhGenEl->Fill(genEntry);
2867
2868     genEntry[kGenSelectSpecies] = 2;
2869     genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePi[n];
2870     fhGenEl->Fill(genEntry);
2871
2872     genEntry[kGenSelectSpecies] = 3;
2873     genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePr[n];
2874     fhGenEl->Fill(genEntry);
2875
2876     // Kaons
2877     genEntry[kGenSelectSpecies] = 0;
2878     genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeEl[n];
2879     fhGenKa->Fill(genEntry);
2880
2881     genEntry[kGenSelectSpecies] = 1;
2882     genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeKa[n];
2883     fhGenKa->Fill(genEntry);
2884
2885     genEntry[kGenSelectSpecies] = 2;
2886     genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePi[n];
2887     fhGenKa->Fill(genEntry);
2888
2889     genEntry[kGenSelectSpecies] = 3;
2890     genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePr[n];
2891     fhGenKa->Fill(genEntry);
2892
2893     // Pions
2894     genEntry[kGenSelectSpecies] = 0;
2895     genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeEl[n];
2896     fhGenPi->Fill(genEntry);
2897
2898     genEntry[kGenSelectSpecies] = 1;
2899     genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeKa[n];
2900     fhGenPi->Fill(genEntry);
2901
2902     genEntry[kGenSelectSpecies] = 2;
2903     genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePi[n];
2904     fhGenPi->Fill(genEntry);
2905
2906     genEntry[kGenSelectSpecies] = 3;
2907     genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePr[n];
2908     fhGenPi->Fill(genEntry);
2909
2910     if (fTakeIntoAccountMuons) {
2911       // Muons, if desired
2912       genEntry[kGenSelectSpecies] = 0;
2913       genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeEl[n];
2914       fhGenMu->Fill(genEntry);
2915
2916       genEntry[kGenSelectSpecies] = 1;
2917       genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeKa[n];
2918       fhGenMu->Fill(genEntry);
2919
2920       genEntry[kGenSelectSpecies] = 2;
2921       genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePi[n];
2922       fhGenMu->Fill(genEntry);
2923
2924       genEntry[kGenSelectSpecies] = 3;
2925       genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePr[n];
2926       fhGenMu->Fill(genEntry);
2927     }
2928
2929     // Protons
2930     genEntry[kGenSelectSpecies] = 0;
2931     genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeEl[n];
2932     fhGenPr->Fill(genEntry);
2933
2934     genEntry[kGenSelectSpecies] = 1;
2935     genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeKa[n];
2936     fhGenPr->Fill(genEntry);
2937
2938     genEntry[kGenSelectSpecies] = 2;
2939     genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePi[n];
2940     fhGenPr->Fill(genEntry);
2941
2942     genEntry[kGenSelectSpecies] = 3;
2943     genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePr[n];
2944     fhGenPr->Fill(genEntry);
2945   }
2946
2947   if(fDebug > 2)
2948     printf("File: %s, Line: %d: ProcessTrack -> Done\n", (char*)__FILE__, __LINE__);
2949
2950   return kTRUE;
2951 }
2952
2953
2954 //_____________________________________________________________________________
2955 Bool_t AliAnalysisTaskPID::SetConvolutedGaussLambdaParameter(Double_t lambda)
2956 {
2957   // Set the lambda parameter of the convolution to the desired value. Automatically
2958   // calculates the parameters for the transition (restricted) gauss -> convoluted gauss.
2959   fConvolutedGaussTransitionPars[2] = lambda;
2960
2961   // Save old parameters and settings of function to restore them later:
2962   Double_t* oldFuncParams = new Double_t[fkConvolutedGausNPar];
2963   for (Int_t i = 0; i < fkConvolutedGausNPar; i++)
2964     oldFuncParams[i] = fConvolutedGausDeltaPrime->GetParameter(i);
2965   Int_t oldFuncNpx = fConvolutedGausDeltaPrime->GetNpx();
2966   Double_t oldFuncRangeLow = 0, oldFuncRangeUp = 100;
2967   fConvolutedGausDeltaPrime->GetRange(oldFuncRangeLow, oldFuncRangeUp);
2968
2969   // Choose some sufficiently large range
2970   const Double_t rangeStart = 0.5;
2971   const Double_t rangeEnd = 2.0;
2972
2973   // To get the parameters for the transition, just choose arbitrarily input parameters for mu and sigma
2974   // (it makes sense to choose typical values). The ratio sigma_gauss / sigma_convolution is independent
2975   // of mu and as well the difference mu_gauss - mu_convolution.
2976   Double_t muInput = 1.0;
2977   Double_t sigmaInput = fgkSigmaReferenceForTransitionPars;
2978
2979
2980   // Step 1: Generate distribution with input parameters
2981   const Int_t nPar = 3;
2982   Double_t inputPar[nPar] = { muInput, sigmaInput, lambda };
2983
2984   TH1D* hInput = new TH1D("hInput", "Input distribution", 2000, rangeStart, rangeEnd);
2985
2986   fConvolutedGausDeltaPrime->SetParameters(inputPar);
2987   fConvolutedGausDeltaPrime->SetRange(rangeStart, rangeEnd);
2988   fConvolutedGausDeltaPrime->SetNpx(2000);
2989
2990   /*OLD
2991   // The resolution and mean of the AliPIDResponse are determined in finite intervals
2992   // of ncl (also second order effects due to finite dEdx and tanTheta).
2993   // Take this into account for the transition parameters via assuming an approximately flat
2994   // distritubtion in ncl in this interval.
2995   // NOTE: The ncl interval should be the same as the one used for the sigma map creation!
2996   const Int_t nclStart = 151;
2997   const Int_t nclEnd = 160;
2998   const Int_t nclSteps = (nclEnd - nclStart) + 1;
2999   for (Int_t ncl = nclStart; ncl <= nclEnd; ncl++) {
3000     // Resolution scales with 1/sqrt(ncl)
3001     fConvolutedGausDeltaPrime->SetParameter(1, inputPar[1] * sqrt(nclEnd) / sqrt(ncl));
3002     TH1* hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
3003
3004     for (Int_t i = 0; i < 50000000 / nclSteps; i++)
3005       hInput->Fill(hProbDensity->GetRandom());
3006   }
3007   */
3008
3009   TH1* hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
3010
3011   for (Int_t i = 0; i < 50000000; i++)
3012     hInput->Fill(hProbDensity->GetRandom());
3013
3014   // Step 2: Fit generated distribution with restricted gaussian
3015   Int_t maxBin = hInput->GetMaximumBin();
3016   Double_t max = hInput->GetBinContent(maxBin);
3017
3018   UChar_t usedBins = 1;
3019   if (maxBin > 1) {
3020     max += hInput->GetBinContent(maxBin - 1);
3021     usedBins++;
3022   }
3023   if (maxBin < hInput->GetNbinsX()) {
3024     max += hInput->GetBinContent(maxBin + 1);
3025     usedBins++;
3026   }
3027   max /= usedBins;
3028
3029   // NOTE: The range (<-> fraction of maximum) should be the same
3030   // as for the spline and eta maps creation
3031   const Double_t lowThreshold = hInput->GetXaxis()->GetBinLowEdge(hInput->FindFirstBinAbove(0.1 * max));
3032   const Double_t highThreshold = hInput->GetXaxis()->GetBinUpEdge(hInput->FindLastBinAbove(0.1 * max));
3033
3034   TFitResultPtr fitResGaussFirstStep = hInput->Fit("gaus", "QNRS", "", lowThreshold, highThreshold);
3035
3036   TFitResultPtr fitResGauss;
3037
3038   if ((Int_t)fitResGaussFirstStep == 0) {
3039     TF1 fGauss("fGauss", "[0]*TMath::Gaus(x, [1], [2], 1)", rangeStart, rangeEnd);
3040     fGauss.SetParameter(0, fitResGaussFirstStep->GetParams()[0]);
3041     fGauss.SetParError(0, fitResGaussFirstStep->GetErrors()[0]);
3042     fGauss.SetParameter(2, fitResGaussFirstStep->GetParams()[2]);
3043     fGauss.SetParError(2, fitResGaussFirstStep->GetErrors()[2]);
3044
3045     fGauss.FixParameter(1, fitResGaussFirstStep->GetParams()[1]);
3046     fitResGauss = hInput->Fit(&fGauss, "QNS", "s", rangeStart, rangeEnd);
3047   }
3048   else {
3049     fitResGauss = hInput->Fit("gaus", "QNRS", "same", rangeStart, rangeEnd);
3050   }
3051   //OLD TFitResultPtr fitResGauss = hInput->Fit("gaus", "QNRS", "", hInput->GetXaxis()->GetBinLowEdge(hInput->FindFirstBinAbove(0.1 * max)),
3052   //                                        hInput->GetXaxis()->GetBinUpEdge(hInput->FindLastBinAbove(0.1 * max)));
3053
3054
3055   // Step 3: Use parameters from gaussian fit to obtain parameters for the transition "restricted gauss" -> "convoluted gauss"
3056
3057   // 3.1 The ratio sigmaInput / sigma_gaussFit ONLY depends on lambda (which is fixed per period) -> Calculate this first
3058   // for an arbitrary (here: typical) sigma. The ratio is then ~the same for ALL sigma for given lambda!
3059   if ((Int_t)fitResGauss != 0) {
3060     AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Gauss Fit failed!\n");
3061     fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
3062     fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
3063     fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
3064
3065     delete hInput;
3066     delete [] oldFuncParams;
3067
3068     return kFALSE;
3069   }
3070
3071   if (fitResGauss->GetParams()[2] <= 0) {
3072     AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Sigma of gauss fit <= 0!\n");
3073     fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
3074     fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
3075     fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
3076
3077     delete hInput;
3078     delete [] oldFuncParams;
3079
3080     return kFALSE;
3081   }
3082
3083   // sigma correction factor
3084   fConvolutedGaussTransitionPars[1] = sigmaInput / fitResGauss->GetParams()[2];
3085
3086   // 3.2 Now that sigma und lambda are determined, one can calculate mu by shifting the maximum to the desired position,
3087   // i.e. the maximum of the input distribution should coincide with that of the re-calculated distribution,
3088   // which is achieved by getting the same mu for the same sigma.
3089   // NOTE: For fixed lambda, the shift is proportional to sigma and independent of mu!
3090   // So, one can calculate the shift for an arbitrary fixed (here: typical)
3091   // sigma and then simply use this shift for any other sigma by scaling it correspondingly!!!
3092
3093   // Mu shift correction:
3094   // Shift in mu (difference between mean of gauss and mean of convolution) is proportional to sigma!
3095   // Thus, choose a reference sigma (typical -> 0.05), such that for arbitrary sigma one can simple scale
3096   // this shift correction with sigma / referenceSigma.
3097   fConvolutedGaussTransitionPars[0] = (fitResGauss->GetParams()[1] - muInput);
3098
3099
3100   /*Changed on 03.07.2013
3101
3102   // Maximum of fConvolutedGausDeltaPrime should agree with maximum of input
3103   Double_t par[nPar] = { fitResGauss->GetParams()[1], // just as a guess of the maximum position
3104                          sigmaInput,
3105                          lambda };
3106
3107   fConvolutedGausDeltaPrime->SetParameters(par);
3108
3109   Double_t maxXInput = fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001, muInput - 3. * sigmaInput),
3110                                                               muInput + 10. * sigmaInput,
3111                                                               0.0001);
3112
3113   // Maximum shifts proportional to sigma and is linear in mu (~mean of gauss)
3114   // Maximum should be typically located within [gaussMean, gaussMean + 3 gaussSigma].
3115   // -> Larger search range for safety reasons (also: sigma and/or mean might be not 100% accurate).
3116   Double_t maxXconvoluted = fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001,
3117                                                                               fitResGauss->GetParams()[1] - 3. * fitResGauss->GetParams()[2]),
3118                                                                    fitResGauss->GetParams()[1] + 10. * fitResGauss->GetParams()[2],
3119                                                                    0.0001);
3120   if (maxXconvoluted <= 0) {
3121     AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Maximum of fConvolutedGausDeltaPrime <= 0!\n");
3122
3123     fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
3124     fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
3125     fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
3126
3127     delete hInput;
3128     delete [] oldFuncParams;
3129
3130     return kFALSE;
3131   }
3132
3133   // maxX perfectly shifts as par[0] (scaled by sigma) -> Can shift maxX to input value.
3134   // Mu shift correction:
3135   fConvolutedGaussTransitionPars[0] = maxXconvoluted - maxXInput;
3136   */
3137
3138
3139
3140   fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
3141   fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
3142   fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
3143
3144   delete hInput;
3145   delete [] oldFuncParams;
3146
3147   return kTRUE;
3148 }
3149
3150
3151 //_____________________________________________________________________________
3152 AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::SetParamsForConvolutedGaus(Double_t gausMean, Double_t gausSigma)
3153 {
3154   // Set parameters for convoluted gauss using parameters for a pure gaussian.
3155   // If SetConvolutedGaussLambdaParameter has not been called before to initialise the translation parameters,
3156   // some default parameters will be used and an error will show up.
3157
3158   if(fDebug > 1)
3159     printf("File: %s, Line: %d: SetParamsForConvolutedGaus: mean %e, sigma %e\n", (char*)__FILE__, __LINE__, gausMean, gausSigma);
3160
3161   if (fConvolutedGaussTransitionPars[1] < -998) {
3162     AliError("Transition parameters not initialised! Default parameters will be used. Please call SetConvolutedGaussLambdaParameter(...) before any calculations!");
3163     SetConvolutedGaussLambdaParameter(2.0);
3164     AliError(Form("Parameters set to:\n[0]: %f\n[1]: %f\n[2]: %f\n", fConvolutedGaussTransitionPars[0],
3165                   fConvolutedGaussTransitionPars[1], fConvolutedGaussTransitionPars[2]));
3166   }
3167
3168   Double_t par[fkConvolutedGausNPar];
3169   par[2] = fConvolutedGaussTransitionPars[2];
3170   par[1] = fConvolutedGaussTransitionPars[1] * gausSigma;
3171   // maxX perfectly shifts as par[0] (scaled by sigma) -> Can shift maxX so that it sits at the right place.
3172   par[0] = gausMean - fConvolutedGaussTransitionPars[0] * par[1] / fgkSigmaReferenceForTransitionPars;
3173
3174   ErrorCode errCode = kNoErrors;
3175   fConvolutedGausDeltaPrime->SetParameters(par);
3176
3177   if(fDebug > 3)
3178     printf("File: %s, Line: %d: SetParamsForConvolutedGaus -> Parameters set to: %e, %e, %e (transition pars: %e, %e, %e, %e)\n",
3179            (char*)__FILE__, __LINE__, par[0], par[1], par[2], fConvolutedGaussTransitionPars[0], fConvolutedGaussTransitionPars[1],
3180            fConvolutedGaussTransitionPars[2], fgkSigmaReferenceForTransitionPars);
3181
3182   fConvolutedGausDeltaPrime->SetNpx(20); // Small value speeds up following algorithm (valid, since extrema far apart)
3183
3184   // Accuracy of 10^-5 is enough to get 0.1% precise peak for MIPS w.r.t. to dEdx = 2000 of protons
3185   // (should boost up the algorithm, because 10^-10 is the default value!)
3186   Double_t maxX= fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001, gausMean - 2. * gausSigma),
3187                                                         gausMean + 6. * gausSigma, 1.0E-5);
3188
3189   const Double_t maximum = fConvolutedGausDeltaPrime->Eval(maxX);
3190   const Double_t maximumFraction = maximum * fAccuracyNonGaussianTail;
3191
3192   // Estimate lower boundary for subsequent search:
3193   Double_t lowBoundSearchBoundLow = TMath::Max(1e-4, maxX - 5. * gausSigma);
3194   Double_t lowBoundSearchBoundUp = maxX;
3195
3196   Bool_t lowerBoundaryFixedAtZero = kFALSE;
3197
3198   while (fConvolutedGausDeltaPrime->Eval(lowBoundSearchBoundLow) >= maximumFraction) {
3199     if (lowBoundSearchBoundLow <= 0) {
3200       // This should only happen to low dEdx particles with very few clusters and therefore large sigma, such that the gauss goes below zero deltaPrime
3201       if (maximum <= 0) { // Something is weired
3202         printf("Error generating signal: maximum is <= 0!\n");
3203         return kError;
3204       }
3205       else {
3206         const Double_t valueAtZero = fConvolutedGausDeltaPrime->Eval(0);
3207         if (valueAtZero / maximum > 0.05) {
3208           // Too large fraction below zero deltaPrime. Signal generation cannot be reliable in this case
3209           printf("Error generating signal: Too large fraction below zero deltaPrime: convGauss(0) / convGauss(max) =  %e / %e = %e!\n",
3210                  valueAtZero, maximum, valueAtZero / maximum);
3211           return kError;
3212         }
3213       }
3214
3215       /*
3216       printf("Warning: LowBoundSearchBoundLow gets smaller zero -> Set left boundary to zero! Debug output: maximumFraction * fAccuracyNonGaussianTail = %e * %e = %e maxX %f, par[0] %f, par[1] %f, par[2] %f, gausMean %f, gausSigma %f\n",
3217              fConvolutedGausDeltaPrime->Eval(maxX), fAccuracyNonGaussianTail, maximumFraction, maxX, par[0], par[1], par[2], gausMean, gausSigma);
3218       */
3219
3220       lowerBoundaryFixedAtZero = kTRUE;
3221
3222       if (errCode != kError)
3223         errCode = kWarning;
3224
3225       break;
3226     }
3227
3228     lowBoundSearchBoundUp -= gausSigma;
3229     lowBoundSearchBoundLow -= gausSigma;
3230
3231     if (lowBoundSearchBoundLow < 0) {
3232       lowBoundSearchBoundLow = 0;
3233       lowBoundSearchBoundUp += gausSigma;
3234     }
3235   }
3236
3237   // Determine lower boundary inside estimated range. For small values of the maximum: Need more precision, since finer binning!
3238   Double_t rangeStart = lowerBoundaryFixedAtZero ? 0 :
3239                         fConvolutedGausDeltaPrime->GetX(maximumFraction, lowBoundSearchBoundLow, lowBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
3240
3241   // .. and the same for the upper boundary
3242   Double_t rangeEnd = 0;
3243   // If distribution starts beyond upper boundary, everything ends up in the overflow bin. So, just reduce range and Npx to minimum
3244   if (rangeStart > fkDeltaPrimeUpLimit) {
3245     rangeEnd = rangeStart + 0.00001;
3246     fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
3247     fConvolutedGausDeltaPrime->SetNpx(4);
3248   }
3249   else {
3250     // Estimate upper boundary for subsequent search:
3251     Double_t upBoundSearchBoundUp = maxX + 5 * gausSigma;
3252     Double_t upBoundSearchBoundLow = maxX;
3253     while (fConvolutedGausDeltaPrime->Eval(upBoundSearchBoundUp) >= maximumFraction) {
3254       upBoundSearchBoundUp += gausSigma;
3255       upBoundSearchBoundLow += gausSigma;
3256     }
3257
3258     //  For small values of the maximum: Need more precision, since finer binning!
3259     rangeEnd = fConvolutedGausDeltaPrime->GetX(maximumFraction, upBoundSearchBoundLow, upBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
3260
3261     fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
3262     fConvolutedGausDeltaPrime->SetNpx(fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindBin(rangeEnd)
3263                                       - fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindBin(rangeStart) + 1);
3264     //fConvolutedGausDeltaPrime->SetNpx((rangeEnd - rangeStart) / fDeltaPrimeBinWidth + 1);
3265   }
3266
3267   if(fDebug > 3)
3268     printf("File: %s, Line: %d: SetParamsForConvolutedGaus -> range %f - %f, error code %d\n", (char*)__FILE__, __LINE__,
3269            rangeStart, rangeEnd, errCode);
3270
3271   return errCode;
3272 }
3273
3274
3275 //________________________________________________________________________
3276 void AliAnalysisTaskPID::SetUpGenHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const
3277 {
3278   // Sets bin limits for axes which are not standard binned and the axes titles.
3279
3280   hist->SetBinEdges(kGenPt, binsPt);
3281   hist->SetBinEdges(kGenDeltaPrimeSpecies, binsDeltaPrime);
3282   hist->SetBinEdges(kGenCentrality, binsCent);
3283
3284   if (fStoreAdditionalJetInformation)
3285     hist->SetBinEdges(kGenJetPt, binsJetPt);
3286
3287   // Set axes titles
3288   hist->GetAxis(kGenMCID)->SetTitle("MC PID");
3289   hist->GetAxis(kGenMCID)->SetBinLabel(1, "e");
3290   hist->GetAxis(kGenMCID)->SetBinLabel(2, "K");
3291   hist->GetAxis(kGenMCID)->SetBinLabel(3, "#mu");
3292   hist->GetAxis(kGenMCID)->SetBinLabel(4, "#pi");
3293   hist->GetAxis(kGenMCID)->SetBinLabel(5, "p");
3294
3295   hist->GetAxis(kGenSelectSpecies)->SetTitle("Select Species");
3296   hist->GetAxis(kGenSelectSpecies)->SetBinLabel(1, "e");
3297   hist->GetAxis(kGenSelectSpecies)->SetBinLabel(2, "K");
3298   hist->GetAxis(kGenSelectSpecies)->SetBinLabel(3, "#pi");
3299   hist->GetAxis(kGenSelectSpecies)->SetBinLabel(4, "p");
3300
3301   hist->GetAxis(kGenPt)->SetTitle("p_{T} (GeV/c)");
3302
3303   hist->GetAxis(kGenDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
3304
3305   hist->GetAxis(kGenCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
3306
3307   if (fStoreAdditionalJetInformation) {
3308     hist->GetAxis(kGenJetPt)->SetTitle("p_{T}^{jet} (GeV/c)");
3309
3310     hist->GetAxis(kGenZ)->SetTitle("z = p_{T}^{track} / p_{T}^{jet}");
3311
3312     hist->GetAxis(kGenXi)->SetTitle("#xi = ln(p_{T}^{jet} / p_{T}^{track})");
3313   }
3314
3315   hist->GetAxis(GetIndexOfChargeAxisGen())->SetTitle("Charge (e_{0})");
3316
3317   hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetTitle("TOF PID Info");
3318   // Offset is (TMath::Abs(kNoTOFinfo) + 1), such that bin 1 (= first label) corresponds to kNoTOFinfo (< 0)
3319   hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kNoTOFinfo + (TMath::Abs(kNoTOFinfo) + 1), "No TOF Info");
3320   hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kNoTOFpid + (TMath::Abs(kNoTOFinfo) + 1), "No TOF PID");
3321   hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kTOFpion + (TMath::Abs(kNoTOFinfo) + 1), "#pi");
3322   hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kTOFkaon + (TMath::Abs(kNoTOFinfo) + 1), "K");
3323   hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kTOFproton + (TMath::Abs(kNoTOFinfo) + 1), "p");
3324 }
3325
3326
3327 //________________________________________________________________________
3328 void AliAnalysisTaskPID::SetUpGenYieldHist(THnSparse* hist, Double_t* binsPt, Double_t* binsCent, Double_t* binsJetPt) const
3329 {
3330   // Sets bin limits for axes which are not standard binned and the axes titles.
3331
3332   hist->SetBinEdges(kGenYieldPt, binsPt);
3333   hist->SetBinEdges(kGenYieldCentrality, binsCent);
3334   if (fStoreAdditionalJetInformation)
3335     hist->SetBinEdges(kGenYieldJetPt, binsJetPt);
3336
3337   for (Int_t i = 0; i < 5; i++)
3338     hist->GetAxis(kGenYieldMCID)->SetBinLabel(i + 1, AliPID::ParticleLatexName(i));
3339
3340   // Set axes titles
3341   hist->GetAxis(kGenYieldMCID)->SetTitle("MC PID");
3342   hist->GetAxis(kGenYieldPt)->SetTitle("p_{T}^{gen} (GeV/c)");
3343   hist->GetAxis(kGenYieldCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
3344
3345   if (fStoreAdditionalJetInformation) {
3346     hist->GetAxis(kGenYieldJetPt)->SetTitle("p_{T}^{jet, gen} (GeV/c)");
3347
3348     hist->GetAxis(kGenYieldZ)->SetTitle("z = p_{T}^{track} / p_{T}^{jet}");
3349
3350     hist->GetAxis(kGenYieldXi)->SetTitle("#xi = ln(p_{T}^{jet} / p_{T}^{track})");
3351   }
3352
3353   hist->GetAxis(GetIndexOfChargeAxisGenYield())->SetTitle("Charge (e_{0})");
3354 }
3355
3356
3357 //________________________________________________________________________
3358 void AliAnalysisTaskPID::SetUpHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const
3359 {
3360   // Sets bin limits for axes which are not standard binned and the axes titles.
3361
3362   hist->SetBinEdges(kDataPt, binsPt);
3363   hist->SetBinEdges(kDataDeltaPrimeSpecies, binsDeltaPrime);
3364   hist->SetBinEdges(kDataCentrality, binsCent);
3365
3366   if (fStoreAdditionalJetInformation)
3367     hist->SetBinEdges(kDataJetPt, binsJetPt);
3368
3369   // Set axes titles
3370   hist->GetAxis(kDataMCID)->SetTitle("MC PID");
3371   hist->GetAxis(kDataMCID)->SetBinLabel(1, "e");
3372   hist->GetAxis(kDataMCID)->SetBinLabel(2, "K");
3373   hist->GetAxis(kDataMCID)->SetBinLabel(3, "#mu");
3374   hist->GetAxis(kDataMCID)->SetBinLabel(4, "#pi");
3375   hist->GetAxis(kDataMCID)->SetBinLabel(5, "p");
3376
3377   hist->GetAxis(kDataSelectSpecies)->SetTitle("Select Species");
3378   hist->GetAxis(kDataSelectSpecies)->SetBinLabel(1, "e");
3379   hist->GetAxis(kDataSelectSpecies)->SetBinLabel(2, "K");
3380   hist->GetAxis(kDataSelectSpecies)->SetBinLabel(3, "#pi");
3381   hist->GetAxis(kDataSelectSpecies)->SetBinLabel(4, "p");
3382
3383   hist->GetAxis(kDataPt)->SetTitle("p_{T} (GeV/c)");
3384
3385   hist->GetAxis(kDataDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
3386
3387   hist->GetAxis(kDataCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
3388
3389   if (fStoreAdditionalJetInformation) {
3390     hist->GetAxis(kDataJetPt)->SetTitle("p_{T}^{jet} (GeV/c)");
3391
3392     hist->GetAxis(kDataZ)->SetTitle("z = p_{T}^{track} / p_{T}^{jet}");
3393
3394     hist->GetAxis(kDataXi)->SetTitle("#xi = ln(p_{T}^{jet} / p_{T}^{track})");
3395   }
3396
3397   hist->GetAxis(GetIndexOfChargeAxisData())->SetTitle("Charge (e_{0})");
3398
3399   hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetTitle("TOF PID Info");
3400   // Offset is (TMath::Abs(kNoTOFinfo) + 1), such that bin 1 (= first label) corresponds to kNoTOFinfo (< 0)
3401   hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kNoTOFinfo + (TMath::Abs(kNoTOFinfo) + 1), "No TOF Info");
3402   hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kNoTOFpid + (TMath::Abs(kNoTOFinfo) + 1), "No TOF PID");
3403   hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kTOFpion + (TMath::Abs(kNoTOFinfo) + 1), "#pi");
3404   hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kTOFkaon + (TMath::Abs(kNoTOFinfo) + 1), "K");
3405   hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kTOFproton + (TMath::Abs(kNoTOFinfo) + 1), "p");
3406 }
3407
3408
3409 //________________________________________________________________________
3410 void AliAnalysisTaskPID::SetUpPtResHist(THnSparse* hist, Double_t* binsPt, Double_t* binsJetPt, Double_t* binsCent) const
3411 {
3412   // Sets bin limits for axes which are not standard binned and the axes titles.
3413
3414   hist->SetBinEdges(kPtResJetPt, binsJetPt);
3415   hist->SetBinEdges(kPtResGenPt, binsPt);
3416   hist->SetBinEdges(kPtResRecPt, binsPt);
3417   hist->SetBinEdges(kPtResCentrality, binsCent);
3418
3419   // Set axes titles
3420   hist->GetAxis(kPtResJetPt)->SetTitle("p_{T}^{jet, rec} (GeV/c)");
3421   hist->GetAxis(kPtResGenPt)->SetTitle("p_{T}^{gen} (GeV/c)");
3422   hist->GetAxis(kPtResRecPt)->SetTitle("p_{T}^{rec} (GeV/c)");
3423
3424   hist->GetAxis(kPtResCharge)->SetTitle("Charge (e_{0})");
3425   hist->GetAxis(kPtResCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
3426 }