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