PWGJE/UserTasks/AliAnalysisTaskPID.cxx

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