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     Double_t dEdxTPC = fPIDResponse->IsTunedOnData() ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
 909     if (dEdxTPC <= 0)
 910       continue;
 911
 912     if(fTrackFilter && !fTrackFilter->IsSelected(track))
 913       continue;
 914
 915     if (GetUseTPCCutMIGeo()) {
 916       if (!TPCCutMIGeo(track, fEvent))
 917         continue;
 918     }
 919     else if (GetUseTPCnclCut()) {
 920       if (!TPCnclCut(track))
 921         continue;
 922     }
 923
 924     if(fUsePhiCut) {
 925       if (!PhiPrimeCut(track, magField))
 926         continue; // reject track
 927     }
 928
 929     Int_t pdg =  0; // = 0 indicates data for the moment
 930     AliMCParticle* mcTrack = 0x0;
 931     Int_t mcID = AliPID::kUnknown;
 932     Int_t label = 0;
 933
 934     if (fMC) {
 935       label = track->GetLabel();
 936
 937       //if (label < 0)
 938       //  continue;
 939
 940       mcTrack = dynamic_cast<AliMCParticle*>(fMC->GetTrack(TMath::Abs(label)));
 941       if (!mcTrack) {
 942         Printf("ERROR: Could not retrieve mcTrack with label %d for track %d", label, iTracks);
 943         continue;
 944       }
 945
 946       pdg = mcTrack->PdgCode();
 947       mcID = PDGtoMCID(mcTrack->PdgCode());
 948
 949       if (fDoEfficiency) {
 950         // For efficiency: Reconstructed track has survived all cuts on the detector level (excluding acceptance)
 951         // and has an associated MC track which is a physical primary and was generated inside the acceptance
 952         if (fMC->IsPhysicalPrimary(TMath::Abs(label)) &&
 953             (TMath::Abs(mcTrack->Eta()) >= fEtaAbsCutLow && TMath::Abs(mcTrack->Eta()) <= fEtaAbsCutUp)) {
 954
 955           // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
 956           Double_t value[kEffNumAxes] = { mcID, mcTrack->Pt(), mcTrack->Eta(), mcTrack->Charge() / 3., centralityPercentile,
 957                                           -1, -1, -1 };
 958           fContainerEff->Fill(value, kStepRecWithGenCuts);
 959
 960           Double_t valueMeas[kEffNumAxes] = { mcID, track->Pt(), track->Eta(), track->Charge(), centralityPercentile,
 961                                               -1, -1, -1 };
 962           fContainerEff->Fill(valueMeas, kStepRecWithGenCutsMeasuredObs);
 963         }
 964       }
 965     }
 966
 967     // Only process tracks inside the desired eta window
 968     if (TMath::Abs(track->Eta()) < fEtaAbsCutLow || TMath::Abs(track->Eta()) > fEtaAbsCutUp)  continue;
 969
 970     if (fDoPID)
 971       ProcessTrack(track, pdg, centralityPercentile, -1); // No jet information in this case -> Set jet pT to -1
 972
 973     if (fDoPtResolution) {
 974       if (mcTrack && fMC->IsPhysicalPrimary(TMath::Abs(label))) {
 975         // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
 976         Double_t valuePtRes[kPtResNumAxes] = { -1, mcTrack->Pt(), track->Pt(), mcTrack->Charge() / 3., centralityPercentile };
 977         fPtResolution[mcID]->Fill(valuePtRes);
 978       }
 979     }
 980
 981     if (fDoEfficiency) {
 982       if (mcTrack) {
 983         Double_t valueRecAllCuts[kEffNumAxes] = { mcID, track->Pt(), track->Eta(), track->Charge(), centralityPercentile,
 984                                                   -1, -1, -1 };
 985         fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObs);
 986
 987         Double_t weight = IsSecondaryWithStrangeMotherMC(fMC, TMath::Abs(label)) ?
 988                                                                 GetMCStrangenessFactorCMS(fMC, mcTrack) : 1.0;
 989         fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObsStrangenessScaled, weight);
 990
 991         // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
 992         Double_t valueGenAllCuts[kEffNumAxes] = { mcID, mcTrack->Pt(), mcTrack->Eta(), mcTrack->Charge() / 3.,
 993                                                   centralityPercentile, -1, -1, -1 };
 994         if (fMC->IsPhysicalPrimary(TMath::Abs(label))) {
 995           fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObsPrimaries);
 996           fContainerEff->Fill(valueGenAllCuts, kStepRecWithRecCutsPrimaries);
 997         }
 998       }
 999     }
1000   } //track loop
1001
1002   IncrementEventsProcessed(centralityPercentile);
1003
1004   if(fDebug > 2)
1005     printf("File: %s, Line: %d: UserExec -> Processing done\n", (char*)__FILE__, __LINE__);
1006
1007   PostOutputData();
1008
1009   if(fDebug > 2)
1010     printf("File: %s, Line: %d: UserExec -> Done\n", (char*)__FILE__, __LINE__);
1011 }
1012
1013 //________________________________________________________________________
1014 void AliAnalysisTaskPID::Terminate(const Option_t *)
1015 {
1016   // Draw result to the screen
1017   // Called once at the end of the query
1018 }
1019
1020
1021 //_____________________________________________________________________________
1022 void AliAnalysisTaskPID::CheckDoAnyStematicStudiesOnTheExpectedSignal()
1023 {
1024   // Check whether at least one scale factor indicates the ussage of systematic studies
1025   // and set internal flag accordingly.
1026
1027   fDoAnySystematicStudiesOnTheExpectedSignal = kFALSE;
1028
1029   if (TMath::Abs(fSystematicScalingSplines - 1.0) > fgkEpsilon) {
1030     fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1031     return;
1032   }
1033
1034   if ((TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon) ||
1035       (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon)) {
1036     fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1037     return;
1038   }
1039
1040   if (TMath::Abs(fSystematicScalingEtaSigmaPara - 1.0) > fgkEpsilon) {
1041     fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1042     return;
1043   }
1044
1045   if (TMath::Abs(fSystematicScalingMultCorrection - 1.0) > fgkEpsilon) {
1046     fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1047     return;
1048   }
1049 }
1050
1051
1052 //_____________________________________________________________________________
1053 Int_t AliAnalysisTaskPID::PDGtoMCID(Int_t pdg)
1054 {
1055   // Returns the corresponding AliPID index to the given pdg code.
1056   // Returns AliPID::kUnkown if pdg belongs to a not considered species.
1057
1058   Int_t absPDGcode = TMath::Abs(pdg);
1059   if (absPDGcode == 211) {//Pion
1060     return AliPID::kPion;
1061   }
1062   else if (absPDGcode == 321) {//Kaon
1063     return AliPID::kKaon;
1064   }
1065   else if (absPDGcode == 2212) {//Proton
1066     return AliPID::kProton;
1067   }
1068   else if (absPDGcode == 11) {//Electron
1069     return AliPID::kElectron;
1070   }
1071   else if (absPDGcode == 13) {//Muon
1072     return AliPID::kMuon;
1073   }
1074
1075   return AliPID::kUnknown;
1076 }
1077
1078
1079 //_____________________________________________________________________________
1080 void AliAnalysisTaskPID::GetJetTrackObservables(Double_t trackPt, Double_t jetPt, Double_t& z, Double_t& xi)
1081 {
1082   // Uses trackPt and jetPt to obtain z and xi.
1083
1084   z = (jetPt > 0 && trackPt >= 0) ? (trackPt / jetPt) : -1;
1085   xi = (z > 0) ? TMath::Log(1. / z) : -1;
1086
1087   if(trackPt > (1. - 1e-06) * jetPt && trackPt < (1. + 1e-06) * jetPt) { // case z=1 : move entry to last histo bin <1
1088     z  = 1. - 1e-06;
1089     xi = 1e-06;
1090   }
1091 }
1092
1093
1094 //_____________________________________________________________________________
1095 void AliAnalysisTaskPID::CleanupParticleFractionHistos()
1096 {
1097   // Delete histos with particle fractions
1098
1099   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1100     delete fFractionHists[i];
1101     fFractionHists[i] = 0x0;
1102
1103     delete fFractionSysErrorHists[i];
1104     fFractionSysErrorHists[i] = 0x0;
1105   }
1106 }
1107
1108
1109 //_____________________________________________________________________________
1110 Double_t AliAnalysisTaskPID::ConvolutedGaus(const Double_t* xx, const Double_t* par) const
1111 {
1112   // Convolutes gauss with an exponential tail which describes dEdx-response better than pure gaussian
1113
1114   const Double_t mean = par[0];
1115   const Double_t sigma = par[1];
1116   const Double_t lambda = par[2];
1117
1118   if(fDebug > 5)
1119     printf("File: %s, Line: %d: ConvolutedGaus: mean %e, sigma %e, lambda %e\n", (char*)__FILE__, __LINE__, mean, sigma, lambda);
1120
1121   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);
1122 }
1123
1124
1125 //_____________________________________________________________________________
1126 inline Double_t AliAnalysisTaskPID::FastGaus(Double_t x, Double_t mean, Double_t sigma) const
1127 {
1128   // Calculate an unnormalised gaussian function with mean and sigma.
1129
1130   if (sigma < fgkEpsilon)
1131     return 1.e30;
1132
1133   const Double_t arg = (x - mean) / sigma;
1134   return exp(-0.5 * arg * arg);
1135 }
1136
1137
1138 //_____________________________________________________________________________
1139 inline Double_t AliAnalysisTaskPID::FastNormalisedGaus(Double_t x, Double_t mean, Double_t sigma) const
1140 {
1141   // Calculate a normalised (divided by sqrt(2*Pi)*sigma) gaussian function with mean and sigma.
1142
1143   if (sigma < fgkEpsilon)
1144     return 1.e30;
1145
1146   const Double_t arg = (x - mean) / sigma;
1147   const Double_t res = exp(-0.5 * arg * arg);
1148   return res / (2.50662827463100024 * sigma); //sqrt(2*Pi)=2.50662827463100024
1149 }
1150
1151
1152 //_____________________________________________________________________________
1153 Int_t AliAnalysisTaskPID::FindBinWithinRange(TAxis* axis, Double_t value) const
1154 {
1155   // Find the corresponding bin of the axis. Values outside the range (also under and overflow) will be set to the first/last
1156   // available bin
1157
1158   Int_t bin = axis->FindFixBin(value);
1159
1160   if (bin <= 0)
1161     bin = 1;
1162   if (bin > axis->GetNbins())
1163     bin = axis->GetNbins();
1164
1165   return bin;
1166 }
1167
1168
1169 //_____________________________________________________________________________
1170 Int_t AliAnalysisTaskPID::FindFirstBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yBin,
1171                                                       Int_t zBin) const
1172 {
1173   // Kind of projects a TH3 to 1 bin combination in y and z
1174   // and looks for the first x bin above a threshold for this projection.
1175   // If no such bin is found, -1 is returned.
1176
1177   if (!hist)
1178     return -1;
1179
1180   Int_t nBinsX = hist->GetNbinsX();
1181   for (Int_t xBin = 1; xBin <= nBinsX; xBin++) {
1182     if (hist->GetBinContent(xBin, yBin, zBin) > threshold)
1183       return xBin;
1184   }
1185
1186   return -1;
1187 }
1188
1189
1190 //_____________________________________________________________________________
1191 Int_t AliAnalysisTaskPID::FindLastBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yBin,
1192                                                      Int_t zBin) const
1193 {
1194   // Kind of projects a TH3 to 1 bin combination in y and z
1195   // and looks for the last x bin above a threshold for this projection.
1196   // If no such bin is found, -1 is returned.
1197
1198   if (!hist)
1199     return -1;
1200
1201   Int_t nBinsX = hist->GetNbinsX();
1202   for (Int_t xBin = nBinsX; xBin >= 1; xBin--) {
1203     if (hist->GetBinContent(xBin, yBin, zBin) > threshold)
1204       return xBin;
1205   }
1206
1207   return -1;
1208 }
1209
1210
1211 //_____________________________________________________________________________
1212 Bool_t AliAnalysisTaskPID::GetParticleFraction(Double_t trackPt, Double_t jetPt, Double_t centralityPercentile,
1213                                                AliPID::EParticleType species,
1214                                                Double_t& fraction, Double_t& fractionErrorStat, Double_t& fractionErrorSys) const
1215 {
1216   // Computes the particle fraction for the corresponding species for the given trackPt, jetPt and centrality.
1217   // Use jetPt = -1 for inclusive spectra and centralityPercentile = -1 for pp.
1218   // On success (return value kTRUE), fraction contains the particle fraction, fractionErrorStat(Sys) the sigma of its
1219   // statistical (systematic) error
1220
1221   fraction = -999.;
1222   fractionErrorStat = 999.;
1223   fractionErrorSys = 999.;
1224
1225   if (species > AliPID::kProton || species < AliPID::kElectron) {
1226     AliError(Form("Only fractions for species index %d to %d availabe, but not for the requested one: %d", 0, AliPID::kProton, species));
1227     return kFALSE;
1228   }
1229
1230   if (!fFractionHists[species]) {
1231     AliError(Form("Histo with particle fractions for species %d not loaded!", species));
1232
1233     return kFALSE;
1234   }
1235
1236   Int_t jetPtBin = FindBinWithinRange(fFractionHists[species]->GetYaxis(), jetPt);
1237   Int_t centBin  = FindBinWithinRange(fFractionHists[species]->GetZaxis(), centralityPercentile);
1238
1239   // The following interpolation takes the bin content of the first/last available bin,
1240   // if requested point lies beyond bin center of first/last bin.
1241   // The interpolation is only done for the x-axis (track pT), i.e. jetPtBin and centBin are fix,
1242   // because the analysis will anyhow be run in bins of jetPt and centrality and
1243   // it is not desired to correlate different jetPt bins via interpolation.
1244
1245   // The same procedure is used for the error of the fraction
1246   TAxis* xAxis = fFractionHists[species]->GetXaxis();
1247
1248   // No interpolation to values beyond the centers of the first/last bins (we don't know exactly where the spectra start or stop,
1249   // thus, search for the first and last bin above 0.0 to constrain the range
1250   Int_t firstBin = TMath::Max(1, FindFirstBinAboveIn3dSubset(fFractionHists[species], 0.0, jetPtBin, centBin));
1251   Int_t lastBin  = TMath::Min(fFractionHists[species]->GetNbinsX(),
1252                               FindLastBinAboveIn3dSubset(fFractionHists[species], 0.0, jetPtBin, centBin));
1253
1254   if (trackPt <= xAxis->GetBinCenter(firstBin)) {
1255     fraction = fFractionHists[species]->GetBinContent(firstBin, jetPtBin, centBin);
1256     fractionErrorStat = fFractionHists[species]->GetBinError(firstBin, jetPtBin, centBin);
1257     fractionErrorSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(firstBin, jetPtBin, centBin) : 0.;
1258   }
1259   else if (trackPt >= xAxis->GetBinCenter(lastBin)) {
1260     fraction = fFractionHists[species]->GetBinContent(lastBin, jetPtBin, centBin);
1261     fractionErrorStat = fFractionHists[species]->GetBinError(lastBin, jetPtBin, centBin);
1262     fractionErrorSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(lastBin, jetPtBin, centBin) : 0.;
1263   }
1264   else {
1265     Double_t x0 = 0., x1 = 0., y0 = 0., y1 = 0.;
1266     Double_t y0errStat = 0., y1errStat = 0., y0errSys = 0., y1errSys = 0.;
1267     Int_t trackPtBin = xAxis->FindBin(trackPt);
1268
1269     // Linear interpolation between nearest neighbours in trackPt
1270     if (trackPt <= xAxis->GetBinCenter(trackPtBin)) {
1271         y0 = fFractionHists[species]->GetBinContent(trackPtBin - 1, jetPtBin, centBin);
1272         y0errStat = fFractionHists[species]->GetBinError(trackPtBin - 1, jetPtBin, centBin);
1273         y0errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin - 1, jetPtBin, centBin)
1274                                                    : 0.;
1275         x0 = xAxis->GetBinCenter(trackPtBin - 1);
1276         y1 = fFractionHists[species]->GetBinContent(trackPtBin, jetPtBin, centBin);
1277         y1errStat = fFractionHists[species]->GetBinError(trackPtBin, jetPtBin, centBin);
1278         y1errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin, jetPtBin, centBin)
1279                                                    : 0.;
1280         x1 = xAxis->GetBinCenter(trackPtBin);
1281     }
1282     else {
1283         y0 = fFractionHists[species]->GetBinContent(trackPtBin, jetPtBin, centBin);
1284         y0errStat = fFractionHists[species]->GetBinError(trackPtBin, jetPtBin, centBin);
1285         y0errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin, jetPtBin, centBin)
1286                                                    : 0.;
1287         x0 = xAxis->GetBinCenter(trackPtBin);
1288         y1 = fFractionHists[species]->GetBinContent(trackPtBin + 1, jetPtBin, centBin);
1289         y1errStat = fFractionHists[species]->GetBinError(trackPtBin + 1, jetPtBin, centBin);
1290         y1errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin + 1, jetPtBin, centBin)
1291                                                    : 0.;
1292         x1 = xAxis->GetBinCenter(trackPtBin + 1);
1293     }
1294
1295     // Per construction: x0 < trackPt < x1
1296     fraction = y0 + (trackPt - x0) * ((y1 - y0) / (x1 - x0));
1297     fractionErrorStat = y0errStat + (trackPt - x0) * ((y1errStat - y0errStat) / (x1 - x0));
1298     fractionErrorSys = fFractionSysErrorHists[species] ? (y0errSys + (trackPt - x0) * ((y1errSys - y0errSys) / (x1 - x0))) : 0.;
1299   }
1300
1301   return kTRUE;
1302 }
1303
1304
1305 //_____________________________________________________________________________
1306 Bool_t AliAnalysisTaskPID::GetParticleFractions(Double_t trackPt, Double_t jetPt, Double_t centralityPercentile,
1307                                                 Double_t* prob, Int_t smearSpeciesByError,
1308                                                 Int_t takeIntoAccountSpeciesSysError, Bool_t uniformSystematicError) const
1309 {
1310   // Fills the particle fractions for the given trackPt, jetPt and centrality into "prob".
1311   // Use jetPt = -1 for inclusive spectra and centralityPercentile = -1 for pp.
1312   // If smearSpeciesByError is >= 0 && < AliPID::kSPECIES, the returned fractions will be a random number distributed
1313   // with a gauss with mean being the corresponding particle fraction and sigma it's error for the considered species
1314   // "smearSpeciesByError".
1315   // Note that in this case the fractions for all species will NOT sum up to 1!
1316   // Thus, all other species fractions will be re-scaled weighted with their corresponding statistical error.
1317   // A similar procedure is used for "takeIntoAccountSpeciesSysError":  The systematic error of the corresponding species
1318   // is used to generate a random number with uniform distribution in [mean - sysError, mean + sysError] for the new mean
1319   // (in cace of uniformSystematicError = kTRUE, otherwise it will be a gaus(mean, sysError)),
1320   // then the other species will be re-scaled according to their systematic errors.
1321   // First, the systematic error uncertainty procedure will be performed (that is including re-scaling), then the statistical
1322   // uncertainty procedure.
1323   // On success, kTRUE is returned.
1324
1325   if (!prob || smearSpeciesByError >= AliPID::kSPECIES || takeIntoAccountSpeciesSysError >= AliPID::kSPECIES)
1326     return kFALSE;
1327
1328   Double_t probTemp[AliPID::kSPECIES];
1329   Double_t probErrorStat[AliPID::kSPECIES];
1330   Double_t probErrorSys[AliPID::kSPECIES];
1331
1332   Bool_t success = kTRUE;
1333   success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kElectron,
1334                                            probTemp[AliPID::kElectron], probErrorStat[AliPID::kElectron],
1335                                            probErrorSys[AliPID::kElectron]);
1336   success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kMuon,
1337                                            probTemp[AliPID::kMuon], probErrorStat[AliPID::kMuon], probErrorSys[AliPID::kMuon]);
1338   success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kPion,
1339                                            probTemp[AliPID::kPion], probErrorStat[AliPID::kPion], probErrorSys[AliPID::kPion]);
1340   success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kKaon,
1341                                            probTemp[AliPID::kKaon], probErrorStat[AliPID::kKaon], probErrorSys[AliPID::kKaon]);
1342   success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kProton,
1343                                            probTemp[AliPID::kProton], probErrorStat[AliPID::kProton], probErrorSys[AliPID::kProton]);
1344
1345   if (!success)
1346     return kFALSE;
1347
1348   // If desired, take into account the systematic error of the corresponding species and re-generate probTemp accordingly
1349   if (takeIntoAccountSpeciesSysError >= 0) {
1350     // Generate random fraction of the considered species "smearSpeciesByError" according to mean and sigma
1351     Double_t generatedFraction = uniformSystematicError
1352                                    ? fRandom->Rndm() * 2. * probErrorSys[takeIntoAccountSpeciesSysError]
1353                                      - probErrorSys[takeIntoAccountSpeciesSysError]
1354                                      + probTemp[takeIntoAccountSpeciesSysError]
1355                                    : fRandom->Gaus(probTemp[takeIntoAccountSpeciesSysError],
1356                                                    probErrorSys[takeIntoAccountSpeciesSysError]);
1357
1358     // Catch cases with invalid fraction (can happen for large errors), i.e. fraction < 0 or > 1
1359     if (generatedFraction < 0.)
1360       generatedFraction = 0.;
1361     else if (generatedFraction > 1.)
1362       generatedFraction = 1.;
1363
1364     // Calculate difference from original fraction (original fractions sum up to 1!)
1365     Double_t deltaFraction = generatedFraction - probTemp[takeIntoAccountSpeciesSysError];
1366
1367     // Fractions must (including errors) lie inside [0,1] -> Adapt weights accordingly by setting the errors
1368     if (deltaFraction > 0) {
1369       // Some part will be SUBTRACTED from the other fractions
1370       for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1371         if (probTemp[i] - probErrorSys[i] < 0)
1372           probErrorSys[i] = probTemp[i];
1373       }
1374     }
1375     else {
1376       // Some part will be ADDED to the other fractions
1377       for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1378         if (probTemp[i] + probErrorSys[i] > 1)
1379           probErrorSys[i] = 1. - probTemp[i];
1380       }
1381     }
1382
1383     // Compute summed weight of all fractions except for the considered one
1384     Double_t summedWeight = 0.;
1385     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1386       if (i != takeIntoAccountSpeciesSysError)
1387         summedWeight += probErrorSys[i];
1388     }
1389
1390     // Compute the weight for the other species
1391     /*
1392     if (summedWeight <= 1e-13) {
1393       // If this happens for some reason (it should not!), just assume flat weight
1394       printf("Error: summedWeight (sys error) ~ 0 for trackPt %f, jetPt %f, centralityPercentile %f. Setting flat weight!\n",
1395              trackPt, jetPt, centralityPercentile);
1396     }*/
1397
1398     Double_t weight[AliPID::kSPECIES];
1399     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1400       if (i != takeIntoAccountSpeciesSysError) {
1401         if (summedWeight > 1e-13)
1402           weight[i] = probErrorSys[i] / summedWeight;
1403         else
1404           weight[i] = probErrorSys[i] / (AliPID::kSPECIES - 1);
1405       }
1406     }
1407
1408     // For the final generated fractions, set the generated value for the considered species
1409     // and the generated value minus delta times statistical weight
1410     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1411       if (i != takeIntoAccountSpeciesSysError)
1412         probTemp[i] = probTemp[i] - weight[i] * deltaFraction;
1413       else
1414         probTemp[i] = generatedFraction;
1415     }
1416   }
1417
1418   // Using the values of probTemp (either the original ones or those after taking into account the systematic error),
1419   // calculate the final fractions - if the statistical error is to be taken into account, smear the corresponding
1420   // fraction. If not, just write probTemp to the final result array.
1421   if (smearSpeciesByError >= 0) {
1422     // Generate random fraction of the considered species "smearSpeciesByError" according to mean and sigma
1423     Double_t generatedFraction = fRandom->Gaus(probTemp[smearSpeciesByError], probErrorStat[smearSpeciesByError]);
1424
1425     // Catch cases with invalid fraction (can happen for large errors), i.e. fraction < 0 or > 1
1426     if (generatedFraction < 0.)
1427       generatedFraction = 0.;
1428     else if (generatedFraction > 1.)
1429       generatedFraction = 1.;
1430
1431     // Calculate difference from original fraction (original fractions sum up to 1!)
1432     Double_t deltaFraction = generatedFraction - probTemp[smearSpeciesByError];
1433
1434     // Fractions must (including errors) lie inside [0,1] -> Adapt weights accordingly by setting the errors
1435     if (deltaFraction > 0) {
1436       // Some part will be SUBTRACTED from the other fractions
1437       for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1438         if (probTemp[i] - probErrorStat[i] < 0)
1439           probErrorStat[i] = probTemp[i];
1440       }
1441     }
1442     else {
1443       // Some part will be ADDED to the other fractions
1444       for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1445         if (probTemp[i] + probErrorStat[i] > 1)
1446           probErrorStat[i] = 1. - probTemp[i];
1447       }
1448     }
1449
1450     // Compute summed weight of all fractions except for the considered one
1451     Double_t summedWeight = 0.;
1452     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1453       if (i != smearSpeciesByError)
1454         summedWeight += probErrorStat[i];
1455     }
1456
1457     // Compute the weight for the other species
1458     /*
1459     if (summedWeight <= 1e-13) {
1460       // If this happens for some reason (it should not!), just assume flat weight
1461       printf("Error: summedWeight (stat error) ~ 0 for trackPt %f, jetPt %f, centralityPercentile %f. Setting flat weight!\n",
1462              trackPt, jetPt, centralityPercentile);
1463     }*/
1464
1465     Double_t weight[AliPID::kSPECIES];
1466     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1467       if (i != smearSpeciesByError) {
1468         if (summedWeight > 1e-13)
1469           weight[i] = probErrorStat[i] / summedWeight;
1470         else
1471           weight[i] = probErrorStat[i] / (AliPID::kSPECIES - 1);
1472       }
1473     }
1474
1475     // For the final generated fractions, set the generated value for the considered species
1476     // and the generated value minus delta times statistical weight
1477     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1478       if (i != smearSpeciesByError)
1479         prob[i] = probTemp[i] - weight[i] * deltaFraction;
1480       else
1481         prob[i] = generatedFraction;
1482     }
1483   }
1484   else {
1485     // Just take the generated values
1486     for (Int_t i = 0; i < AliPID::kSPECIES; i++)
1487       prob[i] = probTemp[i];
1488   }
1489
1490
1491   // Should already be normalised, but make sure that it really is:
1492   Double_t probSum = 0.;
1493   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1494     probSum += prob[i];
1495   }
1496
1497   if (probSum <= 0)
1498     return kFALSE;
1499
1500   if (TMath::Abs(probSum - 1.0) > 1e-4) {
1501     printf("Warning: Re-normalising sum of fractions: Sum is %e\n", probSum);
1502     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1503       prob[i] /= probSum;
1504     }
1505   }
1506
1507   return kTRUE;
1508 }
1509
1510
1511 //_____________________________________________________________________________
1512 const TH3D* AliAnalysisTaskPID::GetParticleFractionHisto(Int_t species, Bool_t sysError) const
1513 {
1514   if (species < AliPID::kElectron || species > AliPID::kProton)
1515     return 0x0;
1516
1517   return sysError ? fFractionSysErrorHists[species] : fFractionHists[species];
1518 }
1519
1520
1521 //_____________________________________________________________________________
1522 Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(Int_t motherPDG, Double_t motherGenPt)
1523 {
1524   // Strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
1525   // -> Based on function in PWGJE/AliAnalysisTaskFragmentationFunction, which uses
1526   // the following data from CMS pp @ 7 TeV inclusive (JHEP 05 (2011) 064)
1527
1528   Double_t fac = 1;
1529
1530   const Int_t absMotherPDG = TMath::Abs(motherPDG);
1531
1532   if (absMotherPDG == 310 || absMotherPDG == 321) { // K0s / K+ / K-
1533     if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.768049;
1534     else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.732933;
1535     else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.650298;
1536     else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.571332;
1537     else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.518734;
1538     else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.492543;
1539     else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.482704;
1540     else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.488056;
1541     else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.488861;
1542     else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.492862;
1543     else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.504332;
1544     else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.501858;
1545     else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.512970;
1546     else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.524131;
1547     else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.539130;
1548     else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.554101;
1549     else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.560348;
1550     else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.568869;
1551     else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.583310;
1552     else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.604818;
1553     else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.632630;
1554     else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.710070;
1555     else if(6.00 <= motherGenPt && motherGenPt < 8.00) fac = 0.736365;
1556     else if(8.00 <= motherGenPt && motherGenPt < 10.00) fac = 0.835865;
1557   }
1558
1559   if (absMotherPDG == 3122) { // Lambda
1560     if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.645162;
1561     else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.627431;
1562     else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.457136;
1563     else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.384369;
1564     else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.330597;
1565     else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.309571;
1566     else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.293620;
1567     else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.283709;
1568     else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.282047;
1569     else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.277261;
1570     else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.275772;
1571     else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.280726;
1572     else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.288540;
1573     else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.288315;
1574     else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.296619;
1575     else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.302993;
1576     else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.338121;
1577     else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.349800;
1578     else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.356802;
1579     else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.391202;
1580     else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.422573;
1581     else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.573815;
1582     else if(6.00 <= motherGenPt && motherGenPt < 8.00) fac = 0.786984;
1583     else if(8.00 <= motherGenPt && motherGenPt < 10.00) fac = 1.020021;
1584   }
1585
1586   if (absMotherPDG == 3312 || absMotherPDG == 3322) { // xi
1587     if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.666620;
1588     else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.575908;
1589     else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.433198;
1590     else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.340901;
1591     else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.290896;
1592     else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.236074;
1593     else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.218681;
1594     else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.207763;
1595     else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.222848;
1596     else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.208806;
1597     else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.197275;
1598     else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.183645;
1599     else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.188788;
1600     else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.188282;
1601     else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.207442;
1602     else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.240388;
1603     else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.241916;
1604     else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.208276;
1605     else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.234550;
1606     else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.251689;
1607     else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.310204;
1608     else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.343492;
1609   }
1610
1611   const Double_t weight = 1. / fac;
1612
1613   return weight;
1614 }
1615
1616
1617 //_____________________________________________________________________________
1618 Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(AliMCEvent* mcEvent, AliMCParticle* daughter)
1619 {
1620   // Strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
1621   // -> Based on function in PWGJE/AliAnalysisTaskFragmentationFunction
1622
1623   if (!mcEvent)
1624     return 1.;
1625
1626   AliMCParticle* currentMother   = daughter;
1627   AliMCParticle* currentDaughter = daughter;
1628
1629
1630   // find first primary mother K0s, Lambda or Xi
1631   while(1) {
1632     Int_t daughterPDG = currentDaughter->PdgCode();
1633
1634     Int_t motherLabel = currentDaughter->GetMother();
1635     if(motherLabel >= mcEvent->GetNumberOfTracks()){ // protection
1636       currentMother = currentDaughter;
1637       break;
1638     }
1639
1640     currentMother = (AliMCParticle*)mcEvent->GetTrack(motherLabel);
1641
1642     if (!currentMother) {
1643       currentMother = currentDaughter;
1644       break;
1645     }
1646
1647     Int_t motherPDG = currentMother->PdgCode();
1648
1649     // phys. primary found ?
1650     if (mcEvent->IsPhysicalPrimary(motherLabel))
1651       break;
1652
1653     if (TMath::Abs(daughterPDG) == 321) {
1654       // K+/K- e.g. from phi (ref data not feeddown corrected)
1655       currentMother = currentDaughter;
1656       break;
1657     }
1658     if (TMath::Abs(motherPDG) == 310) {
1659       // K0s e.g. from phi (ref data not feeddown corrected)
1660       break;
1661     }
1662     if (TMath::Abs(motherPDG) == 3212 && TMath::Abs(daughterPDG) == 3122) {
1663       // Mother Sigma0, daughter Lambda (this case not included in feeddown corr.)
1664       currentMother = currentDaughter;
1665       break;
1666     }
1667
1668     currentDaughter = currentMother;
1669   }
1670
1671
1672   Int_t motherPDG = currentMother->PdgCode();
1673   Double_t motherGenPt = currentMother->Pt();
1674
1675   return GetMCStrangenessFactorCMS(motherPDG, motherGenPt);
1676 }
1677
1678
1679 // _________________________________________________________________________________
1680 Bool_t AliAnalysisTaskPID::IsSecondaryWithStrangeMotherMC(AliMCEvent* mcEvent, Int_t partLabel)
1681 {
1682   // Check whether particle is a secondary with strange mother, i.e. returns kTRUE if a strange mother is found
1683   // and the particle is NOT a physical primary. In all other cases kFALSE is returned
1684
1685   if (!mcEvent || partLabel < 0)
1686     return kFALSE;
1687
1688   AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(partLabel);
1689
1690   if (!part)
1691     return kFALSE;
1692
1693   if (mcEvent->IsPhysicalPrimary(partLabel))
1694     return kFALSE;
1695
1696   Int_t iMother = part->GetMother();
1697   if (iMother < 0)
1698     return kFALSE;
1699
1700
1701   AliMCParticle* partM = (AliMCParticle*)mcEvent->GetTrack(iMother);
1702   if (!partM)
1703    return kFALSE;
1704
1705   Int_t codeM = TMath::Abs(partM->PdgCode());
1706   Int_t mfl = Int_t(codeM / TMath::Power(10, Int_t(TMath::Log10(codeM))));
1707   if (mfl == 3 && codeM != 3) // codeM = 3 is for s quark
1708     return kTRUE;
1709
1710   return kFALSE;
1711 }
1712
1713
1714 //_____________________________________________________________________________
1715 Bool_t AliAnalysisTaskPID::SetParticleFractionHisto(const TH3D* hist, Int_t species, Bool_t sysError)
1716 {
1717   // Store a clone of hist (containing the particle fractions of the corresponding species with statistical error (sysError = kFALSE)
1718   // or systematic error (sysError = kTRUE), respectively), internally
1719
1720   if (species < AliPID::kElectron || species > AliPID::kProton) {
1721     AliError(Form("Only fractions for species index %d to %d can be set, but not for the requested one: %d", 0,
1722                   AliPID::kProton, species));
1723     return kFALSE;
1724   }
1725
1726   if (sysError) {
1727     delete fFractionSysErrorHists[species];
1728
1729     fFractionSysErrorHists[species] = new TH3D(*hist);
1730   }
1731   else {
1732     delete fFractionHists[species];
1733
1734     fFractionHists[species] = new TH3D(*hist);
1735   }
1736
1737   return kTRUE;
1738 }
1739
1740
1741 //_____________________________________________________________________________
1742 Bool_t AliAnalysisTaskPID::SetParticleFractionHistosFromFile(const TString filePathName, Bool_t sysError)
1743 {
1744   // Loads particle fractions for all species from the desired file and returns kTRUE on success.
1745   // The maps are assumed to be of Type TH3D, to sit in the main directory and to have names
1746   // Form("hFraction_%e", AliPID::ParticleName(i)) for sysError = kFALSE and
1747   // Form("hFractionSysError_%e", AliPID::ParticleName(i)) for sysError = kTRUE.
1748
1749   TFile* f = TFile::Open(filePathName.Data());
1750   if (!f)  {
1751     std::cout << "Failed to open file with particle fractions \"" << filePathName.Data() << "\"!" << std::endl;
1752     return kFALSE;
1753   }
1754
1755   TH3D* hist = 0x0;
1756   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1757     TString histName = Form("hFraction%s_%s", sysError ? "SysError" : "", AliPID::ParticleName(i));
1758     hist = dynamic_cast<TH3D*>(f->Get(histName.Data()));
1759     if (!hist) {
1760       std::cout << "Failed to load particle fractions for " << histName.Data() << "!";
1761       std::cout << std::endl << "Cleaning up particle fraction histos!" << std::endl;
1762       CleanupParticleFractionHistos();
1763       return kFALSE;
1764     }
1765
1766     if (!SetParticleFractionHisto(hist, i, sysError)) {
1767       std::cout << "Failed to load particle fractions for " << histName.Data() << "!";
1768       std::cout << std::endl << "Cleaning up particle fraction histos!" << std::endl;
1769       CleanupParticleFractionHistos();
1770       return kFALSE;
1771     }
1772   }
1773
1774   delete hist;
1775
1776   return kTRUE;
1777
1778 }
1779
1780
1781 //_____________________________________________________________________________
1782 Int_t AliAnalysisTaskPID::GetRandomParticleTypeAccordingToParticleFractions(Double_t trackPt, Double_t jetPt,
1783                                                                             Double_t centralityPercentile,
1784                                                                             Bool_t smearByError,
1785                                                                             Bool_t takeIntoAccountSysError) const
1786 {
1787   // Uses the stored histograms with the particle fractions to generate a random particle type according to these fractions.
1788   // In case of problems (e.g. histo missing), AliPID::kUnknown is returned.
1789   // If smearByError is kTRUE, the used fractions will be random numbers distributed with a gauss with mean
1790   // being the corresponding particle fraction and sigma it's error.
1791   // Note that in this case only the fraction of a random species is varied in this way. The other fractions
1792   // will be re-normalised according their statistical errors.
1793   // The same holds for the systematic error of species "takeIntoAccountSpeciesSysError", but the random number will be
1794   // uniformly distributed within [mean - sys, mean + sys] and the re-normalisation will be weighted with the systematic errors.
1795   // Note that the fractions will be calculated first with only the systematic error taken into account (if desired), including
1796   // re-normalisation. Then, the resulting fractions will be used to calculate the final fractions - either with statistical error
1797   // or without. The species, for which the error will be used for smearing, is the same for sys and stat error.
1798
1799   Double_t prob[AliPID::kSPECIES];
1800   Int_t randomSpecies = (smearByError || takeIntoAccountSysError) ? (Int_t)(fRandom->Rndm() * AliPID::kSPECIES) : -1;
1801   Bool_t success = GetParticleFractions(trackPt, jetPt, centralityPercentile, prob, randomSpecies, randomSpecies);
1802
1803   if (!success)
1804     return AliPID::kUnknown;
1805
1806   Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
1807
1808   if (rnd <= prob[AliPID::kPion])
1809     return AliPID::kPion;
1810   else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon])
1811     return AliPID::kKaon;
1812   else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon] + prob[AliPID::kProton])
1813     return AliPID::kProton;
1814   else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon] + prob[AliPID::kProton] + prob[AliPID::kElectron])
1815     return AliPID::kElectron;
1816
1817   return AliPID::kMuon;  //else it must be a muon (only species left)
1818 }
1819
1820
1821 //_____________________________________________________________________________
1822 AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::GenerateDetectorResponse(AliAnalysisTaskPID::ErrorCode errCode,
1823                                                                            Double_t mean, Double_t sigma,
1824                                                                            Double_t* responses, Int_t nResponses,
1825                                                                            Bool_t usePureGaus)
1826 {
1827   // Generate detector response. If a previous generation was not successful or there is something wrong with this signal generation,
1828   // the function will return kFALSE
1829   if (!responses)
1830     return kError;
1831
1832   // Reset response array
1833   for (Int_t i = 0; i < nResponses; i++)
1834     responses[i] = -999;
1835
1836   if (errCode == kError)
1837     return kError;
1838
1839   ErrorCode ownErrCode = kNoErrors;
1840
1841   if (fUseConvolutedGaus && !usePureGaus) {
1842     // In case of convoluted gauss, calculate the probability density only once to save a lot of time!
1843
1844     TH1* hProbDensity = 0x0;
1845     ownErrCode = SetParamsForConvolutedGaus(mean, sigma);
1846     if (ownErrCode == kError)
1847       return kError;
1848
1849     hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
1850
1851     for (Int_t i = 0; i < nResponses; i++) {
1852       responses[i] = hProbDensity->GetRandom();
1853       //responses[i] fConvolutedGausDeltaPrime->GetRandom(); // MUCH slower than using the binned version via the histogram
1854     }
1855   }
1856   else {
1857     for (Int_t i = 0; i < nResponses; i++) {
1858       responses[i] = fRandom->Gaus(mean, sigma);
1859     }
1860   }
1861
1862   // If forwarded error code was a warning (error case has been handled before), return a warning
1863   if (errCode == kWarning)
1864     return kWarning;
1865
1866   return ownErrCode; // Forward success/warning
1867 }
1868
1869
1870 //_____________________________________________________________________________
1871 void AliAnalysisTaskPID::PrintSettings(Bool_t printSystematicsSettings) const
1872 {
1873   // Print current settings.
1874
1875   printf("\n\nSettings for task %s:\n", GetName());
1876   printf("Is pPb/Pbp: %d -> %s\n", GetIsPbpOrpPb(), GetIsPbpOrpPb() ? "Adapting vertex cuts" : "Using standard vertex cuts");
1877   printf("Track cuts: %s\n", fTrackFilter ? fTrackFilter->GetTitle() : "-");
1878   printf("Eta cut: %.2f <= |eta| <= %.2f\n", GetEtaAbsCutLow(), GetEtaAbsCutUp());
1879   printf("Phi' cut: %d\n", GetUsePhiCut());
1880   printf("TPCCutMIGeo: %d\n", GetUseTPCCutMIGeo());
1881   if (GetUseTPCCutMIGeo()) {
1882     printf("GetCutGeo: %f\n", GetCutGeo());
1883     printf("GetCutNcr: %f\n", GetCutNcr());
1884     printf("GetCutNcl: %f\n", GetCutNcl());
1885   }
1886   printf("TPCnclCut: %d\n", GetUseTPCnclCut());
1887   if (GetUseTPCnclCut()) {
1888     printf("GetCutPureNcl: %d\n", GetCutPureNcl());
1889   }
1890
1891   printf("\n");
1892
1893   printf("Centrality estimator: %s\n", GetCentralityEstimator().Data());
1894
1895   printf("\n");
1896
1897   printf("Use MC-ID for signal generation: %d\n", GetUseMCidForGeneration());
1898   printf("Use ITS: %d\n", GetUseITS());
1899   printf("Use TOF: %d\n", GetUseTOF());
1900   printf("Use priors: %d\n", GetUsePriors());
1901   printf("Use TPC default priors: %d\n", GetUseTPCDefaultPriors());
1902   printf("Use convoluted Gauss: %d\n", GetUseConvolutedGaus());
1903   printf("Accuracy of non-Gaussian tail: %e\n", GetAccuracyNonGaussianTail());
1904   printf("Take into account muons: %d\n", GetTakeIntoAccountMuons());
1905   printf("\nParams for transition from gauss to asymmetric shape:\n");
1906   printf("[0]: %e\n", GetConvolutedGaussTransitionPar(0));
1907   printf("[1]: %e\n", GetConvolutedGaussTransitionPar(1));
1908   printf("[2]: %e\n", GetConvolutedGaussTransitionPar(2));
1909
1910   printf("\n");
1911
1912   printf("Do PID: %d\n", fDoPID);
1913   printf("Do Efficiency: %d\n", fDoEfficiency);
1914   printf("Do PtResolution: %d\n", fDoPtResolution);
1915
1916   printf("\n");
1917
1918   printf("Input from other task: %d\n", GetInputFromOtherTask());
1919   printf("Store additional jet information: %d\n", GetStoreAdditionalJetInformation());
1920   printf("Store centrality percentile: %d", GetStoreCentralityPercentile());
1921
1922   if (printSystematicsSettings)
1923     PrintSystematicsSettings();
1924   else
1925     printf("\n\n\n");
1926 }
1927
1928
1929 //_____________________________________________________________________________
1930 void AliAnalysisTaskPID::PrintSystematicsSettings() const
1931 {
1932   // Print current settings for systematic studies.
1933
1934   printf("\n\nSettings for systematics for task %s:\n", GetName());
1935   printf("Splines:\t%f\n", GetSystematicScalingSplines());
1936   printf("EtaCorrMomThr:\t%f\n", GetSystematicScalingEtaCorrectionMomentumThr());
1937   printf("EtaCorrLowP:\t%f\n", GetSystematicScalingEtaCorrectionLowMomenta());
1938   printf("EtaCorrHighP:\t%f\n", GetSystematicScalingEtaCorrectionHighMomenta());
1939   printf("SigmaPara:\t%f\n", GetSystematicScalingEtaSigmaPara());
1940   printf("MultCorr:\t%f\n", GetSystematicScalingMultCorrection());
1941
1942   printf("\n\n");
1943 }
1944
1945
1946 //_____________________________________________________________________________
1947 Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePDGcode, Double_t centralityPercentile,
1948                                         Double_t jetPt)
1949 {
1950   // Process the track (generate expected response, fill histos, etc.).
1951   // particlePDGcode == 0 means data. Otherwise, the corresponding MC ID will be assumed.
1952
1953   //Printf("Debug: Task %s is starting to process track: dEdx %f, pTPC %f, eta %f, ncl %d\n", GetName(), track->GetTPCsignal(), track->GetTPCmomentum(),
1954   //       track->Eta(), track->GetTPCsignalN());
1955
1956   if(fDebug > 1)
1957     printf("File: %s, Line: %d: ProcessTrack\n", (char*)__FILE__, __LINE__);
1958
1959   if (!fDoPID)
1960     return kFALSE;
1961
1962   if(fDebug > 2)
1963     printf("File: %s, Line: %d: ProcessTrack -> Processing started\n", (char*)__FILE__, __LINE__);
1964
1965   const Bool_t isMC = (particlePDGcode == 0) ? kFALSE : kTRUE;
1966
1967   Int_t binMC = -1;
1968
1969   if (isMC) {
1970     if (TMath::Abs(particlePDGcode) == 211) {//Pion
1971       binMC = 3;
1972     }
1973     else if (TMath::Abs(particlePDGcode) == 321) {//Kaon
1974       binMC = 1;
1975     }
1976     else if (TMath::Abs(particlePDGcode) == 2212) {//Proton
1977       binMC = 4;
1978     }
1979     else if (TMath::Abs(particlePDGcode) == 11) {//Electron
1980       binMC = 0;
1981     }
1982     else if (TMath::Abs(particlePDGcode) == 13) {//Muon
1983       binMC = 2;
1984     }
1985     else // In MC-ID case, set to underflow bin such that the response from this track is only used for unidentified signal generation
1986          // or signal generation with PID response and the track is still there (as in data) - e.g. relevant w.r.t. deuterons.
1987          // This is important to be as much as possible consistent with data. And the tracks can still be removed by disabling the
1988          // underflow bin for the projections
1989       binMC = -1;
1990   }
1991
1992   // Momenta
1993   //Double_t p = track->GetP();
1994   //Double_t pTPC = track->GetTPCmomentum();
1995   Double_t pT = track->Pt();
1996
1997   Double_t z = -1, xi = -1;
1998   GetJetTrackObservables(pT, jetPt, z, xi);
1999
2000
2001   Double_t trackCharge = track->Charge();
2002
2003   // TPC signal
2004   Double_t dEdxTPC = fPIDResponse->IsTunedOnData() ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
2005
2006   if (dEdxTPC <= 0) {
2007     Printf("Skipping track with strange dEdx value: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), track->GetTPCmomentum(),
2008            track->Eta(), track->GetTPCsignalN());
2009     return kFALSE;
2010   }
2011
2012
2013
2014
2015   Double_t dEdxEl, dEdxKa, dEdxPi, dEdxMu, dEdxPr;
2016   Double_t sigmaEl, sigmaKa, sigmaPi, sigmaMu, sigmaPr;
2017
2018   if (fDoAnySystematicStudiesOnTheExpectedSignal) {
2019     // Get the uncorrected signal first and the corresponding correction factors.
2020     // Then modify the correction factors and properly recalculate the corrected dEdx
2021
2022     // Get pure spline values for dEdx_expected, without any correction
2023     dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2024     dEdxKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2025     dEdxPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2026     dEdxMu = !fTakeIntoAccountMuons ? -1 :
2027                             fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2028     dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2029
2030     // Scale splines, if desired
2031     if (TMath::Abs(fSystematicScalingSplines - 1.0) > fgkEpsilon) {
2032       dEdxEl *= fSystematicScalingSplines;
2033       dEdxKa *= fSystematicScalingSplines;
2034       dEdxPi *= fSystematicScalingSplines;
2035       dEdxMu *= fTakeIntoAccountMuons ? fSystematicScalingSplines : 1.;
2036       dEdxPr *= fSystematicScalingSplines;
2037     }
2038
2039     // Get the eta correction factors for the (modified) expected dEdx
2040     Double_t etaCorrEl = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxEl) : 1.;
2041     Double_t etaCorrKa = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxKa) : 1.;
2042     Double_t etaCorrPi = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPi) : 1.;
2043     Double_t etaCorrMu = fTakeIntoAccountMuons && !fPIDResponse->UseTPCEtaCorrection() ?
2044                             fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxMu) : 1.;
2045     Double_t etaCorrPr = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPr) : 1.;
2046
2047     // Scale eta correction factors, if desired (and eta correction maps are to be used, otherwise it is not possible!)
2048     if (fPIDResponse->UseTPCEtaCorrection() &&
2049         (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon ||
2050          TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon)) {
2051       // Since we do not want to scale the splines with this, but only the eta variation, only scale the deviation of the correction factor!
2052       // 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!
2053
2054
2055       // Due to additional azimuthal effects, there is an additional eta dependence for low momenta which is not corrected successfully so far.
2056       // One can assign a different (higher) systematic scale factor for this low-p region and a threshold which separates low- and high-p.
2057       // 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
2058       Double_t usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionHighMomenta;
2059
2060       if (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - fSystematicScalingEtaCorrectionLowMomenta) > fgkEpsilon) {
2061         const Double_t pTPC = track->GetTPCmomentum();
2062         const Double_t fractionHighMomentumScaleFactor = 0.5 * (1. + TMath::Erf((pTPC - fSystematicScalingEtaCorrectionMomentumThr) / 0.1));
2063         usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionLowMomenta * (1 - fractionHighMomentumScaleFactor)
2064                                              + fSystematicScalingEtaCorrectionHighMomenta * fractionHighMomentumScaleFactor;
2065       }
2066
2067       etaCorrEl = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrEl - 1.0);
2068       etaCorrKa = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrKa - 1.0);
2069       etaCorrPi = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPi - 1.0);
2070       etaCorrMu = fTakeIntoAccountMuons ? (1.0 + usedSystematicScalingEtaCorrection * (etaCorrMu - 1.0)) : 1.0;
2071       etaCorrPr = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPr - 1.0);
2072     }
2073
2074     // Get the multiplicity correction factors for the (modified) expected dEdx
2075     const Int_t currEvtMultiplicity = fPIDResponse->GetTPCResponse().GetCurrentEventMultiplicity();
2076
2077     Double_t multiplicityCorrEl = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2078                                     dEdxEl * etaCorrEl, currEvtMultiplicity) : 1.;
2079     Double_t multiplicityCorrKa = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2080                                     dEdxKa * etaCorrKa, currEvtMultiplicity) : 1.;
2081     Double_t multiplicityCorrPi = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2082                                     dEdxPi * etaCorrPi, currEvtMultiplicity) : 1.;
2083     Double_t multiplicityCorrMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCMultiplicityCorrection() ?
2084                                     fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track, dEdxMu * etaCorrMu, currEvtMultiplicity) : 1.;
2085     Double_t multiplicityCorrPr = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2086                                     dEdxPr * etaCorrPr, currEvtMultiplicity) : 1.;
2087
2088     Double_t multiplicityCorrSigmaEl = fPIDResponse->UseTPCMultiplicityCorrection() ?
2089                                         fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxEl * etaCorrEl, currEvtMultiplicity) : 1.;
2090     Double_t multiplicityCorrSigmaKa = fPIDResponse->UseTPCMultiplicityCorrection() ?
2091                                         fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxKa * etaCorrKa, currEvtMultiplicity) : 1.;
2092     Double_t multiplicityCorrSigmaPi = fPIDResponse->UseTPCMultiplicityCorrection() ?
2093                                         fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxPi * etaCorrPi, currEvtMultiplicity) : 1.;
2094     Double_t multiplicityCorrSigmaMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCMultiplicityCorrection() ?
2095                                         fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxMu * etaCorrMu, currEvtMultiplicity) : 1.;
2096     Double_t multiplicityCorrSigmaPr = fPIDResponse->UseTPCMultiplicityCorrection() ?
2097                                         fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxPr * etaCorrPr, currEvtMultiplicity) : 1.;
2098
2099     // Scale multiplicity correction factors, if desired (and multiplicity correction functions are to be used, otherwise it is not possible!)
2100     if (fPIDResponse->UseTPCMultiplicityCorrection() && TMath::Abs(fSystematicScalingMultCorrection - 1.0) > fgkEpsilon) {
2101       // Since we do not want to scale the splines with this, but only the multiplicity variation, only scale the deviation of the correction factor!
2102       // 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!
2103
2104       multiplicityCorrEl = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrEl - 1.0);
2105       multiplicityCorrKa = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrKa - 1.0);
2106       multiplicityCorrPi = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrPi - 1.0);
2107       multiplicityCorrMu = fTakeIntoAccountMuons ? (1.0 + fSystematicScalingMultCorrection * (multiplicityCorrMu - 1.0)) : 1.0;
2108       multiplicityCorrPr = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrPr - 1.0);
2109     }
2110
2111     // eta correction must be enabled in order to use the new sigma parametrisation maps. Since this is the absolute sigma
2112     // for a track calculated with the unscaled paramaters, we have to devide out dEdxExpectedEtaCorrected and then need
2113     // to scale with the multiplicitySigmaCorrFactor * fSystematicScalingEtaSigmaPara. In the end, one has to scale with the
2114     // (modified) dEdx to get the absolute sigma
2115     // This means there is no extra parameter for the multiplicitySigmaCorrFactor, but only for the sigma map itself.
2116     // This is valid, since it appears only as a product. One has to assume a larger systematic shift in case of additional
2117     // multiplicity dependence....
2118     Double_t sigmaRelEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2119                           / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2120                           * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaEl;
2121
2122     Double_t sigmaRelKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2123                           / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2124                           * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaKa;
2125
2126     Double_t sigmaRelPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2127                           / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2128                           * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaPi;
2129
2130     Double_t sigmaRelMu = fTakeIntoAccountMuons ?
2131                             fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2132                             / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2133                             * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaMu
2134                             : 999.;
2135
2136     Double_t sigmaRelPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2137                           / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2138                           * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaPr;
2139
2140     // Now scale the (possibly modified) spline values with the (possibly modified) correction factors
2141     dEdxEl *= etaCorrEl * multiplicityCorrEl;
2142     dEdxKa *= etaCorrKa * multiplicityCorrKa;
2143     dEdxPi *= etaCorrPi * multiplicityCorrPi;
2144     dEdxMu *= etaCorrMu * multiplicityCorrMu;
2145     dEdxPr *= etaCorrPr * multiplicityCorrPr;
2146
2147     // Finally, get the absolute sigma
2148     sigmaEl = sigmaRelEl * dEdxEl;
2149     sigmaKa = sigmaRelKa * dEdxKa;
2150     sigmaPi = sigmaRelPi * dEdxPi;
2151     sigmaMu = sigmaRelMu * dEdxMu;
2152     sigmaPr = sigmaRelPr * dEdxPr;
2153
2154   }
2155   else {
2156     // No systematic studies on expected signal - just take it as it comve from the TPCPIDResponse
2157     dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2158                                                               fPIDResponse->UseTPCEtaCorrection(),
2159                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2160     dEdxKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2161                                                               fPIDResponse->UseTPCEtaCorrection(),
2162                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2163     dEdxPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2164                                                               fPIDResponse->UseTPCEtaCorrection(),
2165                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2166     dEdxMu = !fTakeIntoAccountMuons ? -1 :
2167                             fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2168                                                                              fPIDResponse->UseTPCEtaCorrection(),
2169                                                                              fPIDResponse->UseTPCMultiplicityCorrection());
2170     dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2171                                                               fPIDResponse->UseTPCEtaCorrection(),
2172                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2173
2174     sigmaEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2175                                                               fPIDResponse->UseTPCEtaCorrection(),
2176                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2177     sigmaKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2178                                                               fPIDResponse->UseTPCEtaCorrection(),
2179                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2180     sigmaPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2181                                                               fPIDResponse->UseTPCEtaCorrection(),
2182                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2183     sigmaMu = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2184                                                               fPIDResponse->UseTPCEtaCorrection(),
2185                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2186     sigmaPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2187                                                               fPIDResponse->UseTPCEtaCorrection(),
2188                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2189   }
2190   /*OLD with deltaSpecies
2191   Double_t deltaElectron = dEdxTPC - dEdxEl;
2192   Double_t deltaKaon = dEdxTPC - dEdxKa;
2193   Double_t deltaPion = dEdxTPC - dEdxPi;
2194   Double_t deltaProton = dEdxTPC - dEdxPr;
2195   */
2196
2197   Double_t deltaPrimeElectron = (dEdxEl > 0) ? dEdxTPC / dEdxEl : -1;
2198   if (dEdxEl <= 0)  {
2199     Printf("Error: Expected TPC signal <= 0 for electron hypothesis");
2200     return kFALSE;
2201   }
2202
2203   Double_t deltaPrimeKaon = (dEdxKa > 0) ? dEdxTPC / dEdxKa : -1;
2204   if (dEdxKa <= 0)  {
2205     Printf("Error: Expected TPC signal <= 0 for kaon hypothesis");
2206     return kFALSE;
2207   }
2208
2209   Double_t deltaPrimePion = (dEdxPi > 0) ? dEdxTPC / dEdxPi : -1;
2210   if (dEdxPi <= 0)  {
2211     Printf("Error: Expected TPC signal <= 0 for pion hypothesis");
2212     return kFALSE;
2213   }
2214
2215   Double_t deltaPrimeProton = (dEdxPr > 0) ? dEdxTPC / dEdxPr : -1;
2216   if (dEdxPr <= 0)  {
2217     Printf("Error: Expected TPC signal <= 0 for proton hypothesis");
2218     return kFALSE;
2219   }
2220
2221   /*TODO for TOF
2222   // TOF signal
2223   Double_t times[AliPID::kSPECIES];
2224   track->GetIntegratedTimes(times);
2225   AliTOFPIDResponse tofPIDResponse = fPIDResponse->GetTOFResponse();
2226   Float_t electronDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kElectron);
2227   Float_t pionDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kPion);
2228   Float_t kaonDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kKaon);
2229   Float_t protonDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kProton);
2230   */
2231
2232   if(fDebug > 2)
2233     printf("File: %s, Line: %d: ProcessTrack -> Compute probabilities\n", (char*)__FILE__, __LINE__);
2234
2235   // Use probabilities to weigh the response generation for the different species.
2236   // Also determine the most probable particle type.
2237   Double_t prob[AliPID::kSPECIESC];
2238   for (Int_t i = 0; i < AliPID::kSPECIESC; i++)
2239     prob[i] = 0;
2240
2241   fPIDcombined->ComputeProbabilities(track, fPIDResponse, prob);
2242
2243   // Bug: One can set the number of species for PIDcombined, but PIDcombined will call PIDresponse, which writes without testing
2244   // the probs for kSPECIESC (including light nuclei) into the array.
2245   // In this case, when only kSPECIES are considered, the probabilities have to be rescaled!
2246   for (Int_t i = AliPID::kSPECIES; i < AliPID::kSPECIESC; i++)
2247     prob[i] = 0;
2248
2249   // If muons are not to be taken into account, just set their probability to zero and normalise the remaining probabilities
2250   if (!fTakeIntoAccountMuons)
2251     prob[AliPID::kMuon] = 0;
2252
2253   Double_t probSum = 0.;
2254   for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2255     probSum += prob[i];
2256
2257   if (probSum > 0) {
2258     for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2259       prob[i] /= probSum;
2260   }
2261
2262   if (!isMC) {
2263     // If there is no MC information, take the most probable species for the ID
2264     Float_t max = 0.;
2265     Int_t maxIndex = -1;
2266     for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2267       if (prob[i] > max) {
2268         max = prob[i];
2269         maxIndex = i;
2270       }
2271     }
2272
2273     // Translate from AliPID numbering to numbering of this class
2274     if (max > 0) {
2275       if (maxIndex == AliPID::kElectron)
2276         binMC = 0;
2277       else if (maxIndex == AliPID::kKaon)
2278         binMC = 1;
2279       else if (maxIndex == AliPID::kMuon)
2280         binMC = 2;
2281       else if (maxIndex == AliPID::kPion)
2282         binMC = 3;
2283       else if (maxIndex == AliPID::kProton)
2284         binMC = 4;
2285       else
2286         binMC = -1;
2287     }
2288     else {
2289       // Only take track into account for expectation values, if valid pid response is available.. Otherwise: Set to underflow bin.
2290       binMC = -1;
2291     }
2292   }
2293
2294   /*
2295   //For testing: Swap p<->pT to analyse pure p-dependence => Needs to be removed later
2296   Double_t temp = pT;
2297   pT = pTPC;
2298   pTPC = temp;
2299   */
2300
2301
2302   Double_t entry[fStoreAdditionalJetInformation ? kDataNumAxes : kDataNumAxes - fgkNumJetAxes];
2303   entry[kDataMCID] = binMC;
2304   entry[kDataSelectSpecies] = 0;
2305   entry[kDataPt] = pT;
2306   entry[kDataDeltaPrimeSpecies] = deltaPrimeElectron;
2307   entry[kDataCentrality] = centralityPercentile;
2308
2309   if (fStoreAdditionalJetInformation) {
2310     entry[kDataJetPt] = jetPt;
2311     entry[kDataZ] = z;
2312     entry[kDataXi] = xi;
2313   }
2314
2315   entry[GetIndexOfChargeAxisData()] = trackCharge;
2316
2317   /*OLD with TOF, p_TPC_Inner and p_vertex and deltaSpecies
2318   // MC PID, SelectSpecies, P(TPC_inner), pT, p(Vertex), deltaSpecies, deltaPrimeSpecies, deltaTOFspecies
2319   Double_t entry[8] = { binMC, 0, pTPC, pT, p, deltaElectron, deltaPrimeElectron, electronDeltaTOF };
2320   */
2321
2322   fhPIDdataAll->Fill(entry);
2323
2324   entry[kDataSelectSpecies] = 1;
2325   //OLD with deltaSpecies entry[5] = deltaKaon;
2326   entry[kDataDeltaPrimeSpecies] = deltaPrimeKaon;
2327   //TODO for TOF entry[7] = kaonDeltaTOF;
2328   fhPIDdataAll->Fill(entry);
2329
2330   entry[kDataSelectSpecies] = 2;
2331   //OLD with deltaSpecies entry[5] = deltaPion;
2332   entry[kDataDeltaPrimeSpecies] = deltaPrimePion;
2333   //TODO for TOF entry[7] = pionDeltaTOF;
2334   fhPIDdataAll->Fill(entry);
2335
2336   entry[kDataSelectSpecies] = 3;
2337   //OLD with deltaSpecies entry[5] = deltaProton;
2338   entry[kDataDeltaPrimeSpecies] = deltaPrimeProton;
2339   //TODO for TOF entry[7] = protonDeltaTOF;
2340   fhPIDdataAll->Fill(entry);
2341
2342
2343   // Construct the expected shape for the transition p -> pT
2344
2345   Double_t genEntry[fStoreAdditionalJetInformation ? kGenNumAxes : kGenNumAxes - fgkNumJetAxes];
2346   genEntry[kGenMCID] = binMC;
2347   genEntry[kGenSelectSpecies] = 0;
2348   genEntry[kGenPt] = pT;
2349   genEntry[kGenDeltaPrimeSpecies] = -999;
2350   genEntry[kGenCentrality] = centralityPercentile;
2351
2352   if (fStoreAdditionalJetInformation) {
2353     genEntry[kGenJetPt] = jetPt;
2354     genEntry[kGenZ] = z;
2355     genEntry[kGenXi] = xi;
2356   }
2357
2358   genEntry[GetIndexOfChargeAxisGen()] = trackCharge;
2359
2360   //OLD with deltaSpecies Double_t genEntry[5] = { binMC, 0, pT, -999, -999 }; // MC PID, SelectSpecies, pT, deltaSpecies, deltaPrimeSpecies
2361
2362   // Generate numGenEntries "responses" with fluctuations around the expected values.
2363   // The higher the (transverse) momentum, the more "responses" will be generated in order not to run out of statistics too fast.
2364   Int_t numGenEntries = 10;
2365   if (pT >= 5)
2366     numGenEntries *= 5;
2367   else if (pT >= 2)
2368     numGenEntries *= 2;
2369
2370   // Jets have even worse statistics, therefore, scale numGenEntries further
2371   if (jetPt >= 40)
2372     numGenEntries *= 20;
2373   else if (jetPt >= 20)
2374     numGenEntries *= 10;
2375   else if (jetPt >= 10)
2376     numGenEntries *= 2;
2377
2378
2379   // Do not generate more entries than available in memory!
2380   if (numGenEntries > fgkMaxNumGenEntries)// fgkMaxNumGenEntries = 1000
2381     numGenEntries = fgkMaxNumGenEntries;
2382
2383   ErrorCode errCode = kNoErrors;
2384
2385   if(fDebug > 2)
2386     printf("File: %s, Line: %d: ProcessTrack -> Generate Responses\n", (char*)__FILE__, __LINE__);
2387
2388   // Electrons
2389   errCode = GenerateDetectorResponse(errCode, 1.,              sigmaEl / dEdxEl, fGenRespElDeltaPrimeEl, numGenEntries);
2390   errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxKa, sigmaEl / dEdxKa, fGenRespElDeltaPrimeKa, numGenEntries);
2391   errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxPi, sigmaEl / dEdxPi, fGenRespElDeltaPrimePi, numGenEntries);
2392   errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxPr, sigmaEl / dEdxPr, fGenRespElDeltaPrimePr, numGenEntries);
2393
2394   // Kaons
2395   errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxEl, sigmaKa / dEdxEl, fGenRespKaDeltaPrimeEl, numGenEntries);
2396   errCode = GenerateDetectorResponse(errCode, 1.,              sigmaKa / dEdxKa, fGenRespKaDeltaPrimeKa, numGenEntries);
2397   errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxPi, sigmaKa / dEdxPi, fGenRespKaDeltaPrimePi, numGenEntries);
2398   errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxPr, sigmaKa / dEdxPr, fGenRespKaDeltaPrimePr, numGenEntries);
2399
2400   // Pions
2401   errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxEl, sigmaPi / dEdxEl, fGenRespPiDeltaPrimeEl, numGenEntries);
2402   errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxKa, sigmaPi / dEdxKa, fGenRespPiDeltaPrimeKa, numGenEntries);
2403   errCode = GenerateDetectorResponse(errCode, 1.,              sigmaPi / dEdxPi, fGenRespPiDeltaPrimePi, numGenEntries);
2404   errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxPr, sigmaPi / dEdxPr, fGenRespPiDeltaPrimePr, numGenEntries);
2405
2406   // Muons, if desired
2407   if (fTakeIntoAccountMuons) {
2408     errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxEl, sigmaMu / dEdxEl, fGenRespMuDeltaPrimeEl, numGenEntries);
2409     errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxKa, sigmaMu / dEdxKa, fGenRespMuDeltaPrimeKa, numGenEntries);
2410     errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxPi, sigmaMu / dEdxPi, fGenRespMuDeltaPrimePi, numGenEntries);
2411     errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxPr, sigmaMu / dEdxPr, fGenRespMuDeltaPrimePr, numGenEntries);
2412   }
2413
2414   // Protons
2415   errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxEl, sigmaPr / dEdxEl, fGenRespPrDeltaPrimeEl, numGenEntries);
2416   errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxKa, sigmaPr / dEdxKa, fGenRespPrDeltaPrimeKa, numGenEntries);
2417   errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxPi, sigmaPr / dEdxPi, fGenRespPrDeltaPrimePi, numGenEntries);
2418   errCode = GenerateDetectorResponse(errCode, 1.,              sigmaPr / dEdxPr, fGenRespPrDeltaPrimePr, numGenEntries);
2419
2420
2421   /*OLD with deltaSpecies
2422   // Delta
2423
2424   // Electrons
2425   errCode = GenerateDetectorResponse(errCode, 0.,              sigmaEl, fGenRespElDeltaEl, numGenEntries, usePureGausForDelta);
2426   errCode = GenerateDetectorResponse(errCode, dEdxEl - dEdxKa, sigmaEl, fGenRespElDeltaKa, numGenEntries, usePureGausForDelta);
2427   errCode = GenerateDetectorResponse(errCode, dEdxEl - dEdxPi, sigmaEl, fGenRespElDeltaPi, numGenEntries, usePureGausForDelta);
2428   errCode = GenerateDetectorResponse(errCode, dEdxEl - dEdxPr, sigmaEl, fGenRespElDeltaPr, numGenEntries, usePureGausForDelta);
2429
2430   // Kaons
2431   errCode = GenerateDetectorResponse(errCode, dEdxKa - dEdxEl, sigmaKa, fGenRespKaDeltaEl, numGenEntries, usePureGausForDelta);
2432   errCode = GenerateDetectorResponse(errCode, 0.,              sigmaKa, fGenRespKaDeltaKa, numGenEntries, usePureGausForDelta);
2433   errCode = GenerateDetectorResponse(errCode, dEdxKa - dEdxPi, sigmaKa, fGenRespKaDeltaPi, numGenEntries, usePureGausForDelta);
2434   errCode = GenerateDetectorResponse(errCode, dEdxKa - dEdxPr, sigmaKa, fGenRespKaDeltaPr, numGenEntries, usePureGausForDelta);
2435
2436   // Pions
2437   errCode = GenerateDetectorResponse(errCode, dEdxPi - dEdxEl, sigmaPi, fGenRespPiDeltaEl, numGenEntries, usePureGausForDelta);
2438   errCode = GenerateDetectorResponse(errCode, dEdxPi - dEdxKa, sigmaPi, fGenRespPiDeltaKa, numGenEntries, usePureGausForDelta);
2439   errCode = GenerateDetectorResponse(errCode, 0.,              sigmaPi, fGenRespPiDeltaPi, numGenEntries, usePureGausForDelta);
2440   errCode = GenerateDetectorResponse(errCode, dEdxPi - dEdxPr, sigmaPi, fGenRespPiDeltaPr, numGenEntries, usePureGausForDelta);
2441
2442   // Muons
2443   errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxEl, sigmaMu, fGenRespMuDeltaEl, numGenEntries, usePureGausForDelta);
2444   errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxKa, sigmaMu, fGenRespMuDeltaKa, numGenEntries, usePureGausForDelta);
2445   errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxPi, sigmaMu, fGenRespMuDeltaPi, numGenEntries, usePureGausForDelta);
2446   errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxPr, sigmaMu, fGenRespMuDeltaPr, numGenEntries, usePureGausForDelta);
2447
2448   // Protons
2449   errCode = GenerateDetectorResponse(errCode, dEdxPr - dEdxEl, sigmaPr, fGenRespPrDeltaEl, numGenEntries, usePureGausForDelta);
2450   errCode = GenerateDetectorResponse(errCode, dEdxPr - dEdxKa, sigmaPr, fGenRespPrDeltaKa, numGenEntries, usePureGausForDelta);
2451   errCode = GenerateDetectorResponse(errCode, dEdxPr - dEdxPi, sigmaPr, fGenRespPrDeltaPi, numGenEntries, usePureGausForDelta);
2452   errCode = GenerateDetectorResponse(errCode, 0.,              sigmaPr, fGenRespPrDeltaPr, numGenEntries, usePureGausForDelta);
2453   */
2454
2455   if (errCode != kNoErrors) {
2456     if (errCode == kWarning) {
2457       //Printf("Warning: Questionable detector response for track, most likely due to very low number of PID clusters! Debug output (dEdx_expected, sigma_expected):");
2458     }
2459     else
2460       Printf("Error: Failed to generate detector response for track - skipped! Debug output (dEdx_expected, sigma_expected):");
2461
2462     /*
2463     Printf("Pr: %e, %e", dEdxPr, sigmaPr);
2464     Printf("Pi: %e, %e", dEdxPi, sigmaPi);
2465     Printf("El: %e, %e", dEdxEl, sigmaEl);
2466     Printf("Mu: %e, %e", dEdxMu, sigmaMu);
2467     Printf("Ka: %e, %e", dEdxKa, sigmaKa);
2468     Printf("track: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), track->GetTPCmomentum(), track->Eta(),
2469            track->GetTPCsignalN());
2470     */
2471
2472     if (errCode != kWarning) {
2473       fhSkippedTracksForSignalGeneration->Fill(track->GetTPCmomentum(), track->GetTPCsignalN());
2474       return kFALSE;// Skip generated response in case of error
2475     }
2476   }
2477
2478   for (Int_t n = 0; n < numGenEntries; n++)  {
2479     if (!isMC || !fUseMCidForGeneration) {
2480       // If no MC info is available or shall not be used, use weighting with priors to generate entries for the different species
2481       Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
2482
2483       // Consider generated response as originating from...
2484       if (rnd <= prob[AliPID::kElectron])
2485         genEntry[kGenMCID] = 0; // ... an electron
2486       else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon])
2487         genEntry[kGenMCID] = 1;  // ... a kaon
2488       else if (rnd <=  prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon])
2489         genEntry[kGenMCID] = 2;  // ... a muon -> NOTE: prob[AliPID::kMuon] = 0 in case of fTakeIntoAccountMuons = kFALSE
2490       else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon] + prob[AliPID::kPion])
2491         genEntry[kGenMCID] = 3;  // ... a pion
2492       else
2493         genEntry[kGenMCID] = 4;  // ... a proton
2494     }
2495
2496     // Electrons
2497     genEntry[kGenSelectSpecies] = 0;
2498     //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaEl[n];
2499     genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeEl[n];
2500     fhGenEl->Fill(genEntry);
2501
2502     genEntry[kGenSelectSpecies] = 1;
2503     //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaKa[n];
2504     genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeKa[n];
2505     fhGenEl->Fill(genEntry);
2506
2507     genEntry[kGenSelectSpecies] = 2;
2508     //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaPi[n];
2509     genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePi[n];
2510     fhGenEl->Fill(genEntry);
2511
2512     genEntry[kGenSelectSpecies] = 3;
2513     //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaPr[n];
2514     genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePr[n];
2515     fhGenEl->Fill(genEntry);
2516
2517     // Kaons
2518     genEntry[kGenSelectSpecies] = 0;
2519     //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaEl[n];
2520     genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeEl[n];
2521     fhGenKa->Fill(genEntry);
2522
2523     genEntry[kGenSelectSpecies] = 1;
2524     //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaKa[n];
2525     genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeKa[n];
2526     fhGenKa->Fill(genEntry);
2527
2528     genEntry[kGenSelectSpecies] = 2;
2529     //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaPi[n];
2530     genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePi[n];
2531     fhGenKa->Fill(genEntry);
2532
2533     genEntry[kGenSelectSpecies] = 3;
2534     //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaPr[n];
2535     genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePr[n];
2536     fhGenKa->Fill(genEntry);
2537
2538     // Pions
2539     genEntry[kGenSelectSpecies] = 0;
2540     //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaEl[n];
2541     genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeEl[n];
2542     fhGenPi->Fill(genEntry);
2543
2544     genEntry[kGenSelectSpecies] = 1;
2545     //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaKa[n];
2546     genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeKa[n];
2547     fhGenPi->Fill(genEntry);
2548
2549     genEntry[kGenSelectSpecies] = 2;
2550     //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaPi[n];
2551     genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePi[n];
2552     fhGenPi->Fill(genEntry);
2553
2554     genEntry[kGenSelectSpecies] = 3;
2555     //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaPr[n];
2556     genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePr[n];
2557     fhGenPi->Fill(genEntry);
2558
2559     if (fTakeIntoAccountMuons) {
2560       // Muons, if desired
2561       genEntry[kGenSelectSpecies] = 0;
2562       //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaEl[n];
2563       genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeEl[n];
2564       fhGenMu->Fill(genEntry);
2565
2566       genEntry[kGenSelectSpecies] = 1;
2567       //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaKa[n];
2568       genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeKa[n];
2569       fhGenMu->Fill(genEntry);
2570
2571       genEntry[kGenSelectSpecies] = 2;
2572       //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaPi[n];
2573       genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePi[n];
2574       fhGenMu->Fill(genEntry);
2575
2576       genEntry[kGenSelectSpecies] = 3;
2577       //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaPr[n];
2578       genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePr[n];
2579       fhGenMu->Fill(genEntry);
2580     }
2581
2582     // Protons
2583     genEntry[kGenSelectSpecies] = 0;
2584     //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaEl[n];
2585     genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeEl[n];
2586     fhGenPr->Fill(genEntry);
2587
2588     genEntry[kGenSelectSpecies] = 1;
2589     //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaKa[n];
2590     genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeKa[n];
2591     fhGenPr->Fill(genEntry);
2592
2593     genEntry[kGenSelectSpecies] = 2;
2594     //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaPi[n];
2595     genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePi[n];
2596     fhGenPr->Fill(genEntry);
2597
2598     genEntry[kGenSelectSpecies] = 3;
2599     //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaPr[n];
2600     genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePr[n];
2601     fhGenPr->Fill(genEntry);
2602   }
2603
2604   if(fDebug > 2)
2605     printf("File: %s, Line: %d: ProcessTrack -> Done\n", (char*)__FILE__, __LINE__);
2606
2607   return kTRUE;
2608 }
2609
2610
2611 //_____________________________________________________________________________
2612 Bool_t AliAnalysisTaskPID::SetConvolutedGaussLambdaParameter(Double_t lambda)
2613 {
2614   // Set the lambda parameter of the convolution to the desired value. Automatically
2615   // calculates the parameters for the transition (restricted) gauss -> convoluted gauss.
2616   fConvolutedGaussTransitionPars[2] = lambda;
2617
2618   // Save old parameters and settings of function to restore them later:
2619   Double_t* oldFuncParams = new Double_t[fkConvolutedGausNPar];
2620   for (Int_t i = 0; i < fkConvolutedGausNPar; i++)
2621     oldFuncParams[i] = fConvolutedGausDeltaPrime->GetParameter(i);
2622   Int_t oldFuncNpx = fConvolutedGausDeltaPrime->GetNpx();
2623   Double_t oldFuncRangeLow = 0, oldFuncRangeUp = 100;
2624   fConvolutedGausDeltaPrime->GetRange(oldFuncRangeLow, oldFuncRangeUp);
2625
2626   // Choose some sufficiently large range
2627   const Double_t rangeStart = 0.5;
2628   const Double_t rangeEnd = 2.0;
2629
2630   // To get the parameters for the transition, just choose arbitrarily input parameters for mu and sigma
2631   // (it makes sense to choose typical values). The ratio sigma_gauss / sigma_convolution is independent
2632   // of mu and as well the difference mu_gauss - mu_convolution.
2633   Double_t muInput = 1.0;
2634   Double_t sigmaInput = fgkSigmaReferenceForTransitionPars;
2635
2636
2637   // Step 1: Generate distribution with input parameters
2638   const Int_t nPar = 3;
2639   Double_t inputPar[nPar] = { muInput, sigmaInput, lambda };
2640
2641   TH1D* hInput = new TH1D("hInput", "Input distribution", 2000, rangeStart, rangeEnd);
2642
2643   fConvolutedGausDeltaPrime->SetParameters(inputPar);
2644   fConvolutedGausDeltaPrime->SetRange(rangeStart, rangeEnd);
2645   fConvolutedGausDeltaPrime->SetNpx(2000);
2646
2647   /*OLD
2648   // The resolution and mean of the AliPIDResponse are determined in finite intervals
2649   // of ncl (also second order effects due to finite dEdx and tanTheta).
2650   // Take this into account for the transition parameters via assuming an approximately flat
2651   // distritubtion in ncl in this interval.
2652   // NOTE: The ncl interval should be the same as the one used for the sigma map creation!
2653   const Int_t nclStart = 151;
2654   const Int_t nclEnd = 160;
2655   const Int_t nclSteps = (nclEnd - nclStart) + 1;
2656   for (Int_t ncl = nclStart; ncl <= nclEnd; ncl++) {
2657     // Resolution scales with 1/sqrt(ncl)
2658     fConvolutedGausDeltaPrime->SetParameter(1, inputPar[1] * sqrt(nclEnd) / sqrt(ncl));
2659     TH1* hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
2660
2661     for (Int_t i = 0; i < 50000000 / nclSteps; i++)
2662       hInput->Fill(hProbDensity->GetRandom());
2663   }
2664   */
2665
2666   TH1* hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
2667
2668   for (Int_t i = 0; i < 50000000; i++)
2669     hInput->Fill(hProbDensity->GetRandom());
2670
2671   // Step 2: Fit generated distribution with restricted gaussian
2672   Int_t maxBin = hInput->GetMaximumBin();
2673   Double_t max = hInput->GetBinContent(maxBin);
2674
2675   UChar_t usedBins = 1;
2676   if (maxBin > 1) {
2677     max += hInput->GetBinContent(maxBin - 1);
2678     usedBins++;
2679   }
2680   if (maxBin < hInput->GetNbinsX()) {
2681     max += hInput->GetBinContent(maxBin + 1);
2682     usedBins++;
2683   }
2684   max /= usedBins;
2685
2686   // NOTE: The range (<-> fraction of maximum) should be the same
2687   // as for the spline and eta maps creation
2688   const Double_t lowThreshold = hInput->GetXaxis()->GetBinLowEdge(hInput->FindFirstBinAbove(0.1 * max));
2689   const Double_t highThreshold = hInput->GetXaxis()->GetBinUpEdge(hInput->FindLastBinAbove(0.1 * max));
2690
2691   TFitResultPtr fitResGaussFirstStep = hInput->Fit("gaus", "QNRS", "", lowThreshold, highThreshold);
2692
2693   TFitResultPtr fitResGauss;
2694
2695   if ((Int_t)fitResGaussFirstStep == 0) {
2696     TF1 fGauss("fGauss", "[0]*TMath::Gaus(x, [1], [2], 1)", rangeStart, rangeEnd);
2697     fGauss.SetParameter(0, fitResGaussFirstStep->GetParams()[0]);
2698     fGauss.SetParError(0, fitResGaussFirstStep->GetErrors()[0]);
2699     fGauss.SetParameter(2, fitResGaussFirstStep->GetParams()[2]);
2700     fGauss.SetParError(2, fitResGaussFirstStep->GetErrors()[2]);
2701
2702     fGauss.FixParameter(1, fitResGaussFirstStep->GetParams()[1]);
2703     fitResGauss = hInput->Fit(&fGauss, "QNS", "s", rangeStart, rangeEnd);
2704   }
2705   else {
2706     fitResGauss = hInput->Fit("gaus", "QNRS", "same", rangeStart, rangeEnd);
2707   }
2708   //OLD TFitResultPtr fitResGauss = hInput->Fit("gaus", "QNRS", "", hInput->GetXaxis()->GetBinLowEdge(hInput->FindFirstBinAbove(0.1 * max)),
2709   //                                        hInput->GetXaxis()->GetBinUpEdge(hInput->FindLastBinAbove(0.1 * max)));
2710
2711
2712   // Step 3: Use parameters from gaussian fit to obtain parameters for the transition "restricted gauss" -> "convoluted gauss"
2713
2714   // 3.1 The ratio sigmaInput / sigma_gaussFit ONLY depends on lambda (which is fixed per period) -> Calculate this first
2715   // for an arbitrary (here: typical) sigma. The ratio is then ~the same for ALL sigma for given lambda!
2716   if ((Int_t)fitResGauss != 0) {
2717     AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Gauss Fit failed!\n");
2718     fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2719     fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2720     fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2721
2722     delete hInput;
2723     delete [] oldFuncParams;
2724
2725     return kFALSE;
2726   }
2727
2728   if (fitResGauss->GetParams()[2] <= 0) {
2729     AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Sigma of gauss fit <= 0!\n");
2730     fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2731     fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2732     fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2733
2734     delete hInput;
2735     delete [] oldFuncParams;
2736
2737     return kFALSE;
2738   }
2739
2740   // sigma correction factor
2741   fConvolutedGaussTransitionPars[1] = sigmaInput / fitResGauss->GetParams()[2];
2742
2743   // 3.2 Now that sigma und lambda are determined, one can calculate mu by shifting the maximum to the desired position,
2744   // i.e. the maximum of the input distribution should coincide with that of the re-calculated distribution,
2745   // which is achieved by getting the same mu for the same sigma.
2746   // NOTE: For fixed lambda, the shift is proportional to sigma and independent of mu!
2747   // So, one can calculate the shift for an arbitrary fixed (here: typical)
2748   // sigma and then simply use this shift for any other sigma by scaling it correspondingly!!!
2749
2750   // Mu shift correction:
2751   // Shift in mu (difference between mean of gauss and mean of convolution) is proportional to sigma!
2752   // Thus, choose a reference sigma (typical -> 0.05), such that for arbitrary sigma one can simple scale
2753   // this shift correction with sigma / referenceSigma.
2754   fConvolutedGaussTransitionPars[0] = (fitResGauss->GetParams()[1] - muInput);
2755
2756
2757   /*Changed on 03.07.2013
2758
2759   // Maximum of fConvolutedGausDeltaPrime should agree with maximum of input
2760   Double_t par[nPar] = { fitResGauss->GetParams()[1], // just as a guess of the maximum position
2761                          sigmaInput,
2762                          lambda };
2763
2764   fConvolutedGausDeltaPrime->SetParameters(par);
2765
2766   Double_t maxXInput = fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001, muInput - 3. * sigmaInput),
2767                                                               muInput + 10. * sigmaInput,
2768                                                               0.0001);
2769
2770   // Maximum shifts proportional to sigma and is linear in mu (~mean of gauss)
2771   // Maximum should be typically located within [gaussMean, gaussMean + 3 gaussSigma].
2772   // -> Larger search range for safety reasons (also: sigma and/or mean might be not 100% accurate).
2773   Double_t maxXconvoluted = fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001,
2774                                                                               fitResGauss->GetParams()[1] - 3. * fitResGauss->GetParams()[2]),
2775                                                                    fitResGauss->GetParams()[1] + 10. * fitResGauss->GetParams()[2],
2776                                                                    0.0001);
2777   if (maxXconvoluted <= 0) {
2778     AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Maximum of fConvolutedGausDeltaPrime <= 0!\n");
2779
2780     fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2781     fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2782     fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2783
2784     delete hInput;
2785     delete [] oldFuncParams;
2786
2787     return kFALSE;
2788   }
2789
2790   // maxX perfectly shifts as par[0] (scaled by sigma) -> Can shift maxX to input value.
2791   // Mu shift correction:
2792   fConvolutedGaussTransitionPars[0] = maxXconvoluted - maxXInput;
2793   */
2794
2795
2796
2797   fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2798   fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2799   fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2800
2801   delete hInput;
2802   delete [] oldFuncParams;
2803
2804   return kTRUE;
2805 }
2806
2807
2808 //_____________________________________________________________________________
2809 AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::SetParamsForConvolutedGaus(Double_t gausMean, Double_t gausSigma)
2810 {
2811   // Set parameters for convoluted gauss using parameters for a pure gaussian.
2812   // If SetConvolutedGaussLambdaParameter has not been called before to initialise the translation parameters,
2813   // some default parameters will be used and an error will show up.
2814
2815   if(fDebug > 1)
2816     printf("File: %s, Line: %d: SetParamsForConvolutedGaus: mean %e, sigma %e\n", (char*)__FILE__, __LINE__, gausMean, gausSigma);
2817
2818   if (fConvolutedGaussTransitionPars[1] < -998) {
2819     AliError("Transition parameters not initialised! Default parameters will be used. Please call SetConvolutedGaussLambdaParameter(...) before any calculations!");
2820     SetConvolutedGaussLambdaParameter(2.0);
2821     AliError(Form("Parameters set to:\n[0]: %f\n[1]: %f\n[2]: %f\n", fConvolutedGaussTransitionPars[0],
2822                   fConvolutedGaussTransitionPars[1], fConvolutedGaussTransitionPars[2]));
2823   }
2824
2825   Double_t par[fkConvolutedGausNPar];
2826   par[2] = fConvolutedGaussTransitionPars[2];
2827   par[1] = fConvolutedGaussTransitionPars[1] * gausSigma;
2828   // maxX perfectly shifts as par[0] (scaled by sigma) -> Can shift maxX so that it sits at the right place.
2829   par[0] = gausMean - fConvolutedGaussTransitionPars[0] * par[1] / fgkSigmaReferenceForTransitionPars;
2830
2831   ErrorCode errCode = kNoErrors;
2832   fConvolutedGausDeltaPrime->SetParameters(par);
2833
2834   if(fDebug > 3)
2835     printf("File: %s, Line: %d: SetParamsForConvolutedGaus -> Parameters set to: %e, %e, %e (transition pars: %e, %e, %e, %e)\n",
2836            (char*)__FILE__, __LINE__, par[0], par[1], par[2], fConvolutedGaussTransitionPars[0], fConvolutedGaussTransitionPars[1],
2837            fConvolutedGaussTransitionPars[2], fgkSigmaReferenceForTransitionPars);
2838
2839   fConvolutedGausDeltaPrime->SetNpx(20); // Small value speeds up following algorithm (valid, since extrema far apart)
2840
2841   // Accuracy of 10^-5 is enough to get 0.1% precise peak for MIPS w.r.t. to dEdx = 2000 of protons
2842   // (should boost up the algorithm, because 10^-10 is the default value!)
2843   Double_t maxX= fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001, gausMean - 2. * gausSigma),
2844                                                         gausMean + 6. * gausSigma, 1.0E-5);
2845
2846   const Double_t maximum = fConvolutedGausDeltaPrime->Eval(maxX);
2847   const Double_t maximumFraction = maximum * fAccuracyNonGaussianTail;
2848
2849   // Estimate lower boundary for subsequent search:
2850   Double_t lowBoundSearchBoundLow = TMath::Max(1e-4, maxX - 5. * gausSigma);
2851   Double_t lowBoundSearchBoundUp = maxX;
2852
2853   Bool_t lowerBoundaryFixedAtZero = kFALSE;
2854
2855   while (fConvolutedGausDeltaPrime->Eval(lowBoundSearchBoundLow) >= maximumFraction) {
2856     if (lowBoundSearchBoundLow <= 0) {
2857       // This should only happen to low dEdx particles with very few clusters and therefore large sigma, such that the gauss goes below zero deltaPrime
2858       if (maximum <= 0) { // Something is weired
2859         printf("Error generating signal: maximum is <= 0!\n");
2860         return kError;
2861       }
2862       else {
2863         const Double_t valueAtZero = fConvolutedGausDeltaPrime->Eval(0);
2864         if (valueAtZero / maximum > 0.05) {
2865           // Too large fraction below zero deltaPrime. Signal generation cannot be reliable in this case
2866           printf("Error generating signal: Too large fraction below zero deltaPrime: convGauss(0) / convGauss(max) =  %e / %e = %e!\n",
2867                  valueAtZero, maximum, valueAtZero / maximum);
2868           return kError;
2869         }
2870       }
2871
2872       /*
2873       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",
2874              fConvolutedGausDeltaPrime->Eval(maxX), fAccuracyNonGaussianTail, maximumFraction, maxX, par[0], par[1], par[2], gausMean, gausSigma);
2875       */
2876
2877       lowerBoundaryFixedAtZero = kTRUE;
2878
2879       if (errCode != kError)
2880         errCode = kWarning;
2881
2882       break;
2883     }
2884
2885     lowBoundSearchBoundUp -= gausSigma;
2886     lowBoundSearchBoundLow -= gausSigma;
2887
2888     if (lowBoundSearchBoundLow < 0) {
2889       lowBoundSearchBoundLow = 0;
2890       lowBoundSearchBoundUp += gausSigma;
2891     }
2892   }
2893
2894   // Determine lower boundary inside estimated range. For small values of the maximum: Need more precision, since finer binning!
2895   Double_t rangeStart = lowerBoundaryFixedAtZero ? 0 :
2896                         fConvolutedGausDeltaPrime->GetX(maximumFraction, lowBoundSearchBoundLow, lowBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
2897
2898   // .. and the same for the upper boundary
2899   Double_t rangeEnd = 0;
2900   // If distribution starts beyond upper boundary, everything ends up in the overflow bin. So, just reduce range and Npx to minimum
2901   if (rangeStart > fkDeltaPrimeUpLimit) {
2902     rangeEnd = rangeStart + 0.00001;
2903     fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
2904     fConvolutedGausDeltaPrime->SetNpx(4);
2905   }
2906   else {
2907     // Estimate upper boundary for subsequent search:
2908     Double_t upBoundSearchBoundUp = maxX + 5 * gausSigma;
2909     Double_t upBoundSearchBoundLow = maxX;
2910     while (fConvolutedGausDeltaPrime->Eval(upBoundSearchBoundUp) >= maximumFraction) {
2911       upBoundSearchBoundUp += gausSigma;
2912       upBoundSearchBoundLow += gausSigma;
2913     }
2914
2915     //  For small values of the maximum: Need more precision, since finer binning!
2916     rangeEnd = fConvolutedGausDeltaPrime->GetX(maximumFraction, upBoundSearchBoundLow, upBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
2917
2918     fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
2919     fConvolutedGausDeltaPrime->SetNpx(fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindBin(rangeEnd)
2920                                       - fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindBin(rangeStart) + 1);
2921     //fConvolutedGausDeltaPrime->SetNpx((rangeEnd - rangeStart) / fDeltaPrimeBinWidth + 1);
2922   }
2923
2924   if(fDebug > 3)
2925     printf("File: %s, Line: %d: SetParamsForConvolutedGaus -> range %f - %f, error code %d\n", (char*)__FILE__, __LINE__,
2926            rangeStart, rangeEnd, errCode);
2927
2928   return errCode;
2929 }
2930
2931
2932 //________________________________________________________________________
2933 void AliAnalysisTaskPID::SetUpGenHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const
2934 {
2935   // Sets bin limits for axes which are not standard binned and the axes titles.
2936
2937   hist->SetBinEdges(kGenPt, binsPt);
2938   hist->SetBinEdges(kGenDeltaPrimeSpecies, binsDeltaPrime);
2939   hist->SetBinEdges(kGenCentrality, binsCent);
2940
2941   if (fStoreAdditionalJetInformation)
2942     hist->SetBinEdges(kGenJetPt, binsJetPt);
2943
2944   // Set axes titles
2945   hist->GetAxis(kGenMCID)->SetTitle("MC PID");
2946   hist->GetAxis(kGenMCID)->SetBinLabel(1, "e");
2947   hist->GetAxis(kGenMCID)->SetBinLabel(2, "K");
2948   hist->GetAxis(kGenMCID)->SetBinLabel(3, "#mu");
2949   hist->GetAxis(kGenMCID)->SetBinLabel(4, "#pi");
2950   hist->GetAxis(kGenMCID)->SetBinLabel(5, "p");
2951
2952   hist->GetAxis(kGenSelectSpecies)->SetTitle("Select Species");
2953   hist->GetAxis(kGenSelectSpecies)->SetBinLabel(1, "e");
2954   hist->GetAxis(kGenSelectSpecies)->SetBinLabel(2, "K");
2955   hist->GetAxis(kGenSelectSpecies)->SetBinLabel(3, "#pi");
2956   hist->GetAxis(kGenSelectSpecies)->SetBinLabel(4, "p");
2957
2958   hist->GetAxis(kGenPt)->SetTitle("P_{T} (GeV/c)");
2959
2960   hist->GetAxis(kGenDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
2961
2962   hist->GetAxis(kGenCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
2963
2964   if (fStoreAdditionalJetInformation) {
2965     hist->GetAxis(kGenJetPt)->SetTitle("P_{T}^{jet} (GeV/c)");
2966
2967     hist->GetAxis(kGenZ)->SetTitle("z = P_{T}^{track} / P_{T}^{jet}");
2968
2969     hist->GetAxis(kGenXi)->SetTitle("#xi = ln(P_{T}^{jet} / P_{T}^{track})");
2970   }
2971
2972   hist->GetAxis(GetIndexOfChargeAxisGen())->SetTitle("Charge (e_{0})");
2973 }
2974
2975
2976 //________________________________________________________________________
2977 void AliAnalysisTaskPID::SetUpGenYieldHist(THnSparse* hist, Double_t* binsPt, Double_t* binsCent, Double_t* binsJetPt) const
2978 {
2979   // Sets bin limits for axes which are not standard binned and the axes titles.
2980
2981   hist->SetBinEdges(kGenYieldPt, binsPt);
2982   hist->SetBinEdges(kGenYieldCentrality, binsCent);
2983   if (fStoreAdditionalJetInformation)
2984     hist->SetBinEdges(kGenYieldJetPt, binsJetPt);
2985
2986   for (Int_t i = 0; i < 5; i++)
2987     hist->GetAxis(kGenYieldMCID)->SetBinLabel(i + 1, AliPID::ParticleLatexName(i));
2988
2989   // Set axes titles
2990   hist->GetAxis(kGenYieldMCID)->SetTitle("MC PID");
2991   hist->GetAxis(kGenYieldPt)->SetTitle("P_{T}^{gen} (GeV/c)");
2992   hist->GetAxis(kGenYieldCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
2993
2994   if (fStoreAdditionalJetInformation) {
2995     hist->GetAxis(kGenYieldJetPt)->SetTitle("P_{T}^{jet, gen} (GeV/c)");
2996
2997     hist->GetAxis(kGenYieldZ)->SetTitle("z = P_{T}^{track} / P_{T}^{jet}");
2998
2999     hist->GetAxis(kGenYieldXi)->SetTitle("#xi = ln(P_{T}^{jet} / P_{T}^{track})");
3000   }
3001
3002   hist->GetAxis(GetIndexOfChargeAxisGenYield())->SetTitle("Charge (e_{0})");
3003 }
3004
3005
3006 //________________________________________________________________________
3007 void AliAnalysisTaskPID::SetUpHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const
3008 {
3009   // Sets bin limits for axes which are not standard binned and the axes titles.
3010
3011   hist->SetBinEdges(kDataPt, binsPt);
3012   hist->SetBinEdges(kDataDeltaPrimeSpecies, binsDeltaPrime);
3013   hist->SetBinEdges(kDataCentrality, binsCent);
3014
3015   if (fStoreAdditionalJetInformation)
3016     hist->SetBinEdges(kDataJetPt, binsJetPt);
3017
3018   // Set axes titles
3019   hist->GetAxis(kDataMCID)->SetTitle("MC PID");
3020   hist->GetAxis(kDataMCID)->SetBinLabel(1, "e");
3021   hist->GetAxis(kDataMCID)->SetBinLabel(2, "K");
3022   hist->GetAxis(kDataMCID)->SetBinLabel(3, "#mu");
3023   hist->GetAxis(kDataMCID)->SetBinLabel(4, "#pi");
3024   hist->GetAxis(kDataMCID)->SetBinLabel(5, "p");
3025
3026   hist->GetAxis(kDataSelectSpecies)->SetTitle("Select Species");
3027   hist->GetAxis(kDataSelectSpecies)->SetBinLabel(1, "e");
3028   hist->GetAxis(kDataSelectSpecies)->SetBinLabel(2, "K");
3029   hist->GetAxis(kDataSelectSpecies)->SetBinLabel(3, "#pi");
3030   hist->GetAxis(kDataSelectSpecies)->SetBinLabel(4, "p");
3031
3032   hist->GetAxis(kDataPt)->SetTitle("P_{T} (GeV/c)");
3033
3034   hist->GetAxis(kDataDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
3035
3036   hist->GetAxis(kDataCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
3037
3038   if (fStoreAdditionalJetInformation) {
3039     hist->GetAxis(kDataJetPt)->SetTitle("P_{T}^{jet} (GeV/c)");
3040
3041     hist->GetAxis(kDataZ)->SetTitle("z = P_{T}^{track} / P_{T}^{jet}");
3042
3043     hist->GetAxis(kDataXi)->SetTitle("#xi = ln(P_{T}^{jet} / P_{T}^{track})");
3044   }
3045
3046   hist->GetAxis(GetIndexOfChargeAxisData())->SetTitle("Charge (e_{0})");
3047
3048   /*OLD with TOF, p_TPC_Inner and p_vertex
3049   // MC PID, SelectSpecies, P(TPC_inner), pT, p(Vertex), deltaSpecies, deltaPrimeSpecies, deltaTOFspecies
3050   hist->SetBinEdges(2, binsPt);
3051   hist->SetBinEdges(3, binsPt);
3052   hist->SetBinEdges(4, binsPt);
3053
3054   // Set axes titles
3055   hist->GetAxis(0)->SetTitle("MC PID");
3056   hist->GetAxis(0)->SetBinLabel(1, "e");
3057   hist->GetAxis(0)->SetBinLabel(2, "K");
3058   hist->GetAxis(0)->SetBinLabel(3, "#mu");
3059   hist->GetAxis(0)->SetBinLabel(4, "#pi");
3060   hist->GetAxis(0)->SetBinLabel(5, "p");
3061
3062   hist->GetAxis(1)->SetTitle("Select Species");
3063   hist->GetAxis(1)->SetBinLabel(1, "e");
3064   hist->GetAxis(1)->SetBinLabel(2, "K");
3065   hist->GetAxis(1)->SetBinLabel(3, "#pi");
3066   hist->GetAxis(1)->SetBinLabel(4, "p");
3067
3068   hist->GetAxis(2)->SetTitle("P_{TPC_inner} (GeV/c)");
3069   hist->GetAxis(3)->SetTitle("P_{T} (GeV/c)");
3070   hist->GetAxis(4)->SetTitle("P_{vertex} (GeV/c)");
3071
3072   hist->GetAxis(5)->SetTitle("TPC #Delta_{species} (arb. unit)");
3073
3074   hist->GetAxis(6)->SetTitle("TPC #Delta'_{species} (arb. unit)");
3075
3076   hist->GetAxis(7)->SetTitle("#Delta TOF_{species} (ps)");
3077   */
3078 }
3079
3080
3081 //________________________________________________________________________
3082 void AliAnalysisTaskPID::SetUpPtResHist(THnSparse* hist, Double_t* binsPt, Double_t* binsJetPt, Double_t* binsCent) const
3083 {
3084   // Sets bin limits for axes which are not standard binned and the axes titles.
3085
3086   hist->SetBinEdges(kPtResJetPt, binsJetPt);
3087   hist->SetBinEdges(kPtResGenPt, binsPt);
3088   hist->SetBinEdges(kPtResRecPt, binsPt);
3089   hist->SetBinEdges(kPtResCentrality, binsCent);
3090
3091   // Set axes titles
3092   hist->GetAxis(kPtResJetPt)->SetTitle("P_{T}^{jet, rec} (GeV/c)");
3093   hist->GetAxis(kPtResGenPt)->SetTitle("P_{T}^{gen} (GeV/c)");
3094   hist->GetAxis(kPtResRecPt)->SetTitle("P_{T}^{rec} (GeV/c)");
3095
3096   hist->GetAxis(kPtResCharge)->SetTitle("Charge (e_{0})");
3097   hist->GetAxis(kPtResCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
3098 }