PWGJE/UserTasks/AliAnalysisTaskPID.cxx

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