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