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