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