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