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