PWGJE/UserTasks/AliAnalysisTaskPID.cxx

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