PWGJE/UserTasks/AliAnalysisTaskPID.cxx

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