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