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 (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.645162;
1881     else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.627431;
1882     else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.457136;
1883     else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.384369;
1884     else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.330597;
1885     else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.309571;
1886     else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.293620;
1887     else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.283709;
1888     else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.282047;
1889     else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.277261;
1890     else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.275772;
1891     else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.280726;
1892     else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.288540;
1893     else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.288315;
1894     else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.296619;
1895     else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.302993;
1896     else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.338121;
1897     else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.349800;
1898     else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.356802;
1899     else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.391202;
1900     else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.422573;
1901     else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.573815;
1902     else if(6.00 <= motherGenPt && motherGenPt < 8.00) fac = 0.786984;
1903     else if(8.00 <= motherGenPt && motherGenPt < 10.00) fac = 1.020021;
1904   }
1905
1906   if (absMotherPDG == 3312 || absMotherPDG == 3322) { // xi
1907     if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.666620;
1908     else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.575908;
1909     else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.433198;
1910     else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.340901;
1911     else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.290896;
1912     else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.236074;
1913     else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.218681;
1914     else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.207763;
1915     else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.222848;
1916     else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.208806;
1917     else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.197275;
1918     else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.183645;
1919     else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.188788;
1920     else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.188282;
1921     else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.207442;
1922     else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.240388;
1923     else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.241916;
1924     else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.208276;
1925     else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.234550;
1926     else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.251689;
1927     else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.310204;
1928     else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.343492;
1929   }
1930
1931   const Double_t weight = 1. / fac;
1932
1933   return weight;
1934 }
1935
1936
1937 //_____________________________________________________________________________
1938 Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(AliMCEvent* mcEvent, AliMCParticle* daughter)
1939 {
1940   // Strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
1941   // -> Based on function in PWGJE/AliAnalysisTaskFragmentationFunction
1942
1943   if (!mcEvent)
1944     return 1.;
1945
1946   AliMCParticle* currentMother   = daughter;
1947   AliMCParticle* currentDaughter = daughter;
1948
1949
1950   // find first primary mother K0s, Lambda or Xi
1951   while(1) {
1952     Int_t daughterPDG = currentDaughter->PdgCode();
1953
1954     Int_t motherLabel = currentDaughter->GetMother();
1955     if(motherLabel >= mcEvent->GetNumberOfTracks()){ // protection
1956       currentMother = currentDaughter;
1957       break;
1958     }
1959
1960     currentMother = (AliMCParticle*)mcEvent->GetTrack(motherLabel);
1961
1962     if (!currentMother) {
1963       currentMother = currentDaughter;
1964       break;
1965     }
1966
1967     Int_t motherPDG = currentMother->PdgCode();
1968
1969     // phys. primary found ?
1970     if (mcEvent->IsPhysicalPrimary(motherLabel))
1971       break;
1972
1973     if (TMath::Abs(daughterPDG) == 321) {
1974       // K+/K- e.g. from phi (ref data not feeddown corrected)
1975       currentMother = currentDaughter;
1976       break;
1977     }
1978     if (TMath::Abs(motherPDG) == 310) {
1979       // K0s e.g. from phi (ref data not feeddown corrected)
1980       break;
1981     }
1982     if (TMath::Abs(motherPDG) == 3212 && TMath::Abs(daughterPDG) == 3122) {
1983       // Mother Sigma0, daughter Lambda (this case not included in feeddown corr.)
1984       currentMother = currentDaughter;
1985       break;
1986     }
1987
1988     currentDaughter = currentMother;
1989   }
1990
1991
1992   Int_t motherPDG = currentMother->PdgCode();
1993   Double_t motherGenPt = currentMother->Pt();
1994
1995   return GetMCStrangenessFactorCMS(motherPDG, motherGenPt);
1996 }
1997
1998
1999 // _________________________________________________________________________________
2000 AliAnalysisTaskPID::TOFpidInfo AliAnalysisTaskPID::GetTOFType(const AliVTrack* track, Int_t tofMode) const
2001 {
2002   // Get the (locally defined) particle type judged by TOF
2003
2004   if (!fPIDResponse) {
2005     Printf("ERROR: fPIDResponse not available -> Cannot determine TOF type!");
2006     return kNoTOFinfo;
2007   }
2008
2009   /*TODO still needs some further thinking how to implement it....
2010   // TOF PID status (kTOFout, kTIME) automatically checked by return value of ComputeTPCProbability;
2011   // also, probability array will be set there (no need to initialise here)
2012   Double_t p[AliPID::kSPECIES];
2013   const AliPIDResponse::EDetPidStatus tofStatus = fPIDResponse->ComputeTPCProbability(track, AliPID::kSPECIES, p);
2014   if (tofStatus != AliPIDResponse::kDetPidOk)
2015     return kNoTOFinfo;
2016
2017   // Do not consider muons
2018   p[AliPID::kMuon] = 0.;
2019
2020   // Probabilities are not normalised yet
2021   Double_t sum = 0.;
2022   for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2023     sum += p[i];
2024
2025   if (sum <= 0.)
2026     return kNoTOFinfo;
2027
2028   for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2029     p[i] /= sum;
2030
2031   Double_t probThreshold = -999.;
2032
2033   // If there is only one distribution, the threshold corresponds to...
2034   if (tofMode == 0) {
2035     probThreshold = ;
2036   }
2037   else if (tofMode == 1) { // default
2038     probThreshold = 0.9973; // a 3-sigma inclusion cut
2039   }
2040   else if (tofMode == 2) {
2041     inclusion = 3.;
2042     exclusion = 3.5;
2043   }
2044   else {
2045     Printf("ERROR: Bad TOF mode: %d!", tofMode);
2046     return kNoTOFinfo;
2047   }
2048
2049   */
2050
2051   ///* OLD: cut with n-sigmas (ATTENTION: Does not take into account properly the TOF tail!)
2052   // Check kTOFout, kTIME, mismatch
2053   const AliPIDResponse::EDetPidStatus tofStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, track);
2054   if (tofStatus != AliPIDResponse::kDetPidOk)
2055     return kNoTOFinfo;
2056
2057   Double_t nsigma[kNumTOFspecies + 1] = { -999., -999., -999., -999. };
2058   nsigma[kTOFpion]   = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion);
2059   nsigma[kTOFkaon]   = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon);
2060   nsigma[kTOFproton] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton);
2061
2062   Double_t inclusion = -999;
2063   Double_t exclusion = -999;
2064
2065   if (tofMode == 0) {
2066     inclusion = 3.;
2067     exclusion = 2.5;
2068   }
2069   else if (tofMode == 1) { // default
2070     inclusion = 3.;
2071     exclusion = 3.;
2072   }
2073   else if (tofMode == 2) {
2074     inclusion = 3.;
2075     exclusion = 3.5;
2076   }
2077   else {
2078     Printf("ERROR: Bad TOF mode: %d!", tofMode);
2079     return kNoTOFinfo;
2080   }
2081
2082   // Do not cut on nSigma electron because this would also remove pions in a large pT range.
2083   // The electron contamination of the pion sample can be treated by TPC dEdx fits afterwards.
2084   if (TMath::Abs(nsigma[kTOFpion]) < inclusion && TMath::Abs(nsigma[kTOFkaon]) > exclusion && TMath::Abs(nsigma[kTOFproton]) > exclusion)
2085     return kTOFpion;
2086   if (TMath::Abs(nsigma[kTOFpion]) > exclusion && TMath::Abs(nsigma[kTOFkaon]) < inclusion && TMath::Abs(nsigma[kTOFproton]) > exclusion)
2087     return kTOFkaon;
2088   if (TMath::Abs(nsigma[kTOFpion]) > exclusion && TMath::Abs(nsigma[kTOFkaon]) > exclusion && TMath::Abs(nsigma[kTOFproton]) < inclusion)
2089     return kTOFproton;
2090   //*/
2091
2092   // There are no TOF electrons selected because the purity is rather bad, even for small momenta
2093   // (also a small mismatch probability significantly affects electrons because their fraction
2094   // is small). There is no need for TOF electrons anyway, since the dEdx distribution of kaons and
2095   // protons in a given pT bin (also at the dEdx crossings) is very different from that of electrons.
2096   // This is due to the steeply falling dEdx of p and K with momentum, whereas the electron dEdx stays
2097   // rather constant.
2098   // As a result, the TPC fit yields a more accurate electron fraction than the TOF selection can do.
2099
2100   return kNoTOFpid;
2101 }
2102
2103
2104 // _________________________________________________________________________________
2105 Bool_t AliAnalysisTaskPID::IsSecondaryWithStrangeMotherMC(AliMCEvent* mcEvent, Int_t partLabel)
2106 {
2107   // Check whether particle is a secondary with strange mother, i.e. returns kTRUE if a strange mother is found
2108   // and the particle is NOT a physical primary. In all other cases kFALSE is returned
2109
2110   if (!mcEvent || partLabel < 0)
2111     return kFALSE;
2112
2113   AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(partLabel);
2114
2115   if (!part)
2116     return kFALSE;
2117
2118   if (mcEvent->IsPhysicalPrimary(partLabel))
2119     return kFALSE;
2120
2121   Int_t iMother = part->GetMother();
2122   if (iMother < 0)
2123     return kFALSE;
2124
2125
2126   AliMCParticle* partM = (AliMCParticle*)mcEvent->GetTrack(iMother);
2127   if (!partM)
2128    return kFALSE;
2129
2130   Int_t codeM = TMath::Abs(partM->PdgCode());
2131   Int_t mfl = Int_t(codeM / TMath::Power(10, Int_t(TMath::Log10(codeM))));
2132   if (mfl == 3 && codeM != 3) // codeM = 3 is for s quark
2133     return kTRUE;
2134
2135   return kFALSE;
2136 }
2137
2138
2139 //_____________________________________________________________________________
2140 Bool_t AliAnalysisTaskPID::SetParticleFractionHisto(const TH3D* hist, Int_t species, Bool_t sysError)
2141 {
2142   // Store a clone of hist (containing the particle fractions of the corresponding species with statistical error (sysError = kFALSE)
2143   // or systematic error (sysError = kTRUE), respectively), internally
2144
2145   if (species < AliPID::kElectron || species > AliPID::kProton) {
2146     AliError(Form("Only fractions for species index %d to %d can be set, but not for the requested one: %d", 0,
2147                   AliPID::kProton, species));
2148     return kFALSE;
2149   }
2150
2151   if (sysError) {
2152     delete fFractionSysErrorHists[species];
2153
2154     fFractionSysErrorHists[species] = new TH3D(*hist);
2155   }
2156   else {
2157     delete fFractionHists[species];
2158
2159     fFractionHists[species] = new TH3D(*hist);
2160   }
2161
2162   return kTRUE;
2163 }
2164
2165
2166 //_____________________________________________________________________________
2167 Bool_t AliAnalysisTaskPID::SetParticleFractionHistosFromFile(const TString filePathName, Bool_t sysError)
2168 {
2169   // Loads particle fractions for all species from the desired file and returns kTRUE on success.
2170   // The maps are assumed to be of Type TH3D, to sit in the main directory and to have names
2171   // Form("hFraction_%e", AliPID::ParticleName(i)) for sysError = kFALSE and
2172   // Form("hFractionSysError_%e", AliPID::ParticleName(i)) for sysError = kTRUE.
2173
2174   TFile* f = TFile::Open(filePathName.Data());
2175   if (!f)  {
2176     std::cout << "Failed to open file with particle fractions \"" << filePathName.Data() << "\"!" << std::endl;
2177     return kFALSE;
2178   }
2179
2180   TH3D* hist = 0x0;
2181   for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2182     TString histName = Form("hFraction%s_%s", sysError ? "SysError" : "", AliPID::ParticleName(i));
2183     hist = dynamic_cast<TH3D*>(f->Get(histName.Data()));
2184     if (!hist) {
2185       std::cout << "Failed to load particle fractions for " << histName.Data() << "!";
2186       std::cout << std::endl << "Cleaning up particle fraction histos!" << std::endl;
2187       CleanupParticleFractionHistos();
2188       return kFALSE;
2189     }
2190
2191     if (!SetParticleFractionHisto(hist, i, sysError)) {
2192       std::cout << "Failed to load particle fractions for " << histName.Data() << "!";
2193       std::cout << std::endl << "Cleaning up particle fraction histos!" << std::endl;
2194       CleanupParticleFractionHistos();
2195       return kFALSE;
2196     }
2197   }
2198
2199   delete hist;
2200
2201   return kTRUE;
2202
2203 }
2204
2205
2206 //_____________________________________________________________________________
2207 Double_t AliAnalysisTaskPID::GetMaxEtaVariation(Double_t dEdxSplines)
2208 {
2209   // Returns the maximum eta variation (i.e. deviation of eta correction factor from unity) for the
2210   // given (spline) dEdx
2211
2212   if (dEdxSplines < 1. || !fhMaxEtaVariation) {
2213     Printf("Error GetMaxEtaVariation: No map or invalid dEdxSplines (%f)!", dEdxSplines);
2214     return 999.;
2215   }
2216
2217   Int_t bin = fhMaxEtaVariation->GetXaxis()->FindFixBin(1. / dEdxSplines);
2218
2219   if (bin == 0)
2220     bin = 1;
2221   if (bin > fhMaxEtaVariation->GetXaxis()->GetNbins())
2222     bin = fhMaxEtaVariation->GetXaxis()->GetNbins();
2223
2224   return fhMaxEtaVariation->GetBinContent(bin);
2225 }
2226
2227
2228 //_____________________________________________________________________________
2229 Int_t AliAnalysisTaskPID::GetRandomParticleTypeAccordingToParticleFractions(Double_t trackPt, Double_t jetPt,
2230                                                                             Double_t centralityPercentile,
2231                                                                             Bool_t smearByError,
2232                                                                             Bool_t takeIntoAccountSysError) const
2233 {
2234   // Uses the stored histograms with the particle fractions to generate a random particle type according to these fractions.
2235   // In case of problems (e.g. histo missing), AliPID::kUnknown is returned.
2236   // If smearByError is kTRUE, the used fractions will be random numbers distributed with a gauss with mean
2237   // being the corresponding particle fraction and sigma it's error.
2238   // Note that in this case only the fraction of a random species is varied in this way. The other fractions
2239   // will be re-normalised according their statistical errors.
2240   // The same holds for the systematic error of species "takeIntoAccountSpeciesSysError", but the random number will be
2241   // uniformly distributed within [mean - sys, mean + sys] and the re-normalisation will be weighted with the systematic errors.
2242   // Note that the fractions will be calculated first with only the systematic error taken into account (if desired), including
2243   // re-normalisation. Then, the resulting fractions will be used to calculate the final fractions - either with statistical error
2244   // or without. The species, for which the error will be used for smearing, is the same for sys and stat error.
2245
2246   Double_t prob[AliPID::kSPECIES];
2247   Int_t randomSpecies = (smearByError || takeIntoAccountSysError) ? (Int_t)(fRandom->Rndm() * AliPID::kSPECIES) : -1;
2248   Bool_t success = GetParticleFractions(trackPt, jetPt, centralityPercentile, prob, randomSpecies, randomSpecies);
2249
2250   if (!success)
2251     return AliPID::kUnknown;
2252
2253   Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
2254
2255   if (rnd <= prob[AliPID::kPion])
2256     return AliPID::kPion;
2257   else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon])
2258     return AliPID::kKaon;
2259   else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon] + prob[AliPID::kProton])
2260     return AliPID::kProton;
2261   else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon] + prob[AliPID::kProton] + prob[AliPID::kElectron])
2262     return AliPID::kElectron;
2263
2264   return AliPID::kMuon;  //else it must be a muon (only species left)
2265 }
2266
2267
2268 //_____________________________________________________________________________
2269 AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::GenerateDetectorResponse(AliAnalysisTaskPID::ErrorCode errCode,
2270                                                                            Double_t mean, Double_t sigma,
2271                                                                            Double_t* responses, Int_t nResponses,
2272                                                                            Bool_t usePureGaus)
2273 {
2274   // Generate detector response. If a previous generation was not successful or there is something wrong with this signal generation,
2275   // the function will return kFALSE
2276   if (!responses)
2277     return kError;
2278
2279   // Reset response array
2280   for (Int_t i = 0; i < nResponses; i++)
2281     responses[i] = -999;
2282
2283   if (errCode == kError)
2284     return kError;
2285
2286   ErrorCode ownErrCode = kNoErrors;
2287
2288   if (fUseConvolutedGaus && !usePureGaus) {
2289     // In case of convoluted gauss, calculate the probability density only once to save a lot of time!
2290
2291     TH1* hProbDensity = 0x0;
2292     ownErrCode = SetParamsForConvolutedGaus(mean, sigma);
2293     if (ownErrCode == kError)
2294       return kError;
2295
2296     hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
2297
2298     for (Int_t i = 0; i < nResponses; i++) {
2299       responses[i] = hProbDensity->GetRandom();
2300       //responses[i] fConvolutedGausDeltaPrime->GetRandom(); // MUCH slower than using the binned version via the histogram
2301     }
2302   }
2303   else {
2304     for (Int_t i = 0; i < nResponses; i++) {
2305       responses[i] = fRandom->Gaus(mean, sigma);
2306     }
2307   }
2308
2309   // If forwarded error code was a warning (error case has been handled before), return a warning
2310   if (errCode == kWarning)
2311     return kWarning;
2312
2313   return ownErrCode; // Forward success/warning
2314 }
2315
2316
2317 //_____________________________________________________________________________
2318 void AliAnalysisTaskPID::PrintSettings(Bool_t printSystematicsSettings) const
2319 {
2320   // Print current settings.
2321
2322   printf("\n\nSettings for task %s:\n", GetName());
2323   printf("Is pPb/Pbp: %d -> %s\n", GetIsPbpOrpPb(), GetIsPbpOrpPb() ? "Adapting vertex cuts" : "Using standard vertex cuts");
2324   printf("Track cuts: %s\n", fTrackFilter ? fTrackFilter->GetTitle() : "-");
2325   printf("Eta cut: %.2f <= |eta| <= %.2f\n", GetEtaAbsCutLow(), GetEtaAbsCutUp());
2326   printf("Phi' cut: %d\n", GetUsePhiCut());
2327   printf("TPCCutMIGeo: %d\n", GetUseTPCCutMIGeo());
2328   if (GetUseTPCCutMIGeo()) {
2329     printf("GetCutGeo: %f\n", GetCutGeo());
2330     printf("GetCutNcr: %f\n", GetCutNcr());
2331     printf("GetCutNcl: %f\n", GetCutNcl());
2332   }
2333   printf("TPCnclCut: %d\n", GetUseTPCnclCut());
2334   if (GetUseTPCnclCut()) {
2335     printf("GetCutPureNcl: %d\n", GetCutPureNcl());
2336   }
2337
2338   printf("\n");
2339
2340   printf("Centrality estimator: %s\n", GetCentralityEstimator().Data());
2341
2342   printf("\n");
2343
2344   printf("Pile up rejection type: %d\n", (Int_t)fPileUpRejectionType);
2345
2346   printf("\n");
2347
2348   printf("Use MC-ID for signal generation: %d\n", GetUseMCidForGeneration());
2349   printf("Use ITS: %d\n", GetUseITS());
2350   printf("Use TOF: %d\n", GetUseTOF());
2351   printf("Use priors: %d\n", GetUsePriors());
2352   printf("Use TPC default priors: %d\n", GetUseTPCDefaultPriors());
2353   printf("Use convoluted Gauss: %d\n", GetUseConvolutedGaus());
2354   printf("Accuracy of non-Gaussian tail: %e\n", GetAccuracyNonGaussianTail());
2355   printf("Take into account muons: %d\n", GetTakeIntoAccountMuons());
2356   printf("TOF mode: %d\n", GetTOFmode());
2357   printf("\nParams for transition from gauss to asymmetric shape:\n");
2358   printf("[0]: %e\n", GetConvolutedGaussTransitionPar(0));
2359   printf("[1]: %e\n", GetConvolutedGaussTransitionPar(1));
2360   printf("[2]: %e\n", GetConvolutedGaussTransitionPar(2));
2361
2362   printf("\n");
2363
2364   printf("Do PID: %d\n", fDoPID);
2365   printf("Do Efficiency: %d\n", fDoEfficiency);
2366   printf("Do PtResolution: %d\n", fDoPtResolution);
2367   printf("Do dEdxCheck: %d\n", fDoDeDxCheck);
2368   printf("Do binZeroStudy: %d\n", fDoBinZeroStudy);
2369
2370   printf("\n");
2371
2372   printf("Input from other task: %d\n", GetInputFromOtherTask());
2373   printf("Store additional jet information: %d\n", GetStoreAdditionalJetInformation());
2374   printf("Store centrality percentile: %d", GetStoreCentralityPercentile());
2375
2376   if (printSystematicsSettings)
2377     PrintSystematicsSettings();
2378   else
2379     printf("\n\n\n");
2380 }
2381
2382
2383 //_____________________________________________________________________________
2384 void AliAnalysisTaskPID::PrintSystematicsSettings() const
2385 {
2386   // Print current settings for systematic studies.
2387
2388   printf("\n\nSettings for systematics for task %s:\n", GetName());
2389   printf("SplinesThr:\t%f\n", GetSystematicScalingSplinesThreshold());
2390   printf("SplinesBelowThr:\t%f\n", GetSystematicScalingSplinesBelowThreshold());
2391   printf("SplinesAboveThr:\t%f\n", GetSystematicScalingSplinesAboveThreshold());
2392   printf("EtaCorrMomThr:\t%f\n", GetSystematicScalingEtaCorrectionMomentumThr());
2393   printf("EtaCorrLowP:\t%f\n", GetSystematicScalingEtaCorrectionLowMomenta());
2394   printf("EtaCorrHighP:\t%f\n", GetSystematicScalingEtaCorrectionHighMomenta());
2395   printf("SigmaParaThr:\t%f\n", GetSystematicScalingEtaSigmaParaThreshold());
2396   printf("SigmaParaBelowThr:\t%f\n", GetSystematicScalingEtaSigmaParaBelowThreshold());
2397   printf("SigmaParaAboveThr:\t%f\n", GetSystematicScalingEtaSigmaParaAboveThreshold());
2398   printf("MultCorr:\t%f\n", GetSystematicScalingMultCorrection());
2399   printf("TOF mode: %d\n", GetTOFmode());
2400
2401   printf("\n\n");
2402 }
2403
2404
2405 //_____________________________________________________________________________
2406 Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePDGcode, Double_t centralityPercentile,
2407                                         Double_t jetPt)
2408 {
2409   // Process the track (generate expected response, fill histos, etc.).
2410   // particlePDGcode == 0 means data. Otherwise, the corresponding MC ID will be assumed.
2411
2412   //Printf("Debug: Task %s is starting to process track: dEdx %f, pTPC %f, eta %f, ncl %d\n", GetName(), track->GetTPCsignal(), track->GetTPCmomentum(),
2413   //       track->Eta(), track->GetTPCsignalN());
2414
2415   if(fDebug > 1)
2416     printf("File: %s, Line: %d: ProcessTrack\n", (char*)__FILE__, __LINE__);
2417
2418   if (!fDoPID && !fDoDeDxCheck && !fDoPtResolution)
2419     return kFALSE;
2420
2421   if(fDebug > 2)
2422     printf("File: %s, Line: %d: ProcessTrack -> Processing started\n", (char*)__FILE__, __LINE__);
2423
2424   const Bool_t isMC = (particlePDGcode == 0) ? kFALSE : kTRUE;
2425
2426   Int_t binMC = -1;
2427
2428   if (isMC) {
2429     if (TMath::Abs(particlePDGcode) == 211) {//Pion
2430       binMC = 3;
2431     }
2432     else if (TMath::Abs(particlePDGcode) == 321) {//Kaon
2433       binMC = 1;
2434     }
2435     else if (TMath::Abs(particlePDGcode) == 2212) {//Proton
2436       binMC = 4;
2437     }
2438     else if (TMath::Abs(particlePDGcode) == 11) {//Electron
2439       binMC = 0;
2440     }
2441     else if (TMath::Abs(particlePDGcode) == 13) {//Muon
2442       binMC = 2;
2443     }
2444     else // In MC-ID case, set to underflow bin such that the response from this track is only used for unidentified signal generation
2445          // or signal generation with PID response and the track is still there (as in data) - e.g. relevant w.r.t. deuterons.
2446          // This is important to be as much as possible consistent with data. And the tracks can still be removed by disabling the
2447          // underflow bin for the projections
2448       binMC = -1;
2449   }
2450
2451   // Momenta
2452   //Double_t p = track->GetP();
2453   const Double_t pTPC = track->GetTPCmomentum();
2454   Double_t pT = track->Pt();
2455
2456   Double_t z = -1, xi = -1;
2457   GetJetTrackObservables(pT, jetPt, z, xi);
2458
2459
2460   Double_t trackCharge = track->Charge();
2461
2462   // TPC signal
2463   const Bool_t tuneOnDataTPC = fPIDResponse->IsTunedOnData() &&
2464                                ((fPIDResponse->GetTunedOnDataMask() & AliPIDResponse::kDetTPC) == AliPIDResponse::kDetTPC);
2465   Double_t dEdxTPC = tuneOnDataTPC ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
2466
2467   if (dEdxTPC <= 0) {
2468     if (fDebug > 1)
2469       Printf("Skipping track with strange dEdx value: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), pTPC,
2470              track->Eta(), track->GetTPCsignalN());
2471     return kFALSE;
2472   }
2473
2474   // Completely remove tracks from light nuclei defined by el and pr <dEdx> as function of pT.
2475   // A very loose cut to be sure not to throw away electrons and/or protons
2476   Double_t nSigmaPr = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton);
2477   Double_t nSigmaEl = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
2478
2479   if ((pTPC >= 0.3 && (nSigmaPr > 10 && nSigmaEl > 10)) ||
2480       (pTPC <  0.3 && (nSigmaPr > 15 && nSigmaEl > 15))) {
2481     if (fDebug > 1)
2482       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",
2483              track->GetTPCsignal(), pTPC, pT, track->Eta(), track->GetTPCsignalN(), nSigmaPr, nSigmaEl);
2484     return kFALSE;
2485   }
2486
2487
2488   Double_t dEdxEl, dEdxKa, dEdxPi, dEdxMu, dEdxPr;
2489   Double_t sigmaEl, sigmaKa, sigmaPi, sigmaMu, sigmaPr;
2490
2491   if (fDoAnySystematicStudiesOnTheExpectedSignal) {
2492     // Get the uncorrected signal first and the corresponding correction factors.
2493     // Then modify the correction factors and properly recalculate the corrected dEdx
2494
2495     // Get pure spline values for dEdx_expected, without any correction
2496     dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2497     dEdxKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2498     dEdxPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2499     dEdxMu = !fTakeIntoAccountMuons ? -1 :
2500                             fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2501     dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2502
2503     // Scale splines, if desired
2504     if ((TMath::Abs(fSystematicScalingSplinesBelowThreshold - 1.0) > fgkEpsilon) ||
2505         (TMath::Abs(fSystematicScalingSplinesAboveThreshold - 1.0) > fgkEpsilon)) {
2506
2507       // Tune turn-on of correction for pions (not so relevant for the other species, since at very large momenta)
2508       Double_t bg = 0;
2509       Double_t scaleFactor = 1.;
2510       Double_t usedSystematicScalingSplines = 1.;
2511
2512       bg = pTPC / AliPID::ParticleMass(AliPID::kElectron);
2513       scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2514       usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2515                                             + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2516       dEdxEl *= usedSystematicScalingSplines;
2517
2518       bg = pTPC / AliPID::ParticleMass(AliPID::kKaon);
2519       scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2520       usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2521                                             + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2522       dEdxKa *= usedSystematicScalingSplines;
2523
2524       bg = pTPC / AliPID::ParticleMass(AliPID::kPion);
2525       scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2526       usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2527                                             + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2528       dEdxPi *= usedSystematicScalingSplines;
2529
2530       if (fTakeIntoAccountMuons) {
2531         bg = pTPC / AliPID::ParticleMass(AliPID::kMuon);
2532         scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2533         usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2534                                               + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2535         dEdxMu *= usedSystematicScalingSplines;
2536       }
2537
2538
2539       bg = pTPC / AliPID::ParticleMass(AliPID::kProton);
2540       scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2541       usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2542                                             + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2543       dEdxPr *= usedSystematicScalingSplines;
2544     }
2545
2546     // Get the eta correction factors for the (modified) expected dEdx
2547     Double_t etaCorrEl = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxEl) : 1.;
2548     Double_t etaCorrKa = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxKa) : 1.;
2549     Double_t etaCorrPi = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPi) : 1.;
2550     Double_t etaCorrMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCEtaCorrection() ?
2551                             fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxMu) : 1.;
2552     Double_t etaCorrPr = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPr) : 1.;
2553
2554     // Scale eta correction factors, if desired (and eta correction maps are to be used, otherwise it is not possible!)
2555     if (fPIDResponse->UseTPCEtaCorrection() &&
2556         (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon ||
2557          TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon)) {
2558       // Since we do not want to scale the splines with this, but only the eta variation, only scale the deviation of the correction factor!
2559       // 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!
2560
2561
2562       // Due to additional azimuthal effects, there is an additional eta dependence for low momenta which is not corrected successfully so far.
2563       // One can assign a different (higher) systematic scale factor for this low-p region and a threshold which separates low- and high-p.
2564       // 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
2565       Double_t usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionHighMomenta;
2566
2567       if (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - fSystematicScalingEtaCorrectionLowMomenta) > fgkEpsilon) {
2568         const Double_t fractionHighMomentumScaleFactor = 0.5 * (1. + TMath::Erf((pTPC - fSystematicScalingEtaCorrectionMomentumThr) / 0.1));
2569         usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionLowMomenta * (1 - fractionHighMomentumScaleFactor)
2570                                              + fSystematicScalingEtaCorrectionHighMomenta * fractionHighMomentumScaleFactor;
2571       }
2572
2573       Double_t maxEtaVariationEl = GetMaxEtaVariation(dEdxEl);
2574       etaCorrEl = etaCorrEl * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrEl - 1.0) / maxEtaVariationEl);
2575
2576       Double_t maxEtaVariationKa = GetMaxEtaVariation(dEdxKa);
2577       etaCorrKa = etaCorrKa * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrKa - 1.0) / maxEtaVariationKa);
2578
2579       Double_t maxEtaVariationPi = GetMaxEtaVariation(dEdxPi);
2580       etaCorrPi = etaCorrPi * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrPi - 1.0) / maxEtaVariationPi);
2581
2582       if (fTakeIntoAccountMuons) {
2583         Double_t maxEtaVariationMu = GetMaxEtaVariation(dEdxMu);
2584         etaCorrMu = etaCorrMu * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrMu - 1.0) / maxEtaVariationMu);
2585       }
2586       else
2587         etaCorrMu = 1.0;
2588
2589       Double_t maxEtaVariationPr = GetMaxEtaVariation(dEdxPr);
2590       etaCorrPr = etaCorrPr * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrPr - 1.0) / maxEtaVariationPr);
2591
2592
2593       /*OLD
2594       etaCorrEl = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrEl - 1.0);
2595       etaCorrKa = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrKa - 1.0);
2596       etaCorrPi = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPi - 1.0);
2597       etaCorrMu = fTakeIntoAccountMuons ? (1.0 + usedSystematicScalingEtaCorrection * (etaCorrMu - 1.0)) : 1.0;
2598       etaCorrPr = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPr - 1.0);
2599       */
2600     }
2601
2602     // Get the multiplicity correction factors for the (modified) expected dEdx
2603     const Int_t currEvtMultiplicity = fPIDResponse->GetTPCResponse().GetCurrentEventMultiplicity();
2604
2605     Double_t multiplicityCorrEl = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2606                                     dEdxEl * etaCorrEl, currEvtMultiplicity) : 1.;
2607     Double_t multiplicityCorrKa = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2608                                     dEdxKa * etaCorrKa, currEvtMultiplicity) : 1.;
2609     Double_t multiplicityCorrPi = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2610                                     dEdxPi * etaCorrPi, currEvtMultiplicity) : 1.;
2611     Double_t multiplicityCorrMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCMultiplicityCorrection() ?
2612                                     fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track, dEdxMu * etaCorrMu, currEvtMultiplicity) : 1.;
2613     Double_t multiplicityCorrPr = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2614                                     dEdxPr * etaCorrPr, currEvtMultiplicity) : 1.;
2615
2616     Double_t multiplicityCorrSigmaEl = fPIDResponse->UseTPCMultiplicityCorrection() ?
2617                                         fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxEl * etaCorrEl, currEvtMultiplicity) : 1.;
2618     Double_t multiplicityCorrSigmaKa = fPIDResponse->UseTPCMultiplicityCorrection() ?
2619                                         fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxKa * etaCorrKa, currEvtMultiplicity) : 1.;
2620     Double_t multiplicityCorrSigmaPi = fPIDResponse->UseTPCMultiplicityCorrection() ?
2621                                         fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxPi * etaCorrPi, currEvtMultiplicity) : 1.;
2622     Double_t multiplicityCorrSigmaMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCMultiplicityCorrection() ?
2623                                         fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxMu * etaCorrMu, currEvtMultiplicity) : 1.;
2624     Double_t multiplicityCorrSigmaPr = fPIDResponse->UseTPCMultiplicityCorrection() ?
2625                                         fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxPr * etaCorrPr, currEvtMultiplicity) : 1.;
2626
2627     // Scale multiplicity correction factors, if desired (and multiplicity correction functions are to be used, otherwise it is not possible!)
2628     if (fPIDResponse->UseTPCMultiplicityCorrection() && TMath::Abs(fSystematicScalingMultCorrection - 1.0) > fgkEpsilon) {
2629       // Since we do not want to scale the splines with this, but only the multiplicity variation, only scale the deviation of the correction factor!
2630       // 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!
2631
2632       multiplicityCorrEl = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrEl - 1.0);
2633       multiplicityCorrKa = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrKa - 1.0);
2634       multiplicityCorrPi = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrPi - 1.0);
2635       multiplicityCorrMu = fTakeIntoAccountMuons ? (1.0 + fSystematicScalingMultCorrection * (multiplicityCorrMu - 1.0)) : 1.0;
2636       multiplicityCorrPr = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrPr - 1.0);
2637     }
2638
2639     // eta correction must be enabled in order to use the new sigma parametrisation maps. Since this is the absolute sigma
2640     // for a track calculated with the unscaled paramaters, we have to devide out dEdxExpectedEtaCorrected and then need
2641     // to scale with the multiplicitySigmaCorrFactor * fSystematicScalingEtaSigmaPara. In the end, one has to scale with the
2642     // (modified) dEdx to get the absolute sigma
2643     // This means there is no extra parameter for the multiplicitySigmaCorrFactor, but only for the sigma map itself.
2644     // This is valid, since it appears only as a product. One has to assume a larger systematic shift in case of additional
2645     // multiplicity dependence....
2646
2647     Bool_t doSigmaSystematics = (TMath::Abs(fSystematicScalingEtaSigmaParaBelowThreshold - 1.0) > fgkEpsilon) ||
2648                                 (TMath::Abs(fSystematicScalingEtaSigmaParaAboveThreshold - 1.0) > fgkEpsilon);
2649
2650
2651     Double_t dEdxElExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2652                                                                                kTRUE, kFALSE);
2653     Double_t systematicScalingEtaSigmaParaEl = 1.;
2654     if (doSigmaSystematics) {
2655       Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxElExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2656       systematicScalingEtaSigmaParaEl = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2657                                         + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2658     }
2659     Double_t sigmaRelEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2660                                                                           kTRUE, kFALSE)
2661                           / dEdxElExpected * systematicScalingEtaSigmaParaEl * multiplicityCorrSigmaEl;
2662
2663
2664     Double_t dEdxKaExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2665                                                                                kTRUE, kFALSE);
2666     Double_t systematicScalingEtaSigmaParaKa = 1.;
2667     if (doSigmaSystematics) {
2668       Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxKaExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2669       systematicScalingEtaSigmaParaKa = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2670                                         + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2671     }
2672     Double_t sigmaRelKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2673                                                                           kTRUE, kFALSE)
2674                           / dEdxKaExpected * systematicScalingEtaSigmaParaKa * multiplicityCorrSigmaKa;
2675
2676
2677     Double_t dEdxPiExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2678                                                                                kTRUE, kFALSE);
2679     Double_t systematicScalingEtaSigmaParaPi = 1.;
2680     if (doSigmaSystematics) {
2681       Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxPiExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2682       systematicScalingEtaSigmaParaPi = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2683                                         + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2684     }
2685     Double_t sigmaRelPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2686                                                                           kTRUE, kFALSE)
2687                           / dEdxPiExpected * systematicScalingEtaSigmaParaPi * multiplicityCorrSigmaPi;
2688
2689
2690     Double_t sigmaRelMu = 999.;
2691     if (fTakeIntoAccountMuons) {
2692       Double_t dEdxMuExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2693                                                                                  kTRUE, kFALSE);
2694       Double_t systematicScalingEtaSigmaParaMu = 1.;
2695       if (doSigmaSystematics) {
2696         Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxMuExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2697         systematicScalingEtaSigmaParaMu = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2698                                           + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2699       }
2700       sigmaRelMu = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2701                    / dEdxMuExpected * systematicScalingEtaSigmaParaMu * multiplicityCorrSigmaMu;
2702     }
2703
2704
2705     Double_t dEdxPrExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2706                                                                                kTRUE, kFALSE);
2707     Double_t systematicScalingEtaSigmaParaPr = 1.;
2708     if (doSigmaSystematics) {
2709       Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxPrExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2710       systematicScalingEtaSigmaParaPr = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2711                                         + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2712     }
2713     Double_t sigmaRelPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2714                                                                           kTRUE, kFALSE)
2715                           / dEdxPrExpected * systematicScalingEtaSigmaParaPr * multiplicityCorrSigmaPr;
2716
2717     // Now scale the (possibly modified) spline values with the (possibly modified) correction factors
2718     dEdxEl *= etaCorrEl * multiplicityCorrEl;
2719     dEdxKa *= etaCorrKa * multiplicityCorrKa;
2720     dEdxPi *= etaCorrPi * multiplicityCorrPi;
2721     dEdxMu *= etaCorrMu * multiplicityCorrMu;
2722     dEdxPr *= etaCorrPr * multiplicityCorrPr;
2723
2724     // Finally, get the absolute sigma
2725     sigmaEl = sigmaRelEl * dEdxEl;
2726     sigmaKa = sigmaRelKa * dEdxKa;
2727     sigmaPi = sigmaRelPi * dEdxPi;
2728     sigmaMu = sigmaRelMu * dEdxMu;
2729     sigmaPr = sigmaRelPr * dEdxPr;
2730
2731   }
2732   else {
2733     // No systematic studies on expected signal - just take it as it comes from the TPCPIDResponse
2734     dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2735                                                               fPIDResponse->UseTPCEtaCorrection(),
2736                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2737     dEdxKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2738                                                               fPIDResponse->UseTPCEtaCorrection(),
2739                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2740     dEdxPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2741                                                               fPIDResponse->UseTPCEtaCorrection(),
2742                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2743     dEdxMu = !fTakeIntoAccountMuons ? -1 :
2744                             fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2745                                                                              fPIDResponse->UseTPCEtaCorrection(),
2746                                                                              fPIDResponse->UseTPCMultiplicityCorrection());
2747     dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2748                                                               fPIDResponse->UseTPCEtaCorrection(),
2749                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2750
2751     sigmaEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2752                                                               fPIDResponse->UseTPCEtaCorrection(),
2753                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2754     sigmaKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2755                                                               fPIDResponse->UseTPCEtaCorrection(),
2756                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2757     sigmaPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2758                                                               fPIDResponse->UseTPCEtaCorrection(),
2759                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2760     sigmaMu = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2761                                                               fPIDResponse->UseTPCEtaCorrection(),
2762                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2763     sigmaPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2764                                                               fPIDResponse->UseTPCEtaCorrection(),
2765                                                               fPIDResponse->UseTPCMultiplicityCorrection());
2766   }
2767
2768   Double_t deltaPrimeElectron = (dEdxEl > 0) ? dEdxTPC / dEdxEl : -1;
2769   if (dEdxEl <= 0)  {
2770     Printf("Error: Expected TPC signal <= 0 for electron hypothesis");
2771     return kFALSE;
2772   }
2773
2774   Double_t deltaPrimeKaon = (dEdxKa > 0) ? dEdxTPC / dEdxKa : -1;
2775   if (dEdxKa <= 0)  {
2776     Printf("Error: Expected TPC signal <= 0 for kaon hypothesis");
2777     return kFALSE;
2778   }
2779
2780   Double_t deltaPrimePion = (dEdxPi > 0) ? dEdxTPC / dEdxPi : -1;
2781   if (dEdxPi <= 0)  {
2782     Printf("Error: Expected TPC signal <= 0 for pion hypothesis");
2783     return kFALSE;
2784   }
2785
2786   Double_t deltaPrimeProton = (dEdxPr > 0) ? dEdxTPC / dEdxPr : -1;
2787   if (dEdxPr <= 0)  {
2788     Printf("Error: Expected TPC signal <= 0 for proton hypothesis");
2789     return kFALSE;
2790   }
2791
2792   if(fDebug > 2)
2793     printf("File: %s, Line: %d: ProcessTrack -> Compute probabilities\n", (char*)__FILE__, __LINE__);
2794
2795
2796   // Use probabilities to weigh the response generation for the different species.
2797   // Also determine the most probable particle type.
2798   Double_t prob[AliPID::kSPECIESC];
2799   for (Int_t i = 0; i < AliPID::kSPECIESC; i++)
2800     prob[i] = 0;
2801
2802   fPIDcombined->ComputeProbabilities(track, fPIDResponse, prob);
2803
2804   // If muons are not to be taken into account, just set their probability to zero and normalise the remaining probabilities later
2805   if (!fTakeIntoAccountMuons)
2806     prob[AliPID::kMuon] = 0;
2807
2808   // Priors might cause problems in case of mismatches, i.e. mismatch particles will be identified as pions.
2809   // Can be easily caught for low p_TPC by just setting pion (and muon) probs to zero, if dEdx is larger than
2810   // expected dEdx of electrons and kaons
2811   if (pTPC < 1. && dEdxTPC > dEdxEl && dEdxTPC > dEdxKa) {
2812     prob[AliPID::kMuon] = 0;
2813     prob[AliPID::kPion] = 0;
2814   }
2815
2816
2817   // Bug: One can set the number of species for PIDcombined, but PIDcombined will call PIDresponse, which writes without testing
2818   // the probs for kSPECIESC (including light nuclei) into the array.
2819   // In this case, when only kSPECIES are considered, the probabilities have to be rescaled!
2820
2821   // Exceptions:
2822   // 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
2823   //    high-pT region (also contribution there completely negligable!)
2824   // 2. If dEdx is larger than the expected dEdx + 5 sigma of el and pr, it is most likely a light nucleus
2825   //    (if dEdx is more than 5 sigma away from expectation, flat TPC probs are set... and the light nuclei splines are not
2826   //     too accurate)
2827   // In these cases:
2828   // The highest expected dEdx is either for pr or for el. Assign 100% probability to the species with highest expected dEdx then.
2829   // In all other cases: Throw away light nuclei probs and rescale other probs to 100%.
2830
2831   // Find most probable species for the ID
2832   Int_t maxIndex = TMath::LocMax(AliPID::kSPECIESC, prob);
2833
2834   if ((prob[maxIndex] > 0 && maxIndex >= AliPID::kSPECIES && dEdxTPC > dEdxEl + 5. * sigmaEl) ||
2835       (dEdxTPC > dEdxEl + 5. * sigmaEl && dEdxTPC > dEdxPr + 5. * sigmaPr)) {
2836     for (Int_t i = 0; i < AliPID::kSPECIESC; i++)
2837       prob[i] = 0;
2838
2839     if (dEdxEl > dEdxPr) {
2840       prob[AliPID::kElectron] = 1.;
2841       maxIndex = AliPID::kElectron;
2842     }
2843     else {
2844       prob[AliPID::kProton] = 1.;
2845       maxIndex = AliPID::kProton;
2846     }
2847   }
2848   else {
2849     for (Int_t i = AliPID::kSPECIES; i < AliPID::kSPECIESC; i++)
2850       prob[i] = 0;
2851
2852     Double_t probSum = 0.;
2853     for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2854       probSum += prob[i];
2855
2856     if (probSum > 0) {
2857       for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2858         prob[i] /= probSum;
2859     }
2860
2861     // If all probabilities are equal (should only happen if no priors are used), set most probable PID to pion
2862     // because LocMax returns just the first maximal value (i.e. AliPID::kElectron)
2863     if (TMath::Abs(prob[maxIndex] - 1. / AliPID::kSPECIES) < 1e-5)
2864       maxIndex = AliPID::kPion;
2865
2866     // In all other cases, the maximum remains untouched from the scaling (and is < AliPID::kSPECIES by construction)
2867   }
2868
2869   if (fDoDeDxCheck) {
2870     // Generate the expected responses in DeltaPrime and translate to real dEdx. Only take responses for exactly that species
2871     // (i.e. with pre-PID)
2872
2873     Int_t numGenEntries = fgkMaxNumGenEntries; // fgkMaxNumGenEntries = 500
2874
2875     ErrorCode errCode = kNoErrors;
2876
2877     if(fDebug > 2)
2878       printf("File: %s, Line: %d: ProcessTrack -> Generate Responses for dEdx check\n", (char*)__FILE__, __LINE__);
2879
2880     // Electrons
2881     errCode = GenerateDetectorResponse(errCode, 1., sigmaEl / dEdxEl, fGenRespElDeltaPrimeEl, numGenEntries);
2882
2883     // Kaons
2884     errCode = GenerateDetectorResponse(errCode, 1., sigmaKa / dEdxKa, fGenRespKaDeltaPrimeKa, numGenEntries);
2885
2886     // Pions
2887     errCode = GenerateDetectorResponse(errCode, 1., sigmaPi / dEdxPi, fGenRespPiDeltaPrimePi, numGenEntries);
2888
2889     // Protons
2890     errCode = GenerateDetectorResponse(errCode, 1., sigmaPr / dEdxPr, fGenRespPrDeltaPrimePr, numGenEntries);
2891
2892     if (errCode == kNoErrors) {
2893       Double_t genEntry[kDeDxCheckNumAxes];
2894       genEntry[kDeDxCheckJetPt] = jetPt;
2895       genEntry[kDeDxCheckEtaAbs] = TMath::Abs(track->Eta());
2896       genEntry[kDeDxCheckP] = pTPC;
2897
2898
2899       for (Int_t n = 0; n < numGenEntries; n++)  {
2900         // If no MC info is available or shall not be used, use weighting with priors to generate entries for the different species
2901         Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
2902
2903         // Consider generated response as originating from...
2904         if (rnd <= prob[AliPID::kElectron]) {
2905           genEntry[kDeDxCheckPID] = 0; // ... an electron
2906           genEntry[kDeDxCheckDeDx] = fGenRespElDeltaPrimeEl[n] * dEdxEl;
2907         }
2908         else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon]) {
2909           genEntry[kDeDxCheckPID] = 1;  // ... a kaon
2910           genEntry[kDeDxCheckDeDx] = fGenRespKaDeltaPrimeKa[n] * dEdxKa;
2911         }
2912         else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon] + prob[AliPID::kPion]) {
2913           genEntry[kDeDxCheckPID] = 2;  // ... a pion (or a muon)
2914           genEntry[kDeDxCheckDeDx] = fGenRespPiDeltaPrimePi[n] * dEdxPi;
2915         }
2916         else {
2917           genEntry[kDeDxCheckPID] = 3;  // ... a proton
2918           genEntry[kDeDxCheckDeDx] = fGenRespPrDeltaPrimePr[n] * dEdxPr;
2919         }
2920
2921         fDeDxCheck->Fill(genEntry);
2922       }
2923     }
2924
2925     if(fDebug > 2)
2926       printf("File: %s, Line: %d: ProcessTrack -> Generate Responses for dEdx check done\n", (char*)__FILE__, __LINE__);
2927   }
2928
2929   if (fDoPtResolution) {
2930     // Check shared clusters, which is done together with the pT resolution
2931     Double_t qaEntry[kQASharedClsNumAxes];
2932     qaEntry[kQASharedClsJetPt] = jetPt;
2933     qaEntry[kQASharedClsPt] = pT;
2934     qaEntry[kDeDxCheckP] = pTPC;
2935     qaEntry[kQASharedClsNumSharedCls] = track->GetTPCSharedMapPtr()->CountBits();
2936
2937     Int_t iRowInd = -1;
2938     // iRowInd == -1 for "all rows w/o multiple counting"
2939     qaEntry[kQASharedClsPadRow] = iRowInd;
2940     fQASharedCls->Fill(qaEntry);
2941
2942     // Fill hist for every pad row with shared cluster
2943     for (iRowInd = 0; iRowInd < 159; iRowInd++) {
2944       if (track->GetTPCSharedMapPtr()->TestBitNumber(iRowInd)) {
2945         qaEntry[kQASharedClsPadRow] = iRowInd;
2946         fQASharedCls->Fill(qaEntry);
2947       }
2948     }
2949   }
2950
2951   if (!fDoPID)
2952     return kTRUE;
2953
2954   if (!isMC) {
2955     // Clean up the most probable PID for low momenta (not an issue for the probabilities, since fractions small,
2956     // but to get clean (and nice looking) templates (if most probable PID is used instead of generated responses), it should be done).
2957     // Problem: If more than 5 sigma away (and since priors also included), most probable PID for dEdx around protons could be pions.
2958     // Idea: Only clean selection for K and p possible. So, just take for the most probable PID the probs from TPC only calculated
2959     // by hand without restriction to 5 sigma and without priors. This will not affect the template generation, but afterwards one
2960     // can replace the K and p templates with the most probable PID distributions, so all other species templates remain the same.
2961     // I.e. it doesn't hurt if the most probable PID is distributed incorrectly between el, pi, mu.
2962     Bool_t maxIndexSetManually = kFALSE;
2963
2964     if (pTPC < 1.) {
2965       Double_t probManualTPC[AliPID::kSPECIES];
2966        for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2967         probManualTPC[i] = 0;
2968
2969       probManualTPC[AliPID::kElectron] = TMath::Exp(-0.5*(dEdxTPC-dEdxEl)*(dEdxTPC-dEdxEl)/(sigmaEl*sigmaEl))/sigmaEl;
2970       // Muons are not important here, just ignore and save processing time
2971       probManualTPC[AliPID::kMuon] = 0.;//TMath::Exp(-0.5*(dEdxTPC-dEdxMu)*(dEdxTPC-dEdxMu)/(sigmaMu*sigmaMu))/sigmaMu;
2972       probManualTPC[AliPID::kPion] = TMath::Exp(-0.5*(dEdxTPC-dEdxPi)*(dEdxTPC-dEdxPi)/(sigmaPi*sigmaPi))/sigmaPi;
2973       probManualTPC[AliPID::kKaon] = TMath::Exp(-0.5*(dEdxTPC-dEdxKa)*(dEdxTPC-dEdxKa)/(sigmaKa*sigmaKa))/sigmaKa;
2974       probManualTPC[AliPID::kProton] = TMath::Exp(-0.5*(dEdxTPC-dEdxPr)*(dEdxTPC-dEdxPr)/(sigmaPr*sigmaPr))/sigmaPr;
2975
2976       const Int_t maxIndexManualTPC = TMath::LocMax(AliPID::kSPECIES, probManualTPC);
2977       // Should always be > 0, but in case the dEdx is far off such that the machine accuracy sets every prob to zero,
2978       // better take the "old" result
2979       if (probManualTPC[maxIndexManualTPC] > 0.)
2980         maxIndex = maxIndexManualTPC;
2981
2982       maxIndexSetManually = kTRUE;
2983     }
2984
2985
2986     // Translate from AliPID numbering to numbering of this class
2987     if (prob[maxIndex] > 0 || maxIndexSetManually) {
2988       if (maxIndex == AliPID::kElectron)
2989         binMC = 0;
2990       else if (maxIndex == AliPID::kKaon)
2991         binMC = 1;
2992       else if (maxIndex == AliPID::kMuon)
2993         binMC = 2;
2994       else if (maxIndex == AliPID::kPion)
2995         binMC = 3;
2996       else if (maxIndex == AliPID::kProton)
2997         binMC = 4;
2998       else
2999         binMC = -1;
3000     }
3001     else {
3002       // Only take track into account for expectation values, if valid pid response is available.. Otherwise: Set to underflow bin.
3003       binMC = -1;
3004     }
3005   }
3006
3007   /*
3008   //For testing: Swap p<->pT to analyse pure p-dependence => Needs to be removed later
3009   Double_t temp = pT;
3010   pT = pTPC;
3011   pTPC = temp;
3012   */
3013
3014   TOFpidInfo tofPIDinfo = GetTOFType(track, fTOFmode);
3015
3016   Double_t entry[fStoreAdditionalJetInformation ? kDataNumAxes : kDataNumAxes - fgkNumJetAxes];
3017   entry[kDataMCID] = binMC;
3018   entry[kDataSelectSpecies] = 0;
3019   entry[kDataPt] = pT;
3020   entry[kDataDeltaPrimeSpecies] = deltaPrimeElectron;
3021   entry[kDataCentrality] = centralityPercentile;
3022
3023   if (fStoreAdditionalJetInformation) {
3024     entry[kDataJetPt] = jetPt;
3025     entry[kDataZ] = z;
3026     entry[kDataXi] = xi;
3027   }
3028
3029   entry[GetIndexOfChargeAxisData()] = trackCharge;
3030   entry[GetIndexOfTOFpidInfoAxisData()] = tofPIDinfo;
3031
3032   fhPIDdataAll->Fill(entry);
3033
3034   entry[kDataSelectSpecies] = 1;
3035   entry[kDataDeltaPrimeSpecies] = deltaPrimeKaon;
3036   fhPIDdataAll->Fill(entry);
3037
3038   entry[kDataSelectSpecies] = 2;
3039   entry[kDataDeltaPrimeSpecies] = deltaPrimePion;
3040   fhPIDdataAll->Fill(entry);
3041
3042   entry[kDataSelectSpecies] = 3;
3043   entry[kDataDeltaPrimeSpecies] = deltaPrimeProton;
3044   fhPIDdataAll->Fill(entry);
3045
3046
3047   // Construct the expected shape for the transition p -> pT
3048
3049   Double_t genEntry[fStoreAdditionalJetInformation ? kGenNumAxes : kGenNumAxes - fgkNumJetAxes];
3050   genEntry[kGenMCID] = binMC;
3051   genEntry[kGenSelectSpecies] = 0;
3052   genEntry[kGenPt] = pT;
3053   genEntry[kGenDeltaPrimeSpecies] = -999;
3054   genEntry[kGenCentrality] = centralityPercentile;
3055
3056   if (fStoreAdditionalJetInformation) {
3057     genEntry[kGenJetPt] = jetPt;
3058     genEntry[kGenZ] = z;
3059     genEntry[kGenXi] = xi;
3060   }
3061
3062   genEntry[GetIndexOfChargeAxisGen()] = trackCharge;
3063   genEntry[GetIndexOfTOFpidInfoAxisGen()] = tofPIDinfo;
3064
3065   // Generate numGenEntries "responses" with fluctuations around the expected values.
3066   // fgkMaxNumGenEntries = 500 turned out to give reasonable templates even for highest track pT in 15-20 GeV/c jet pT bin.
3067   Int_t numGenEntries = fgkMaxNumGenEntries; // fgkMaxNumGenEntries = 500
3068
3069   /*OLD: Different number of responses depending on pT and jet pT for fgkMaxNumGenEntries = 1000
3070    * => Problem: If threshold to higher number of responses inside a bin (or after rebinning), then the template
3071    * is biased to the higher pT.
3072   // Generate numGenEntries "responses" with fluctuations around the expected values.
3073   // The higher the (transverse) momentum, the more "responses" will be generated in order not to run out of statistics too fast.
3074   Int_t numGenEntries = 10;
3075
3076   // Jets have even worse statistics, therefore, scale numGenEntries further
3077   if (jetPt >= 40)
3078     numGenEntries *= 20;
3079   else if (jetPt >= 20)
3080     numGenEntries *= 10;
3081   else if (jetPt >= 10)
3082     numGenEntries *= 2;
3083
3084
3085
3086   // Do not generate more entries than available in memory!
3087   if (numGenEntries > fgkMaxNumGenEntries)// fgkMaxNumGenEntries = 1000
3088     numGenEntries = fgkMaxNumGenEntries;
3089   */
3090
3091
3092   ErrorCode errCode = kNoErrors;
3093
3094   if(fDebug > 2)
3095     printf("File: %s, Line: %d: ProcessTrack -> Generate Responses\n", (char*)__FILE__, __LINE__);
3096
3097   // Electrons
3098   errCode = GenerateDetectorResponse(errCode, 1.,              sigmaEl / dEdxEl, fGenRespElDeltaPrimeEl, numGenEntries);
3099   errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxKa, sigmaEl / dEdxKa, fGenRespElDeltaPrimeKa, numGenEntries);
3100   errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxPi, sigmaEl / dEdxPi, fGenRespElDeltaPrimePi, numGenEntries);
3101   errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxPr, sigmaEl / dEdxPr, fGenRespElDeltaPrimePr, numGenEntries);
3102
3103   // Kaons
3104   errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxEl, sigmaKa / dEdxEl, fGenRespKaDeltaPrimeEl, numGenEntries);
3105   errCode = GenerateDetectorResponse(errCode, 1.,              sigmaKa / dEdxKa, fGenRespKaDeltaPrimeKa, numGenEntries);
3106   errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxPi, sigmaKa / dEdxPi, fGenRespKaDeltaPrimePi, numGenEntries);
3107   errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxPr, sigmaKa / dEdxPr, fGenRespKaDeltaPrimePr, numGenEntries);
3108
3109   // Pions
3110   errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxEl, sigmaPi / dEdxEl, fGenRespPiDeltaPrimeEl, numGenEntries);
3111   errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxKa, sigmaPi / dEdxKa, fGenRespPiDeltaPrimeKa, numGenEntries);
3112   errCode = GenerateDetectorResponse(errCode, 1.,              sigmaPi / dEdxPi, fGenRespPiDeltaPrimePi, numGenEntries);
3113   errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxPr, sigmaPi / dEdxPr, fGenRespPiDeltaPrimePr, numGenEntries);
3114
3115   // Muons, if desired
3116   if (fTakeIntoAccountMuons) {
3117     errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxEl, sigmaMu / dEdxEl, fGenRespMuDeltaPrimeEl, numGenEntries);
3118     errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxKa, sigmaMu / dEdxKa, fGenRespMuDeltaPrimeKa, numGenEntries);
3119     errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxPi, sigmaMu / dEdxPi, fGenRespMuDeltaPrimePi, numGenEntries);
3120     errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxPr, sigmaMu / dEdxPr, fGenRespMuDeltaPrimePr, numGenEntries);
3121   }
3122
3123   // Protons
3124   errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxEl, sigmaPr / dEdxEl, fGenRespPrDeltaPrimeEl, numGenEntries);
3125   errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxKa, sigmaPr / dEdxKa, fGenRespPrDeltaPrimeKa, numGenEntries);
3126   errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxPi, sigmaPr / dEdxPi, fGenRespPrDeltaPrimePi, numGenEntries);
3127   errCode = GenerateDetectorResponse(errCode, 1.,              sigmaPr / dEdxPr, fGenRespPrDeltaPrimePr, numGenEntries);
3128
3129   if (errCode != kNoErrors) {
3130     if (errCode == kWarning && fDebug > 1) {
3131       Printf("Warning: Questionable detector response for track, most likely due to very low number of PID clusters! Debug output (dEdx_expected, sigma_expected):");
3132     }
3133     else
3134       Printf("Error: Failed to generate detector response for track - skipped! Debug output (dEdx_expected, sigma_expected):");
3135
3136     if (fDebug > 1) {
3137       Printf("Pr: %e, %e", dEdxPr, sigmaPr);
3138       Printf("Pi: %e, %e", dEdxPi, sigmaPi);
3139       Printf("El: %e, %e", dEdxEl, sigmaEl);
3140       Printf("Mu: %e, %e", dEdxMu, sigmaMu);
3141       Printf("Ka: %e, %e", dEdxKa, sigmaKa);
3142       Printf("track: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), track->GetTPCmomentum(), track->Eta(),
3143             track->GetTPCsignalN());
3144     }
3145
3146     if (errCode != kWarning) {
3147       return kFALSE;// Skip generated response in case of error
3148     }
3149   }
3150
3151   for (Int_t n = 0; n < numGenEntries; n++)  {
3152     if (!isMC || !fUseMCidForGeneration) {
3153       // If no MC info is available or shall not be used, use weighting with priors to generate entries for the different species
3154       Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
3155
3156       // Consider generated response as originating from...
3157       if (rnd <= prob[AliPID::kElectron])
3158         genEntry[kGenMCID] = 0; // ... an electron
3159       else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon])
3160         genEntry[kGenMCID] = 1;  // ... a kaon
3161       else if (rnd <=  prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon])
3162         genEntry[kGenMCID] = 2;  // ... a muon -> NOTE: prob[AliPID::kMuon] = 0 in case of fTakeIntoAccountMuons = kFALSE
3163       else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon] + prob[AliPID::kPion])
3164         genEntry[kGenMCID] = 3;  // ... a pion
3165       else
3166         genEntry[kGenMCID] = 4;  // ... a proton
3167     }
3168
3169     // Electrons
3170     genEntry[kGenSelectSpecies] = 0;
3171     genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeEl[n];
3172     fhGenEl->Fill(genEntry);
3173
3174     genEntry[kGenSelectSpecies] = 1;
3175     genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeKa[n];
3176     fhGenEl->Fill(genEntry);
3177
3178     genEntry[kGenSelectSpecies] = 2;
3179     genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePi[n];
3180     fhGenEl->Fill(genEntry);
3181
3182     genEntry[kGenSelectSpecies] = 3;
3183     genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePr[n];
3184     fhGenEl->Fill(genEntry);
3185
3186     // Kaons
3187     genEntry[kGenSelectSpecies] = 0;
3188     genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeEl[n];
3189     fhGenKa->Fill(genEntry);
3190
3191     genEntry[kGenSelectSpecies] = 1;
3192     genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeKa[n];
3193     fhGenKa->Fill(genEntry);
3194
3195     genEntry[kGenSelectSpecies] = 2;
3196     genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePi[n];
3197     fhGenKa->Fill(genEntry);
3198
3199     genEntry[kGenSelectSpecies] = 3;
3200     genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePr[n];
3201     fhGenKa->Fill(genEntry);
3202
3203     // Pions
3204     genEntry[kGenSelectSpecies] = 0;
3205     genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeEl[n];
3206     fhGenPi->Fill(genEntry);
3207
3208     genEntry[kGenSelectSpecies] = 1;
3209     genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeKa[n];
3210     fhGenPi->Fill(genEntry);
3211
3212     genEntry[kGenSelectSpecies] = 2;
3213     genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePi[n];
3214     fhGenPi->Fill(genEntry);
3215
3216     genEntry[kGenSelectSpecies] = 3;
3217     genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePr[n];
3218     fhGenPi->Fill(genEntry);
3219
3220     if (fTakeIntoAccountMuons) {
3221       // Muons, if desired
3222       genEntry[kGenSelectSpecies] = 0;
3223       genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeEl[n];
3224       fhGenMu->Fill(genEntry);
3225
3226       genEntry[kGenSelectSpecies] = 1;
3227       genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeKa[n];
3228       fhGenMu->Fill(genEntry);
3229
3230       genEntry[kGenSelectSpecies] = 2;
3231       genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePi[n];
3232       fhGenMu->Fill(genEntry);
3233
3234       genEntry[kGenSelectSpecies] = 3;
3235       genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePr[n];
3236       fhGenMu->Fill(genEntry);
3237     }
3238
3239     // Protons
3240     genEntry[kGenSelectSpecies] = 0;
3241     genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeEl[n];
3242     fhGenPr->Fill(genEntry);
3243
3244     genEntry[kGenSelectSpecies] = 1;
3245     genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeKa[n];
3246     fhGenPr->Fill(genEntry);
3247
3248     genEntry[kGenSelectSpecies] = 2;
3249     genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePi[n];
3250     fhGenPr->Fill(genEntry);
3251
3252     genEntry[kGenSelectSpecies] = 3;
3253     genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePr[n];
3254     fhGenPr->Fill(genEntry);
3255   }
3256
3257   if(fDebug > 2)
3258     printf("File: %s, Line: %d: ProcessTrack -> Done\n", (char*)__FILE__, __LINE__);
3259
3260   return kTRUE;
3261 }
3262
3263
3264 //_____________________________________________________________________________
3265 Bool_t AliAnalysisTaskPID::SetConvolutedGaussLambdaParameter(Double_t lambda)
3266 {
3267   // Set the lambda parameter of the convolution to the desired value. Automatically
3268   // calculates the parameters for the transition (restricted) gauss -> convoluted gauss.
3269   fConvolutedGaussTransitionPars[2] = lambda;
3270
3271   // Save old parameters and settings of function to restore them later:
3272   Double_t* oldFuncParams = new Double_t[fkConvolutedGausNPar];
3273   for (Int_t i = 0; i < fkConvolutedGausNPar; i++)
3274     oldFuncParams[i] = fConvolutedGausDeltaPrime->GetParameter(i);
3275   Int_t oldFuncNpx = fConvolutedGausDeltaPrime->GetNpx();
3276   Double_t oldFuncRangeLow = 0, oldFuncRangeUp = 100;
3277   fConvolutedGausDeltaPrime->GetRange(oldFuncRangeLow, oldFuncRangeUp);
3278
3279   // Choose some sufficiently large range
3280   const Double_t rangeStart = 0.5;
3281   const Double_t rangeEnd = 2.0;
3282
3283   // To get the parameters for the transition, just choose arbitrarily input parameters for mu and sigma
3284   // (it makes sense to choose typical values). The ratio sigma_gauss / sigma_convolution is independent
3285   // of mu and as well the difference mu_gauss - mu_convolution.
3286   Double_t muInput = 1.0;
3287   Double_t sigmaInput = fgkSigmaReferenceForTransitionPars;
3288
3289
3290   // Step 1: Generate distribution with input parameters
3291   const Int_t nPar = 3;
3292   Double_t inputPar[nPar] = { muInput, sigmaInput, lambda };
3293
3294   TH1D* hInput = new TH1D("hInput", "Input distribution", 2000, rangeStart, rangeEnd);
3295
3296   fConvolutedGausDeltaPrime->SetParameters(inputPar);
3297   fConvolutedGausDeltaPrime->SetRange(rangeStart, rangeEnd);
3298   fConvolutedGausDeltaPrime->SetNpx(2000);
3299
3300   /*OLD
3301   // The resolution and mean of the AliPIDResponse are determined in finite intervals
3302   // of ncl (also second order effects due to finite dEdx and tanTheta).
3303   // Take this into account for the transition parameters via assuming an approximately flat
3304   // distritubtion in ncl in this interval.
3305   // NOTE: The ncl interval should be the same as the one used for the sigma map creation!
3306   const Int_t nclStart = 151;
3307   const Int_t nclEnd = 160;
3308   const Int_t nclSteps = (nclEnd - nclStart) + 1;
3309   for (Int_t ncl = nclStart; ncl <= nclEnd; ncl++) {
3310     // Resolution scales with 1/sqrt(ncl)
3311     fConvolutedGausDeltaPrime->SetParameter(1, inputPar[1] * sqrt(nclEnd) / sqrt(ncl));
3312     TH1* hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
3313
3314     for (Int_t i = 0; i < 50000000 / nclSteps; i++)
3315       hInput->Fill(hProbDensity->GetRandom());
3316   }
3317   */
3318
3319   TH1* hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
3320
3321   for (Int_t i = 0; i < 50000000; i++)
3322     hInput->Fill(hProbDensity->GetRandom());
3323
3324   // Step 2: Fit generated distribution with restricted gaussian
3325   Int_t maxBin = hInput->GetMaximumBin();
3326   Double_t max = hInput->GetBinContent(maxBin);
3327
3328   UChar_t usedBins = 1;
3329   if (maxBin > 1) {
3330     max += hInput->GetBinContent(maxBin - 1);
3331     usedBins++;
3332   }
3333   if (maxBin < hInput->GetNbinsX()) {
3334     max += hInput->GetBinContent(maxBin + 1);
3335     usedBins++;
3336   }
3337   max /= usedBins;
3338
3339   // NOTE: The range (<-> fraction of maximum) should be the same
3340   // as for the spline and eta maps creation
3341   const Double_t lowThreshold = hInput->GetXaxis()->GetBinLowEdge(hInput->FindFirstBinAbove(0.1 * max));
3342   const Double_t highThreshold = hInput->GetXaxis()->GetBinUpEdge(hInput->FindLastBinAbove(0.1 * max));
3343
3344   TFitResultPtr fitResGaussFirstStep = hInput->Fit("gaus", "QNRS", "", lowThreshold, highThreshold);
3345
3346   TFitResultPtr fitResGauss;
3347
3348   if ((Int_t)fitResGaussFirstStep == 0) {
3349     TF1 fGauss("fGauss", "[0]*TMath::Gaus(x, [1], [2], 1)", rangeStart, rangeEnd);
3350     fGauss.SetParameter(0, fitResGaussFirstStep->GetParams()[0]);
3351     fGauss.SetParError(0, fitResGaussFirstStep->GetErrors()[0]);
3352     fGauss.SetParameter(2, fitResGaussFirstStep->GetParams()[2]);
3353     fGauss.SetParError(2, fitResGaussFirstStep->GetErrors()[2]);
3354
3355     fGauss.FixParameter(1, fitResGaussFirstStep->GetParams()[1]);
3356     fitResGauss = hInput->Fit(&fGauss, "QNS", "s", rangeStart, rangeEnd);
3357   }
3358   else {
3359     fitResGauss = hInput->Fit("gaus", "QNRS", "same", rangeStart, rangeEnd);
3360   }
3361   //OLD TFitResultPtr fitResGauss = hInput->Fit("gaus", "QNRS", "", hInput->GetXaxis()->GetBinLowEdge(hInput->FindFirstBinAbove(0.1 * max)),
3362   //                                        hInput->GetXaxis()->GetBinUpEdge(hInput->FindLastBinAbove(0.1 * max)));
3363
3364
3365   // Step 3: Use parameters from gaussian fit to obtain parameters for the transition "restricted gauss" -> "convoluted gauss"
3366
3367   // 3.1 The ratio sigmaInput / sigma_gaussFit ONLY depends on lambda (which is fixed per period) -> Calculate this first
3368   // for an arbitrary (here: typical) sigma. The ratio is then ~the same for ALL sigma for given lambda!
3369   if ((Int_t)fitResGauss != 0) {
3370     AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Gauss Fit failed!\n");
3371     fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
3372     fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
3373     fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
3374
3375     delete hInput;
3376     delete [] oldFuncParams;
3377
3378     return kFALSE;
3379   }
3380
3381   if (fitResGauss->GetParams()[2] <= 0) {
3382     AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Sigma of gauss fit <= 0!\n");
3383     fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
3384     fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
3385     fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
3386
3387     delete hInput;
3388     delete [] oldFuncParams;
3389
3390     return kFALSE;
3391   }
3392
3393   // sigma correction factor
3394   fConvolutedGaussTransitionPars[1] = sigmaInput / fitResGauss->GetParams()[2];
3395
3396   // 3.2 Now that sigma und lambda are determined, one can calculate mu by shifting the maximum to the desired position,
3397   // i.e. the maximum of the input distribution should coincide with that of the re-calculated distribution,
3398   // which is achieved by getting the same mu for the same sigma.
3399   // NOTE: For fixed lambda, the shift is proportional to sigma and independent of mu!
3400   // So, one can calculate the shift for an arbitrary fixed (here: typical)
3401   // sigma and then simply use this shift for any other sigma by scaling it correspondingly!!!
3402
3403   // Mu shift correction:
3404   // Shift in mu (difference between mean of gauss and mean of convolution) is proportional to sigma!
3405   // Thus, choose a reference sigma (typical -> 0.05), such that for arbitrary sigma one can simple scale
3406   // this shift correction with sigma / referenceSigma.
3407   fConvolutedGaussTransitionPars[0] = (fitResGauss->GetParams()[1] - muInput);
3408
3409
3410   /*Changed on 03.07.2013
3411
3412   // Maximum of fConvolutedGausDeltaPrime should agree with maximum of input
3413   Double_t par[nPar] = { fitResGauss->GetParams()[1], // just as a guess of the maximum position
3414                          sigmaInput,
3415                          lambda };
3416
3417   fConvolutedGausDeltaPrime->SetParameters(par);
3418
3419   Double_t maxXInput = fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001, muInput - 3. * sigmaInput),
3420                                                               muInput + 10. * sigmaInput,
3421                                                               0.0001);
3422
3423   // Maximum shifts proportional to sigma and is linear in mu (~mean of gauss)
3424   // Maximum should be typically located within [gaussMean, gaussMean + 3 gaussSigma].
3425   // -> Larger search range for safety reasons (also: sigma and/or mean might be not 100% accurate).
3426   Double_t maxXconvoluted = fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001,
3427                                                                               fitResGauss->GetParams()[1] - 3. * fitResGauss->GetParams()[2]),
3428                                                                    fitResGauss->GetParams()[1] + 10. * fitResGauss->GetParams()[2],
3429                                                                    0.0001);
3430   if (maxXconvoluted <= 0) {
3431     AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Maximum of fConvolutedGausDeltaPrime <= 0!\n");
3432
3433     fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
3434     fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
3435     fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
3436
3437     delete hInput;
3438     delete [] oldFuncParams;
3439
3440     return kFALSE;
3441   }
3442
3443   // maxX perfectly shifts as par[0] (scaled by sigma) -> Can shift maxX to input value.
3444   // Mu shift correction:
3445   fConvolutedGaussTransitionPars[0] = maxXconvoluted - maxXInput;
3446   */
3447
3448
3449
3450   fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
3451   fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
3452   fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
3453
3454   delete hInput;
3455   delete [] oldFuncParams;
3456
3457   return kTRUE;
3458 }
3459
3460
3461 //_____________________________________________________________________________
3462 AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::SetParamsForConvolutedGaus(Double_t gausMean, Double_t gausSigma)
3463 {
3464   // Set parameters for convoluted gauss using parameters for a pure gaussian.
3465   // If SetConvolutedGaussLambdaParameter has not been called before to initialise the translation parameters,
3466   // some default parameters will be used and an error will show up.
3467
3468   if(fDebug > 1)
3469     printf("File: %s, Line: %d: SetParamsForConvolutedGaus: mean %e, sigma %e\n", (char*)__FILE__, __LINE__, gausMean, gausSigma);
3470
3471   if (fConvolutedGaussTransitionPars[1] < -998) {
3472     AliError("Transition parameters not initialised! Default parameters will be used. Please call SetConvolutedGaussLambdaParameter(...) before any calculations!");
3473     SetConvolutedGaussLambdaParameter(2.0);
3474     AliError(Form("Parameters set to:\n[0]: %f\n[1]: %f\n[2]: %f\n", fConvolutedGaussTransitionPars[0],
3475                   fConvolutedGaussTransitionPars[1], fConvolutedGaussTransitionPars[2]));
3476   }
3477
3478   Double_t par[fkConvolutedGausNPar];
3479   par[2] = fConvolutedGaussTransitionPars[2];
3480   par[1] = fConvolutedGaussTransitionPars[1] * gausSigma;
3481   // maxX perfectly shifts as par[0] (scaled by sigma) -> Can shift maxX so that it sits at the right place.
3482   par[0] = gausMean - fConvolutedGaussTransitionPars[0] * par[1] / fgkSigmaReferenceForTransitionPars;
3483
3484   ErrorCode errCode = kNoErrors;
3485   fConvolutedGausDeltaPrime->SetParameters(par);
3486
3487   if(fDebug > 3)
3488     printf("File: %s, Line: %d: SetParamsForConvolutedGaus -> Parameters set to: %e, %e, %e (transition pars: %e, %e, %e, %e)\n",
3489            (char*)__FILE__, __LINE__, par[0], par[1], par[2], fConvolutedGaussTransitionPars[0], fConvolutedGaussTransitionPars[1],
3490            fConvolutedGaussTransitionPars[2], fgkSigmaReferenceForTransitionPars);
3491
3492   fConvolutedGausDeltaPrime->SetNpx(20); // Small value speeds up following algorithm (valid, since extrema far apart)
3493
3494   // Accuracy of 10^-5 is enough to get 0.1% precise peak for MIPS w.r.t. to dEdx = 2000 of protons
3495   // (should boost up the algorithm, because 10^-10 is the default value!)
3496   Double_t maxX= fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001, gausMean - 2. * gausSigma),
3497                                                         gausMean + 6. * gausSigma, 1.0E-5);
3498
3499   const Double_t maximum = fConvolutedGausDeltaPrime->Eval(maxX);
3500   const Double_t maximumFraction = maximum * fAccuracyNonGaussianTail;
3501
3502   // Estimate lower boundary for subsequent search:
3503   Double_t lowBoundSearchBoundLow = TMath::Max(1e-4, maxX - 5. * gausSigma);
3504   Double_t lowBoundSearchBoundUp = maxX;
3505
3506   Bool_t lowerBoundaryFixedAtZero = kFALSE;
3507
3508   while (fConvolutedGausDeltaPrime->Eval(lowBoundSearchBoundLow) >= maximumFraction) {
3509     if (lowBoundSearchBoundLow <= 0) {
3510       // This should only happen to low dEdx particles with very few clusters and therefore large sigma, such that the gauss goes below zero deltaPrime
3511       if (maximum <= 0) { // Something is weired
3512         printf("Error generating signal: maximum is <= 0!\n");
3513         return kError;
3514       }
3515       else {
3516         const Double_t valueAtZero = fConvolutedGausDeltaPrime->Eval(0);
3517         if (valueAtZero / maximum > 0.05) {
3518           // Too large fraction below zero deltaPrime. Signal generation cannot be reliable in this case
3519           printf("Error generating signal: Too large fraction below zero deltaPrime: convGauss(0) / convGauss(max) =  %e / %e = %e!\n",
3520                  valueAtZero, maximum, valueAtZero / maximum);
3521           return kError;
3522         }
3523       }
3524
3525       /*
3526       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",
3527              fConvolutedGausDeltaPrime->Eval(maxX), fAccuracyNonGaussianTail, maximumFraction, maxX, par[0], par[1], par[2], gausMean, gausSigma);
3528       */
3529
3530       lowerBoundaryFixedAtZero = kTRUE;
3531
3532       if (errCode != kError)
3533         errCode = kWarning;
3534
3535       break;
3536     }
3537
3538     lowBoundSearchBoundUp -= gausSigma;
3539     lowBoundSearchBoundLow -= gausSigma;
3540
3541     if (lowBoundSearchBoundLow < 0) {
3542       lowBoundSearchBoundLow = 0;
3543       lowBoundSearchBoundUp += gausSigma;
3544     }
3545   }
3546
3547   // Determine lower boundary inside estimated range. For small values of the maximum: Need more precision, since finer binning!
3548   Double_t rangeStart = lowerBoundaryFixedAtZero ? 0 :
3549                         fConvolutedGausDeltaPrime->GetX(maximumFraction, lowBoundSearchBoundLow, lowBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
3550
3551   // .. and the same for the upper boundary
3552   Double_t rangeEnd = 0;
3553   // If distribution starts beyond upper boundary, everything ends up in the overflow bin. So, just reduce range and Npx to minimum
3554   if (rangeStart > fkDeltaPrimeUpLimit) {
3555     rangeEnd = rangeStart + 0.00001;
3556     fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
3557     fConvolutedGausDeltaPrime->SetNpx(4);
3558   }
3559   else {
3560     // Estimate upper boundary for subsequent search:
3561     Double_t upBoundSearchBoundUp = maxX + 5 * gausSigma;
3562     Double_t upBoundSearchBoundLow = maxX;
3563     while (fConvolutedGausDeltaPrime->Eval(upBoundSearchBoundUp) >= maximumFraction) {
3564       upBoundSearchBoundUp += gausSigma;
3565       upBoundSearchBoundLow += gausSigma;
3566     }
3567
3568     //  For small values of the maximum: Need more precision, since finer binning!
3569     rangeEnd = fConvolutedGausDeltaPrime->GetX(maximumFraction, upBoundSearchBoundLow, upBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
3570
3571     fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
3572     fConvolutedGausDeltaPrime->SetNpx(fDeltaPrimeAxis->FindFixBin(rangeEnd) - fDeltaPrimeAxis->FindFixBin(rangeStart) + 1);
3573
3574     fConvolutedGausDeltaPrime->SetNpx(fDeltaPrimeAxis->FindFixBin(rangeEnd) - fDeltaPrimeAxis->FindFixBin(rangeStart) + 1);
3575     //fConvolutedGausDeltaPrime->SetNpx(fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindFixBin(rangeEnd)
3576     //                                  - fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindFixBin(rangeStart) + 1);
3577     //fConvolutedGausDeltaPrime->SetNpx((rangeEnd - rangeStart) / fDeltaPrimeBinWidth + 1);
3578   }
3579
3580   if(fDebug > 3)
3581     printf("File: %s, Line: %d: SetParamsForConvolutedGaus -> range %f - %f, error code %d\n", (char*)__FILE__, __LINE__,
3582            rangeStart, rangeEnd, errCode);
3583
3584   return errCode;
3585 }
3586
3587
3588 //________________________________________________________________________
3589 void AliAnalysisTaskPID::SetUpGenHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const
3590 {
3591   // Sets bin limits for axes which are not standard binned and the axes titles.
3592
3593   hist->SetBinEdges(kGenPt, binsPt);
3594   hist->SetBinEdges(kGenDeltaPrimeSpecies, binsDeltaPrime);
3595   hist->SetBinEdges(kGenCentrality, binsCent);
3596
3597   if (fStoreAdditionalJetInformation)
3598     hist->SetBinEdges(kGenJetPt, binsJetPt);
3599
3600   // Set axes titles
3601   hist->GetAxis(kGenMCID)->SetTitle("MC PID");
3602   hist->GetAxis(kGenMCID)->SetBinLabel(1, "e");
3603   hist->GetAxis(kGenMCID)->SetBinLabel(2, "K");
3604   hist->GetAxis(kGenMCID)->SetBinLabel(3, "#mu");
3605   hist->GetAxis(kGenMCID)->SetBinLabel(4, "#pi");
3606   hist->GetAxis(kGenMCID)->SetBinLabel(5, "p");
3607
3608   hist->GetAxis(kGenSelectSpecies)->SetTitle("Select Species");
3609   hist->GetAxis(kGenSelectSpecies)->SetBinLabel(1, "e");
3610   hist->GetAxis(kGenSelectSpecies)->SetBinLabel(2, "K");
3611   hist->GetAxis(kGenSelectSpecies)->SetBinLabel(3, "#pi");
3612   hist->GetAxis(kGenSelectSpecies)->SetBinLabel(4, "p");
3613
3614   hist->GetAxis(kGenPt)->SetTitle("p_{T} (GeV/c)");
3615
3616   hist->GetAxis(kGenDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
3617
3618   hist->GetAxis(kGenCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
3619
3620   if (fStoreAdditionalJetInformation) {
3621     hist->GetAxis(kGenJetPt)->SetTitle("p_{T}^{jet} (GeV/c)");
3622
3623     hist->GetAxis(kGenZ)->SetTitle("z = p_{T}^{track} / p_{T}^{jet}");
3624
3625     hist->GetAxis(kGenXi)->SetTitle("#xi = ln(p_{T}^{jet} / p_{T}^{track})");
3626   }
3627
3628   hist->GetAxis(GetIndexOfChargeAxisGen())->SetTitle("Charge (e_{0})");
3629
3630   hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetTitle("TOF PID Info");
3631   // Offset is (TMath::Abs(kNoTOFinfo) + 1), such that bin 1 (= first label) corresponds to kNoTOFinfo (< 0)
3632   hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kNoTOFinfo + (TMath::Abs(kNoTOFinfo) + 1), "No TOF Info");
3633   hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kNoTOFpid + (TMath::Abs(kNoTOFinfo) + 1), "No TOF PID");
3634   hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kTOFpion + (TMath::Abs(kNoTOFinfo) + 1), "#pi");
3635   hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kTOFkaon + (TMath::Abs(kNoTOFinfo) + 1), "K");
3636   hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kTOFproton + (TMath::Abs(kNoTOFinfo) + 1), "p");
3637 }
3638
3639
3640 //________________________________________________________________________
3641 void AliAnalysisTaskPID::SetUpGenYieldHist(THnSparse* hist, Double_t* binsPt, Double_t* binsCent, Double_t* binsJetPt) const
3642 {
3643   // Sets bin limits for axes which are not standard binned and the axes titles.
3644
3645   hist->SetBinEdges(kGenYieldPt, binsPt);
3646   hist->SetBinEdges(kGenYieldCentrality, binsCent);
3647   if (fStoreAdditionalJetInformation)
3648     hist->SetBinEdges(kGenYieldJetPt, binsJetPt);
3649
3650   for (Int_t i = 0; i < 5; i++)
3651     hist->GetAxis(kGenYieldMCID)->SetBinLabel(i + 1, AliPID::ParticleLatexName(i));
3652
3653   // Set axes titles
3654   hist->GetAxis(kGenYieldMCID)->SetTitle("MC PID");
3655   hist->GetAxis(kGenYieldPt)->SetTitle("p_{T}^{gen} (GeV/c)");
3656   hist->GetAxis(kGenYieldCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
3657
3658   if (fStoreAdditionalJetInformation) {
3659     hist->GetAxis(kGenYieldJetPt)->SetTitle("p_{T}^{jet, gen} (GeV/c)");
3660
3661     hist->GetAxis(kGenYieldZ)->SetTitle("z = p_{T}^{track} / p_{T}^{jet}");
3662
3663     hist->GetAxis(kGenYieldXi)->SetTitle("#xi = ln(p_{T}^{jet} / p_{T}^{track})");
3664   }
3665
3666   hist->GetAxis(GetIndexOfChargeAxisGenYield())->SetTitle("Charge (e_{0})");
3667 }
3668
3669
3670 //________________________________________________________________________
3671 void AliAnalysisTaskPID::SetUpHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const
3672 {
3673   // Sets bin limits for axes which are not standard binned and the axes titles.
3674
3675   hist->SetBinEdges(kDataPt, binsPt);
3676   hist->SetBinEdges(kDataDeltaPrimeSpecies, binsDeltaPrime);
3677   hist->SetBinEdges(kDataCentrality, binsCent);
3678
3679   if (fStoreAdditionalJetInformation)
3680     hist->SetBinEdges(kDataJetPt, binsJetPt);
3681
3682   // Set axes titles
3683   hist->GetAxis(kDataMCID)->SetTitle("MC PID");
3684   hist->GetAxis(kDataMCID)->SetBinLabel(1, "e");
3685   hist->GetAxis(kDataMCID)->SetBinLabel(2, "K");
3686   hist->GetAxis(kDataMCID)->SetBinLabel(3, "#mu");
3687   hist->GetAxis(kDataMCID)->SetBinLabel(4, "#pi");
3688   hist->GetAxis(kDataMCID)->SetBinLabel(5, "p");
3689
3690   hist->GetAxis(kDataSelectSpecies)->SetTitle("Select Species");
3691   hist->GetAxis(kDataSelectSpecies)->SetBinLabel(1, "e");
3692   hist->GetAxis(kDataSelectSpecies)->SetBinLabel(2, "K");
3693   hist->GetAxis(kDataSelectSpecies)->SetBinLabel(3, "#pi");
3694   hist->GetAxis(kDataSelectSpecies)->SetBinLabel(4, "p");
3695
3696   hist->GetAxis(kDataPt)->SetTitle("p_{T} (GeV/c)");
3697
3698   hist->GetAxis(kDataDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
3699
3700   hist->GetAxis(kDataCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
3701
3702   if (fStoreAdditionalJetInformation) {
3703     hist->GetAxis(kDataJetPt)->SetTitle("p_{T}^{jet} (GeV/c)");
3704
3705     hist->GetAxis(kDataZ)->SetTitle("z = p_{T}^{track} / p_{T}^{jet}");
3706
3707     hist->GetAxis(kDataXi)->SetTitle("#xi = ln(p_{T}^{jet} / p_{T}^{track})");
3708   }
3709
3710   hist->GetAxis(GetIndexOfChargeAxisData())->SetTitle("Charge (e_{0})");
3711
3712   hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetTitle("TOF PID Info");
3713   // Offset is (TMath::Abs(kNoTOFinfo) + 1), such that bin 1 (= first label) corresponds to kNoTOFinfo (< 0)
3714   hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kNoTOFinfo + (TMath::Abs(kNoTOFinfo) + 1), "No TOF Info");
3715   hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kNoTOFpid + (TMath::Abs(kNoTOFinfo) + 1), "No TOF PID");
3716   hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kTOFpion + (TMath::Abs(kNoTOFinfo) + 1), "#pi");
3717   hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kTOFkaon + (TMath::Abs(kNoTOFinfo) + 1), "K");
3718   hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kTOFproton + (TMath::Abs(kNoTOFinfo) + 1), "p");
3719 }
3720
3721
3722 //________________________________________________________________________
3723 void AliAnalysisTaskPID::SetUpPtResHist(THnSparse* hist, Double_t* binsPt, Double_t* binsJetPt, Double_t* binsCent) const
3724 {
3725   // Sets bin limits for axes which are not standard binned and the axes titles.
3726
3727   hist->SetBinEdges(kPtResJetPt, binsJetPt);
3728   hist->SetBinEdges(kPtResGenPt, binsPt);
3729   hist->SetBinEdges(kPtResRecPt, binsPt);
3730   hist->SetBinEdges(kPtResCentrality, binsCent);
3731
3732   // Set axes titles
3733   hist->GetAxis(kPtResJetPt)->SetTitle("p_{T}^{jet, rec} (GeV/c)");
3734   hist->GetAxis(kPtResGenPt)->SetTitle("p_{T}^{gen} (GeV/c)");
3735   hist->GetAxis(kPtResRecPt)->SetTitle("p_{T}^{rec} (GeV/c)");
3736
3737   hist->GetAxis(kPtResCharge)->SetTitle("Charge (e_{0})");
3738   hist->GetAxis(kPtResCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
3739 }
3740
3741
3742 //________________________________________________________________________
3743 void AliAnalysisTaskPID::SetUpSharedClsHist(THnSparse* hist, Double_t* binsPt, Double_t* binsJetPt) const
3744 {
3745   // Sets bin limits for axes which are not standard binned and the axes titles.
3746
3747   hist->SetBinEdges(kQASharedClsJetPt, binsJetPt);
3748   hist->SetBinEdges(kQASharedClsPt, binsPt);
3749
3750   // Set axes titles
3751   hist->GetAxis(kQASharedClsJetPt)->SetTitle("#it{p}_{T}^{jet} (GeV/#it{c})");
3752   hist->GetAxis(kQASharedClsPt)->SetTitle("#it{p}_{T} (GeV/#it{c})");
3753   hist->GetAxis(kQASharedClsNumSharedCls)->SetTitle("#it{N}_{shared}^{cls}");
3754   hist->GetAxis(kQASharedClsPadRow)->SetTitle("Pad row");
3755
3756 }
3757
3758
3759 //________________________________________________________________________
3760 void AliAnalysisTaskPID::SetUpDeDxCheckHist(THnSparse* hist, const Double_t* binsPt, const Double_t* binsJetPt, const Double_t* binsEtaAbs) const
3761 {
3762   // Sets bin limits for axes which are not standard binned and the axes titles.
3763   hist->SetBinEdges(kDeDxCheckP, binsPt);
3764   hist->SetBinEdges(kDeDxCheckJetPt, binsJetPt);
3765   hist->SetBinEdges(kDeDxCheckEtaAbs, binsEtaAbs);
3766
3767
3768   // Set axes titles
3769   hist->GetAxis(kDeDxCheckPID)->SetTitle("PID");
3770   hist->GetAxis(kDeDxCheckPID)->SetBinLabel(1, "e");
3771   hist->GetAxis(kDeDxCheckPID)->SetBinLabel(2, "K");
3772   hist->GetAxis(kDeDxCheckPID)->SetBinLabel(3, "#pi");
3773   hist->GetAxis(kDeDxCheckPID)->SetBinLabel(4, "p");
3774
3775
3776   hist->GetAxis(kDeDxCheckJetPt)->SetTitle("p_{T}^{jet} (GeV/c)");
3777   hist->GetAxis(kDeDxCheckEtaAbs)->SetTitle("|#eta|");
3778   hist->GetAxis(kDeDxCheckP)->SetTitle("p_{TPC} (GeV/c)");
3779   hist->GetAxis(kDeDxCheckDeDx)->SetTitle("TPC dE/dx (arb. unit)");
3780 }
3781
3782
3783 //________________________________________________________________________
3784 void AliAnalysisTaskPID::SetUpBinZeroStudyHist(THnSparse* hist, const Double_t* binsCent, const Double_t* binsPt) const
3785 {
3786   // Sets bin limits for axes which are not standard binned and the axes titles.
3787   hist->SetBinEdges(kBinZeroStudyCentrality, binsCent);
3788   hist->SetBinEdges(kBinZeroStudyGenPt, binsPt);
3789
3790   // Set axes titles
3791   hist->GetAxis(kBinZeroStudyCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
3792   hist->GetAxis(kBinZeroStudyGenEta)->SetTitle("#it{#eta}^{gen}");
3793   hist->GetAxis(kBinZeroStudyGenPt)->SetTitle("#it{p}_{T}^{gen} (GeV/#it{c})");
3794 }