]> 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
Update of the ITS version number
[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
e4351829 767 const Int_t nPtBinsRes = 100;
768 Double_t pTbinsRes[nPtBinsRes + 1];
769
770 const Double_t fromLowPtRes = 0.15;
771 const Double_t toHighPtRes = 50.;
772 const Double_t factorPtRes = TMath::Power(toHighPtRes/fromLowPtRes, 1./nPtBinsRes);
773 // Log binning for whole pT range
774 pTbinsRes[0] = fromLowPtRes;
775 for (Int_t i = 0 + 1; i <= nPtBinsRes; i++) {
776 pTbinsRes[i] = factorPtRes * pTbinsRes[i - 1];
777 }
778
9e95a906 779 const Int_t nBinsPtResolution = kPtResNumAxes;
e4351829 780 Int_t ptResolutionBins[kPtResNumAxes] = { nJetPtBins, nPtBinsRes, nPtBinsRes,
781 nChargeBins, nCentBins };
782 Double_t ptResolutionXmin[kPtResNumAxes] = { binsJetPt[0], pTbinsRes[0], pTbinsRes[0],
783 binsCharge[0], binsCent[0] };
784 Double_t ptResolutionXmax[kPtResNumAxes] = { binsJetPt[nJetPtBins], pTbinsRes[nPtBinsRes], pTbinsRes[nPtBinsRes],
785 binsCharge[nChargeBins], binsCent[nCentBins] };
9e95a906 786
787 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
788 fPtResolution[i] = new THnSparseD(Form("fPtResolution_%s", AliPID::ParticleShortName(i)),
789 Form("Pt resolution for primaries, %s", AliPID::ParticleLatexName(i)),
790 nBinsPtResolution, ptResolutionBins, ptResolutionXmin, ptResolutionXmax);
e4351829 791 SetUpPtResHist(fPtResolution[i], pTbinsRes, binsJetPt, binsCent);
9e95a906 792 fPtResolutionContainer->Add(fPtResolution[i]);
793 }
794 }
795
796 if(fDebug > 2)
797 printf("File: %s, Line: %d: UserCreateOutputObjects -> Posting output data\n", (char*)__FILE__, __LINE__);
798
e131b05f 799 PostOutputData();
9e95a906 800
801 if(fDebug > 2)
802 printf("File: %s, Line: %d: UserCreateOutputObjects -> Done\n", (char*)__FILE__, __LINE__);
e131b05f 803}
804
805
806//________________________________________________________________________
807void AliAnalysisTaskPID::UserExec(Option_t *)
808{
809 // Main loop
810 // Called for each event
9d7ad2e4 811
9e95a906 812 if(fDebug > 1)
813 printf("File: %s, Line: %d: UserExec\n", (char*)__FILE__, __LINE__);
814
e131b05f 815 // No processing of event, if input is fed in directly from another task
816 if (fInputFromOtherTask)
817 return;
9e95a906 818
819 if(fDebug > 1)
820 printf("File: %s, Line: %d: UserExec -> Processing started\n", (char*)__FILE__, __LINE__);
e131b05f 821
822 fEvent = dynamic_cast<AliVEvent*>(InputEvent());
823 if (!fEvent) {
824 Printf("ERROR: fEvent not available");
825 return;
826 }
827
828 fMC = dynamic_cast<AliMCEvent*>(MCEvent());
829
830 if (!fPIDResponse || !fPIDcombined)
831 return;
832
833 if (!GetVertexIsOk(fEvent))
834 return;
835
836 fESD = dynamic_cast<AliESDEvent*>(fEvent);
837 const AliVVertex* primaryVertex = fESD ? fESD->GetPrimaryVertexTracks() : fEvent->GetPrimaryVertex();
838 if (!primaryVertex)
839 return;
840
841 if(primaryVertex->GetNContributors() <= 0)
842 return;
843
844 Double_t magField = fEvent->GetMagneticField();
845
846 //OLD with DeltaSpecies const Bool_t usePureGausForDelta = kTRUE;
847
848
849 Double_t centralityPercentile = -1;
850 if (fStoreCentralityPercentile)
851 centralityPercentile = fEvent->GetCentrality()->GetCentralityPercentile(fCentralityEstimator.Data());
852
853 if (fMC) {
9e95a906 854 if (fDoPID || fDoEfficiency) {
855 for (Int_t iPart = 0; iPart < fMC->GetNumberOfTracks(); iPart++) {
856 AliMCParticle *mcPart = dynamic_cast<AliMCParticle*>(fMC->GetTrack(iPart));
857
858 if (!mcPart)
859 continue;
860
861 // Define clean MC sample with corresponding particle level track cuts:
862 // - MC-track must be in desired eta range
863 // - MC-track must be physical primary
864 // - Species must be one of those in question (everything else goes to the overflow bin of mcID)
865
866 // Geometrie should be the same as on the reconstructed level -> By definition analysis within this eta interval
867 if (TMath::Abs(mcPart->Eta()) < fEtaAbsCutLow || TMath::Abs(mcPart->Eta()) > fEtaAbsCutUp) continue;
868
869 Int_t mcID = PDGtoMCID(mcPart->PdgCode());
870
871 // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
872 Double_t chargeMC = mcPart->Charge() / 3.;
873
874 if (TMath::Abs(chargeMC) < 0.01)
875 continue; // Reject neutral particles (only relevant, if mcID is not used)
876
877 if (!fMC->IsPhysicalPrimary(iPart))
878 continue;
879
880 if (fDoPID) {
881 Double_t valuesGenYield[kGenYieldNumAxes] = { mcID, mcPart->Pt(), centralityPercentile, -1, -1, -1, -1 };
882 valuesGenYield[GetIndexOfChargeAxisGenYield()] = chargeMC;
883
884 fhMCgeneratedYieldsPrimaries->Fill(valuesGenYield);
885 }
886
887
888 if (fDoEfficiency) {
889 Double_t valueEff[kEffNumAxes] = { mcID, mcPart->Pt(), mcPart->Eta(), chargeMC, centralityPercentile,
890 -1, -1, -1 };
891
892 fContainerEff->Fill(valueEff, kStepGenWithGenCuts);
893 }
894 }
e131b05f 895 }
896 }
897
898 // Track loop to fill a Train spectrum
899 for (Int_t iTracks = 0; iTracks < fEvent->GetNumberOfTracks(); iTracks++) {
900 AliVTrack* track = dynamic_cast<AliVTrack*>(fEvent->GetTrack(iTracks));
901 if (!track) {
902 Printf("ERROR: Could not retrieve track %d", iTracks);
903 continue;
904 }
905
906
907 // Apply detector level track cuts
a6852ea8 908 Double_t dEdxTPC = fPIDResponse->IsTunedOnData() ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
909 if (dEdxTPC <= 0)
910 continue;
e131b05f 911
912 if(fTrackFilter && !fTrackFilter->IsSelected(track))
913 continue;
914
493982d9 915 if (GetUseTPCCutMIGeo()) {
a6852ea8 916 if (!TPCCutMIGeo(track, fEvent))
917 continue;
918 }
493982d9
ML		 919 else if (GetUseTPCnclCut()) {
920 if (!TPCnclCut(track))
921 continue;
922 }
e131b05f 923
924 if(fUsePhiCut) {
925 if (!PhiPrimeCut(track, magField))
926 continue; // reject track
927 }
928
e131b05f 929 Int_t pdg = 0; // = 0 indicates data for the moment
930 AliMCParticle* mcTrack = 0x0;
931 Int_t mcID = AliPID::kUnknown;
932 Int_t label = 0;
933
934 if (fMC) {
935 label = track->GetLabel();
936
937 //if (label < 0)
938 // continue;
939
940 mcTrack = dynamic_cast<AliMCParticle*>(fMC->GetTrack(TMath::Abs(label)));
941 if (!mcTrack) {
942 Printf("ERROR: Could not retrieve mcTrack with label %d for track %d", label, iTracks);
943 continue;
944 }
945
946 pdg = mcTrack->PdgCode();
947 mcID = PDGtoMCID(mcTrack->PdgCode());
948
9e95a906 949 if (fDoEfficiency) {
950 // For efficiency: Reconstructed track has survived all cuts on the detector level (excluding acceptance)
951 // and has an associated MC track which is a physical primary and was generated inside the acceptance
952 if (fMC->IsPhysicalPrimary(TMath::Abs(label)) &&
953 (TMath::Abs(mcTrack->Eta()) >= fEtaAbsCutLow && TMath::Abs(mcTrack->Eta()) <= fEtaAbsCutUp)) {
e131b05f 954
9e95a906 955 // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
956 Double_t value[kEffNumAxes] = { mcID, mcTrack->Pt(), mcTrack->Eta(), mcTrack->Charge() / 3., centralityPercentile,
957 -1, -1, -1 };
958 fContainerEff->Fill(value, kStepRecWithGenCuts);
959
960 Double_t valueMeas[kEffNumAxes] = { mcID, track->Pt(), track->Eta(), track->Charge(), centralityPercentile,
961 -1, -1, -1 };
962 fContainerEff->Fill(valueMeas, kStepRecWithGenCutsMeasuredObs);
963 }
e131b05f 964 }
965 }
966
967 // Only process tracks inside the desired eta window
968 if (TMath::Abs(track->Eta()) < fEtaAbsCutLow || TMath::Abs(track->Eta()) > fEtaAbsCutUp) continue;
969
9e95a906 970 if (fDoPID)
971 ProcessTrack(track, pdg, centralityPercentile, -1); // No jet information in this case -> Set jet pT to -1
e131b05f 972
9e95a906 973 if (fDoPtResolution) {
974 if (mcTrack && fMC->IsPhysicalPrimary(TMath::Abs(label))) {
e4351829 975 // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
976 Double_t valuePtRes[kPtResNumAxes] = { -1, mcTrack->Pt(), track->Pt(), mcTrack->Charge() / 3., centralityPercentile };
9e95a906 977 fPtResolution[mcID]->Fill(valuePtRes);
978 }
979 }
980
981 if (fDoEfficiency) {
982 if (mcTrack) {
983 Double_t valueRecAllCuts[kEffNumAxes] = { mcID, track->Pt(), track->Eta(), track->Charge(), centralityPercentile,
984 -1, -1, -1 };
985 fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObs);
986
987 Double_t weight = IsSecondaryWithStrangeMotherMC(fMC, TMath::Abs(label)) ?
988 GetMCStrangenessFactorCMS(fMC, mcTrack) : 1.0;
989 fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObsStrangenessScaled, weight);
990
991 // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
992 Double_t valueGenAllCuts[kEffNumAxes] = { mcID, mcTrack->Pt(), mcTrack->Eta(), mcTrack->Charge() / 3.,
993 centralityPercentile, -1, -1, -1 };
994 if (fMC->IsPhysicalPrimary(TMath::Abs(label))) {
995 fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObsPrimaries);
996 fContainerEff->Fill(valueGenAllCuts, kStepRecWithRecCutsPrimaries);
997 }
998 }
e131b05f 999 }
1000 } //track loop
1001
1002 IncrementEventsProcessed(centralityPercentile);
1003
9e95a906 1004 if(fDebug > 2)
1005 printf("File: %s, Line: %d: UserExec -> Processing done\n", (char*)__FILE__, __LINE__);
1006
e131b05f 1007 PostOutputData();
9e95a906 1008
1009 if(fDebug > 2)
1010 printf("File: %s, Line: %d: UserExec -> Done\n", (char*)__FILE__, __LINE__);
e131b05f 1011}
1012
1013//________________________________________________________________________
1014void AliAnalysisTaskPID::Terminate(const Option_t *)
1015{
1016 // Draw result to the screen
1017 // Called once at the end of the query
1018}
1019
1020
1021//_____________________________________________________________________________
1022void AliAnalysisTaskPID::CheckDoAnyStematicStudiesOnTheExpectedSignal()
1023{
1024 // Check whether at least one scale factor indicates the ussage of systematic studies
1025 // and set internal flag accordingly.
1026
1027 fDoAnySystematicStudiesOnTheExpectedSignal = kFALSE;
1028
1029 if (TMath::Abs(fSystematicScalingSplines - 1.0) > fgkEpsilon) {
1030 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1031 return;
1032 }
1033
1034 if ((TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon) ||
1035 (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon)) {
1036 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1037 return;
1038 }
1039
1040 if (TMath::Abs(fSystematicScalingEtaSigmaPara - 1.0) > fgkEpsilon) {
1041 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1042 return;
1043 }
1044
1045 if (TMath::Abs(fSystematicScalingMultCorrection - 1.0) > fgkEpsilon) {
1046 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1047 return;
1048 }
1049}
1050
1051
1052//_____________________________________________________________________________
1053Int_t AliAnalysisTaskPID::PDGtoMCID(Int_t pdg)
1054{
1055 // Returns the corresponding AliPID index to the given pdg code.
1056 // Returns AliPID::kUnkown if pdg belongs to a not considered species.
1057
1058 Int_t absPDGcode = TMath::Abs(pdg);
1059 if (absPDGcode == 211) {//Pion
1060 return AliPID::kPion;
1061 }
1062 else if (absPDGcode == 321) {//Kaon
1063 return AliPID::kKaon;
1064 }
1065 else if (absPDGcode == 2212) {//Proton
1066 return AliPID::kProton;
1067 }
1068 else if (absPDGcode == 11) {//Electron
1069 return AliPID::kElectron;
1070 }
1071 else if (absPDGcode == 13) {//Muon
1072 return AliPID::kMuon;
1073 }
1074
1075 return AliPID::kUnknown;
1076}
1077
1078
1079//_____________________________________________________________________________
9d7ad2e4 1080void AliAnalysisTaskPID::GetJetTrackObservables(Double_t trackPt, Double_t jetPt, Double_t& z, Double_t& xi)
e131b05f 1081{
1082 // Uses trackPt and jetPt to obtain z and xi.
1083
1084 z = (jetPt > 0 && trackPt >= 0) ? (trackPt / jetPt) : -1;
1085 xi = (z > 0) ? TMath::Log(1. / z) : -1;
1086
1087 if(trackPt > (1. - 1e-06) * jetPt && trackPt < (1. + 1e-06) * jetPt) { // case z=1 : move entry to last histo bin <1
1088 z = 1. - 1e-06;
1089 xi = 1e-06;
1090 }
1091}
1092
1093
1094//_____________________________________________________________________________
1095void AliAnalysisTaskPID::CleanupParticleFractionHistos()
1096{
1097 // Delete histos with particle fractions
1098
1099 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1100 delete fFractionHists[i];
1101 fFractionHists[i] = 0x0;
1102
1103 delete fFractionSysErrorHists[i];
1104 fFractionSysErrorHists[i] = 0x0;
1105 }
1106}
1107
1108
1109//_____________________________________________________________________________
1110Double_t AliAnalysisTaskPID::ConvolutedGaus(const Double_t* xx, const Double_t* par) const
1111{
1112 // Convolutes gauss with an exponential tail which describes dEdx-response better than pure gaussian
1113
1114 const Double_t mean = par[0];
1115 const Double_t sigma = par[1];
1116 const Double_t lambda = par[2];
1117
9e95a906 1118 if(fDebug > 5)
1119 printf("File: %s, Line: %d: ConvolutedGaus: mean %e, sigma %e, lambda %e\n", (char*)__FILE__, __LINE__, mean, sigma, lambda);
1120
e131b05f 1121 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);
1122}
1123
1124
1125//_____________________________________________________________________________
9d7ad2e4 1126inline Double_t AliAnalysisTaskPID::FastGaus(Double_t x, Double_t mean, Double_t sigma) const
e131b05f 1127{
1128 // Calculate an unnormalised gaussian function with mean and sigma.
1129
1130 if (sigma < fgkEpsilon)
1131 return 1.e30;
1132
1133 const Double_t arg = (x - mean) / sigma;
1134 return exp(-0.5 * arg * arg);
1135}
1136
1137
1138//_____________________________________________________________________________
9d7ad2e4 1139inline Double_t AliAnalysisTaskPID::FastNormalisedGaus(Double_t x, Double_t mean, Double_t sigma) const
e131b05f 1140{
1141 // Calculate a normalised (divided by sqrt(2*Pi)*sigma) gaussian function with mean and sigma.
1142
1143 if (sigma < fgkEpsilon)
1144 return 1.e30;
1145
1146 const Double_t arg = (x - mean) / sigma;
1147 const Double_t res = exp(-0.5 * arg * arg);
1148 return res / (2.50662827463100024 * sigma); //sqrt(2*Pi)=2.50662827463100024
1149}
1150
1151
1152//_____________________________________________________________________________
1153Int_t AliAnalysisTaskPID::FindBinWithinRange(TAxis* axis, Double_t value) const
1154{
1155 // Find the corresponding bin of the axis. Values outside the range (also under and overflow) will be set to the first/last
1156 // available bin
1157
1158 Int_t bin = axis->FindFixBin(value);
1159
1160 if (bin <= 0)
1161 bin = 1;
1162 if (bin > axis->GetNbins())
1163 bin = axis->GetNbins();
1164
1165 return bin;
1166}
1167
1168
1169//_____________________________________________________________________________
9d7ad2e4 1170Int_t AliAnalysisTaskPID::FindFirstBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yBin,
1171 Int_t zBin) const
e131b05f 1172{
1173 // Kind of projects a TH3 to 1 bin combination in y and z
1174 // and looks for the first x bin above a threshold for this projection.
1175 // If no such bin is found, -1 is returned.
1176
1177 if (!hist)
1178 return -1;
1179
1180 Int_t nBinsX = hist->GetNbinsX();
1181 for (Int_t xBin = 1; xBin <= nBinsX; xBin++) {
1182 if (hist->GetBinContent(xBin, yBin, zBin) > threshold)
1183 return xBin;
1184 }
1185
1186 return -1;
1187}
1188
1189
1190//_____________________________________________________________________________
9d7ad2e4 1191Int_t AliAnalysisTaskPID::FindLastBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yBin,
1192 Int_t zBin) const
e131b05f 1193{
1194 // Kind of projects a TH3 to 1 bin combination in y and z
1195 // and looks for the last x bin above a threshold for this projection.
1196 // If no such bin is found, -1 is returned.
1197
1198 if (!hist)
1199 return -1;
1200
1201 Int_t nBinsX = hist->GetNbinsX();
1202 for (Int_t xBin = nBinsX; xBin >= 1; xBin--) {
1203 if (hist->GetBinContent(xBin, yBin, zBin) > threshold)
1204 return xBin;
1205 }
1206
1207 return -1;
1208}
1209
1210
1211//_____________________________________________________________________________
9d7ad2e4 1212Bool_t AliAnalysisTaskPID::GetParticleFraction(Double_t trackPt, Double_t jetPt, Double_t centralityPercentile,
1213 AliPID::EParticleType species,
e131b05f 1214 Double_t& fraction, Double_t& fractionErrorStat, Double_t& fractionErrorSys) const
1215{
1216 // Computes the particle fraction for the corresponding species for the given trackPt, jetPt and centrality.
1217 // Use jetPt = -1 for inclusive spectra and centralityPercentile = -1 for pp.
1218 // On success (return value kTRUE), fraction contains the particle fraction, fractionErrorStat(Sys) the sigma of its
1219 // statistical (systematic) error
1220
1221 fraction = -999.;
1222 fractionErrorStat = 999.;
1223 fractionErrorSys = 999.;
1224
1225 if (species > AliPID::kProton || species < AliPID::kElectron) {
1226 AliError(Form("Only fractions for species index %d to %d availabe, but not for the requested one: %d", 0, AliPID::kProton, species));
1227 return kFALSE;
1228 }
1229
1230 if (!fFractionHists[species]) {
1231 AliError(Form("Histo with particle fractions for species %d not loaded!", species));
1232
1233 return kFALSE;
1234 }
1235
1236 Int_t jetPtBin = FindBinWithinRange(fFractionHists[species]->GetYaxis(), jetPt);
1237 Int_t centBin = FindBinWithinRange(fFractionHists[species]->GetZaxis(), centralityPercentile);
1238
1239 // The following interpolation takes the bin content of the first/last available bin,
1240 // if requested point lies beyond bin center of first/last bin.
1241 // The interpolation is only done for the x-axis (track pT), i.e. jetPtBin and centBin are fix,
1242 // because the analysis will anyhow be run in bins of jetPt and centrality and
1243 // it is not desired to correlate different jetPt bins via interpolation.
1244
1245 // The same procedure is used for the error of the fraction
1246 TAxis* xAxis = fFractionHists[species]->GetXaxis();
1247
1248 // No interpolation to values beyond the centers of the first/last bins (we don't know exactly where the spectra start or stop,
1249 // thus, search for the first and last bin above 0.0 to constrain the range
1250 Int_t firstBin = TMath::Max(1, FindFirstBinAboveIn3dSubset(fFractionHists[species], 0.0, jetPtBin, centBin));
1251 Int_t lastBin = TMath::Min(fFractionHists[species]->GetNbinsX(),
1252 FindLastBinAboveIn3dSubset(fFractionHists[species], 0.0, jetPtBin, centBin));
1253
1254 if (trackPt <= xAxis->GetBinCenter(firstBin)) {
1255 fraction = fFractionHists[species]->GetBinContent(firstBin, jetPtBin, centBin);
1256 fractionErrorStat = fFractionHists[species]->GetBinError(firstBin, jetPtBin, centBin);
1257 fractionErrorSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(firstBin, jetPtBin, centBin) : 0.;
1258 }
1259 else if (trackPt >= xAxis->GetBinCenter(lastBin)) {
1260 fraction = fFractionHists[species]->GetBinContent(lastBin, jetPtBin, centBin);
1261 fractionErrorStat = fFractionHists[species]->GetBinError(lastBin, jetPtBin, centBin);
1262 fractionErrorSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(lastBin, jetPtBin, centBin) : 0.;
1263 }
1264 else {
1265 Double_t x0 = 0., x1 = 0., y0 = 0., y1 = 0.;
1266 Double_t y0errStat = 0., y1errStat = 0., y0errSys = 0., y1errSys = 0.;
1267 Int_t trackPtBin = xAxis->FindBin(trackPt);
1268
1269 // Linear interpolation between nearest neighbours in trackPt
1270 if (trackPt <= xAxis->GetBinCenter(trackPtBin)) {
1271 y0 = fFractionHists[species]->GetBinContent(trackPtBin - 1, jetPtBin, centBin);
1272 y0errStat = fFractionHists[species]->GetBinError(trackPtBin - 1, jetPtBin, centBin);
1273 y0errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin - 1, jetPtBin, centBin)
1274 : 0.;
1275 x0 = xAxis->GetBinCenter(trackPtBin - 1);
1276 y1 = fFractionHists[species]->GetBinContent(trackPtBin, jetPtBin, centBin);
1277 y1errStat = fFractionHists[species]->GetBinError(trackPtBin, jetPtBin, centBin);
1278 y1errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin, jetPtBin, centBin)
1279 : 0.;
1280 x1 = xAxis->GetBinCenter(trackPtBin);
1281 }
1282 else {
1283 y0 = fFractionHists[species]->GetBinContent(trackPtBin, jetPtBin, centBin);
1284 y0errStat = fFractionHists[species]->GetBinError(trackPtBin, jetPtBin, centBin);
1285 y0errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin, jetPtBin, centBin)
1286 : 0.;
1287 x0 = xAxis->GetBinCenter(trackPtBin);
1288 y1 = fFractionHists[species]->GetBinContent(trackPtBin + 1, jetPtBin, centBin);
1289 y1errStat = fFractionHists[species]->GetBinError(trackPtBin + 1, jetPtBin, centBin);
1290 y1errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin + 1, jetPtBin, centBin)
1291 : 0.;
1292 x1 = xAxis->GetBinCenter(trackPtBin + 1);
1293 }
1294
1295 // Per construction: x0 < trackPt < x1
1296 fraction = y0 + (trackPt - x0) * ((y1 - y0) / (x1 - x0));
1297 fractionErrorStat = y0errStat + (trackPt - x0) * ((y1errStat - y0errStat) / (x1 - x0));
1298 fractionErrorSys = fFractionSysErrorHists[species] ? (y0errSys + (trackPt - x0) * ((y1errSys - y0errSys) / (x1 - x0))) : 0.;
1299 }
1300
1301 return kTRUE;
1302}
1303
1304
1305//_____________________________________________________________________________
9d7ad2e4 1306Bool_t AliAnalysisTaskPID::GetParticleFractions(Double_t trackPt, Double_t jetPt, Double_t centralityPercentile,
1307 Double_t* prob, Int_t smearSpeciesByError,
1308 Int_t takeIntoAccountSpeciesSysError, Bool_t uniformSystematicError) const
e131b05f 1309{
1310 // Fills the particle fractions for the given trackPt, jetPt and centrality into "prob".
1311 // Use jetPt = -1 for inclusive spectra and centralityPercentile = -1 for pp.
1312 // If smearSpeciesByError is >= 0 && < AliPID::kSPECIES, the returned fractions will be a random number distributed
1313 // with a gauss with mean being the corresponding particle fraction and sigma it's error for the considered species
1314 // "smearSpeciesByError".
1315 // Note that in this case the fractions for all species will NOT sum up to 1!
1316 // Thus, all other species fractions will be re-scaled weighted with their corresponding statistical error.
1317 // A similar procedure is used for "takeIntoAccountSpeciesSysError": The systematic error of the corresponding species
1318 // is used to generate a random number with uniform distribution in [mean - sysError, mean + sysError] for the new mean
1319 // (in cace of uniformSystematicError = kTRUE, otherwise it will be a gaus(mean, sysError)),
1320 // then the other species will be re-scaled according to their systematic errors.
1321 // First, the systematic error uncertainty procedure will be performed (that is including re-scaling), then the statistical
1322 // uncertainty procedure.
1323 // On success, kTRUE is returned.
1324
1325 if (!prob || smearSpeciesByError >= AliPID::kSPECIES || takeIntoAccountSpeciesSysError >= AliPID::kSPECIES)
1326 return kFALSE;
1327
1328 Double_t probTemp[AliPID::kSPECIES];
1329 Double_t probErrorStat[AliPID::kSPECIES];
1330 Double_t probErrorSys[AliPID::kSPECIES];
1331
1332 Bool_t success = kTRUE;
1333 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kElectron,
1334 probTemp[AliPID::kElectron], probErrorStat[AliPID::kElectron],
1335 probErrorSys[AliPID::kElectron]);
1336 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kMuon,
1337 probTemp[AliPID::kMuon], probErrorStat[AliPID::kMuon], probErrorSys[AliPID::kMuon]);
1338 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kPion,
1339 probTemp[AliPID::kPion], probErrorStat[AliPID::kPion], probErrorSys[AliPID::kPion]);
1340 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kKaon,
1341 probTemp[AliPID::kKaon], probErrorStat[AliPID::kKaon], probErrorSys[AliPID::kKaon]);
1342 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kProton,
1343 probTemp[AliPID::kProton], probErrorStat[AliPID::kProton], probErrorSys[AliPID::kProton]);
1344
1345 if (!success)
1346 return kFALSE;
1347
1348 // If desired, take into account the systematic error of the corresponding species and re-generate probTemp accordingly
1349 if (takeIntoAccountSpeciesSysError >= 0) {
1350 // Generate random fraction of the considered species "smearSpeciesByError" according to mean and sigma
1351 Double_t generatedFraction = uniformSystematicError
1352 ? fRandom->Rndm() * 2. * probErrorSys[takeIntoAccountSpeciesSysError]
1353 - probErrorSys[takeIntoAccountSpeciesSysError]
1354 + probTemp[takeIntoAccountSpeciesSysError]
1355 : fRandom->Gaus(probTemp[takeIntoAccountSpeciesSysError],
1356 probErrorSys[takeIntoAccountSpeciesSysError]);
1357
1358 // Catch cases with invalid fraction (can happen for large errors), i.e. fraction < 0 or > 1
1359 if (generatedFraction < 0.)
1360 generatedFraction = 0.;
1361 else if (generatedFraction > 1.)
1362 generatedFraction = 1.;
1363
1364 // Calculate difference from original fraction (original fractions sum up to 1!)
1365 Double_t deltaFraction = generatedFraction - probTemp[takeIntoAccountSpeciesSysError];
1366
1367 // Fractions must (including errors) lie inside [0,1] -> Adapt weights accordingly by setting the errors
1368 if (deltaFraction > 0) {
1369 // Some part will be SUBTRACTED from the other fractions
1370 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1371 if (probTemp[i] - probErrorSys[i] < 0)
1372 probErrorSys[i] = probTemp[i];
1373 }
1374 }
1375 else {
1376 // Some part will be ADDED to the other fractions
1377 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1378 if (probTemp[i] + probErrorSys[i] > 1)
1379 probErrorSys[i] = 1. - probTemp[i];
1380 }
1381 }
1382
1383 // Compute summed weight of all fractions except for the considered one
1384 Double_t summedWeight = 0.;
1385 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1386 if (i != takeIntoAccountSpeciesSysError)
1387 summedWeight += probErrorSys[i];
1388 }
1389
1390 // Compute the weight for the other species
1391 /*
1392 if (summedWeight <= 1e-13) {
1393 // If this happens for some reason (it should not!), just assume flat weight
1394 printf("Error: summedWeight (sys error) ~ 0 for trackPt %f, jetPt %f, centralityPercentile %f. Setting flat weight!\n",
1395 trackPt, jetPt, centralityPercentile);
1396 }*/
1397
1398 Double_t weight[AliPID::kSPECIES];
1399 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1400 if (i != takeIntoAccountSpeciesSysError) {
1401 if (summedWeight > 1e-13)
1402 weight[i] = probErrorSys[i] / summedWeight;
1403 else
1404 weight[i] = probErrorSys[i] / (AliPID::kSPECIES - 1);
1405 }
1406 }
1407
1408 // For the final generated fractions, set the generated value for the considered species
1409 // and the generated value minus delta times statistical weight
1410 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1411 if (i != takeIntoAccountSpeciesSysError)
1412 probTemp[i] = probTemp[i] - weight[i] * deltaFraction;
1413 else
1414 probTemp[i] = generatedFraction;
1415 }
1416 }
1417
1418 // Using the values of probTemp (either the original ones or those after taking into account the systematic error),
1419 // calculate the final fractions - if the statistical error is to be taken into account, smear the corresponding
1420 // fraction. If not, just write probTemp to the final result array.
1421 if (smearSpeciesByError >= 0) {
1422 // Generate random fraction of the considered species "smearSpeciesByError" according to mean and sigma
1423 Double_t generatedFraction = fRandom->Gaus(probTemp[smearSpeciesByError], probErrorStat[smearSpeciesByError]);
1424
1425 // Catch cases with invalid fraction (can happen for large errors), i.e. fraction < 0 or > 1
1426 if (generatedFraction < 0.)
1427 generatedFraction = 0.;
1428 else if (generatedFraction > 1.)
1429 generatedFraction = 1.;
1430
1431 // Calculate difference from original fraction (original fractions sum up to 1!)
1432 Double_t deltaFraction = generatedFraction - probTemp[smearSpeciesByError];
1433
1434 // Fractions must (including errors) lie inside [0,1] -> Adapt weights accordingly by setting the errors
1435 if (deltaFraction > 0) {
1436 // Some part will be SUBTRACTED from the other fractions
1437 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1438 if (probTemp[i] - probErrorStat[i] < 0)
1439 probErrorStat[i] = probTemp[i];
1440 }
1441 }
1442 else {
1443 // Some part will be ADDED to the other fractions
1444 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1445 if (probTemp[i] + probErrorStat[i] > 1)
1446 probErrorStat[i] = 1. - probTemp[i];
1447 }
1448 }
1449
1450 // Compute summed weight of all fractions except for the considered one
1451 Double_t summedWeight = 0.;
1452 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1453 if (i != smearSpeciesByError)
1454 summedWeight += probErrorStat[i];
1455 }
1456
1457 // Compute the weight for the other species
1458 /*
1459 if (summedWeight <= 1e-13) {
1460 // If this happens for some reason (it should not!), just assume flat weight
1461 printf("Error: summedWeight (stat error) ~ 0 for trackPt %f, jetPt %f, centralityPercentile %f. Setting flat weight!\n",
1462 trackPt, jetPt, centralityPercentile);
1463 }*/
1464
1465 Double_t weight[AliPID::kSPECIES];
1466 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1467 if (i != smearSpeciesByError) {
1468 if (summedWeight > 1e-13)
1469 weight[i] = probErrorStat[i] / summedWeight;
1470 else
1471 weight[i] = probErrorStat[i] / (AliPID::kSPECIES - 1);
1472 }
1473 }
1474
1475 // For the final generated fractions, set the generated value for the considered species
1476 // and the generated value minus delta times statistical weight
1477 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1478 if (i != smearSpeciesByError)
1479 prob[i] = probTemp[i] - weight[i] * deltaFraction;
1480 else
1481 prob[i] = generatedFraction;
1482 }
1483 }
1484 else {
1485 // Just take the generated values
1486 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
1487 prob[i] = probTemp[i];
1488 }
1489
1490
1491 // Should already be normalised, but make sure that it really is:
1492 Double_t probSum = 0.;
1493 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1494 probSum += prob[i];
1495 }
1496
1497 if (probSum <= 0)
1498 return kFALSE;
1499
1500 if (TMath::Abs(probSum - 1.0) > 1e-4) {
1501 printf("Warning: Re-normalising sum of fractions: Sum is %e\n", probSum);
1502 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1503 prob[i] /= probSum;
1504 }
1505 }
1506
1507 return kTRUE;
1508}
1509
1510
1511//_____________________________________________________________________________
9d7ad2e4 1512const TH3D* AliAnalysisTaskPID::GetParticleFractionHisto(Int_t species, Bool_t sysError) const
e131b05f 1513{
1514 if (species < AliPID::kElectron || species > AliPID::kProton)
1515 return 0x0;
1516
1517 return sysError ? fFractionSysErrorHists[species] : fFractionHists[species];
1518}
1519
1520
1521//_____________________________________________________________________________
1522Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(Int_t motherPDG, Double_t motherGenPt)
1523{
1524 // Strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
1525 // -> Based on function in PWGJE/AliAnalysisTaskFragmentationFunction, which uses
1526 // the following data from CMS pp @ 7 TeV inclusive (JHEP 05 (2011) 064)
1527
1528 Double_t fac = 1;
1529
1530 const Int_t absMotherPDG = TMath::Abs(motherPDG);
1531
1532 if (absMotherPDG == 310 || absMotherPDG == 321) { // K0s / K+ / K-
1533 if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.768049;
1534 else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.732933;
1535 else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.650298;
1536 else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.571332;
1537 else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.518734;
1538 else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.492543;
1539 else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.482704;
1540 else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.488056;
1541 else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.488861;
1542 else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.492862;
1543 else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.504332;
1544 else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.501858;
1545 else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.512970;
1546 else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.524131;
1547 else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.539130;
1548 else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.554101;
1549 else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.560348;
1550 else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.568869;
1551 else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.583310;
1552 else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.604818;
1553 else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.632630;
1554 else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.710070;
1555 else if(6.00 <= motherGenPt && motherGenPt < 8.00) fac = 0.736365;
1556 else if(8.00 <= motherGenPt && motherGenPt < 10.00) fac = 0.835865;
1557 }
1558
1559 if (absMotherPDG == 3122) { // Lambda
1560 if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.645162;
1561 else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.627431;
1562 else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.457136;
1563 else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.384369;
1564 else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.330597;
1565 else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.309571;
1566 else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.293620;
1567 else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.283709;
1568 else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.282047;
1569 else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.277261;
1570 else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.275772;
1571 else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.280726;
1572 else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.288540;
1573 else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.288315;
1574 else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.296619;
1575 else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.302993;
1576 else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.338121;
1577 else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.349800;
1578 else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.356802;
1579 else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.391202;
1580 else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.422573;
1581 else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.573815;
1582 else if(6.00 <= motherGenPt && motherGenPt < 8.00) fac = 0.786984;
1583 else if(8.00 <= motherGenPt && motherGenPt < 10.00) fac = 1.020021;
1584 }
1585
1586 if (absMotherPDG == 3312 || absMotherPDG == 3322) { // xi
1587 if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.666620;
1588 else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.575908;
1589 else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.433198;
1590 else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.340901;
1591 else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.290896;
1592 else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.236074;
1593 else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.218681;
1594 else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.207763;
1595 else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.222848;
1596 else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.208806;
1597 else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.197275;
1598 else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.183645;
1599 else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.188788;
1600 else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.188282;
1601 else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.207442;
1602 else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.240388;
1603 else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.241916;
1604 else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.208276;
1605 else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.234550;
1606 else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.251689;
1607 else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.310204;
1608 else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.343492;
1609 }
1610
1611 const Double_t weight = 1. / fac;
1612
1613 return weight;
1614}
1615
1616
1617//_____________________________________________________________________________
1618Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(AliMCEvent* mcEvent, AliMCParticle* daughter)
1619{
1620 // Strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
1621 // -> Based on function in PWGJE/AliAnalysisTaskFragmentationFunction
1622
1623 if (!mcEvent)
1624 return 1.;
1625
1626 AliMCParticle* currentMother = daughter;
1627 AliMCParticle* currentDaughter = daughter;
1628
1629
1630 // find first primary mother K0s, Lambda or Xi
1631 while(1) {
1632 Int_t daughterPDG = currentDaughter->PdgCode();
1633
1634 Int_t motherLabel = currentDaughter->GetMother();
1635 if(motherLabel >= mcEvent->GetNumberOfTracks()){ // protection
1636 currentMother = currentDaughter;
1637 break;
1638 }
1639
1640 currentMother = (AliMCParticle*)mcEvent->GetTrack(motherLabel);
1641
1642 if (!currentMother) {
1643 currentMother = currentDaughter;
1644 break;
1645 }
1646
1647 Int_t motherPDG = currentMother->PdgCode();
1648
1649 // phys. primary found ?
1650 if (mcEvent->IsPhysicalPrimary(motherLabel))
1651 break;
1652
1653 if (TMath::Abs(daughterPDG) == 321) {
1654 // K+/K- e.g. from phi (ref data not feeddown corrected)
1655 currentMother = currentDaughter;
1656 break;
1657 }
1658 if (TMath::Abs(motherPDG) == 310) {
1659 // K0s e.g. from phi (ref data not feeddown corrected)
1660 break;
1661 }
1662 if (TMath::Abs(motherPDG) == 3212 && TMath::Abs(daughterPDG) == 3122) {
1663 // Mother Sigma0, daughter Lambda (this case not included in feeddown corr.)
1664 currentMother = currentDaughter;
1665 break;
1666 }
1667
1668 currentDaughter = currentMother;
1669 }
1670
1671
1672 Int_t motherPDG = currentMother->PdgCode();
1673 Double_t motherGenPt = currentMother->Pt();
1674
1675 return GetMCStrangenessFactorCMS(motherPDG, motherGenPt);
1676}
1677
1678
1679// _________________________________________________________________________________
1680Bool_t AliAnalysisTaskPID::IsSecondaryWithStrangeMotherMC(AliMCEvent* mcEvent, Int_t partLabel)
1681{
1682 // Check whether particle is a secondary with strange mother, i.e. returns kTRUE if a strange mother is found
1683 // and the particle is NOT a physical primary. In all other cases kFALSE is returned
1684
1685 if (!mcEvent || partLabel < 0)
1686 return kFALSE;
1687
1688 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(partLabel);
1689
1690 if (!part)
1691 return kFALSE;
1692
1693 if (mcEvent->IsPhysicalPrimary(partLabel))
1694 return kFALSE;
1695
1696 Int_t iMother = part->GetMother();
1697 if (iMother < 0)
1698 return kFALSE;
1699
1700
1701 AliMCParticle* partM = (AliMCParticle*)mcEvent->GetTrack(iMother);
1702 if (!partM)
1703 return kFALSE;
1704
1705 Int_t codeM = TMath::Abs(partM->PdgCode());
1706 Int_t mfl = Int_t(codeM / TMath::Power(10, Int_t(TMath::Log10(codeM))));
1707 if (mfl == 3 && codeM != 3) // codeM = 3 is for s quark
1708 return kTRUE;
1709
1710 return kFALSE;
1711}
1712
1713
1714//_____________________________________________________________________________
9d7ad2e4 1715Bool_t AliAnalysisTaskPID::SetParticleFractionHisto(const TH3D* hist, Int_t species, Bool_t sysError)
e131b05f 1716{
1717 // Store a clone of hist (containing the particle fractions of the corresponding species with statistical error (sysError = kFALSE)
1718 // or systematic error (sysError = kTRUE), respectively), internally
1719
1720 if (species < AliPID::kElectron || species > AliPID::kProton) {
1721 AliError(Form("Only fractions for species index %d to %d can be set, but not for the requested one: %d", 0,
1722 AliPID::kProton, species));
1723 return kFALSE;
1724 }
1725
1726 if (sysError) {
1727 delete fFractionSysErrorHists[species];
1728
1729 fFractionSysErrorHists[species] = new TH3D(*hist);
1730 }
1731 else {
1732 delete fFractionHists[species];
1733
1734 fFractionHists[species] = new TH3D(*hist);
1735 }
1736
1737 return kTRUE;
1738}
1739
1740
1741//_____________________________________________________________________________
9d7ad2e4 1742Bool_t AliAnalysisTaskPID::SetParticleFractionHistosFromFile(const TString filePathName, Bool_t sysError)
e131b05f 1743{
1744 // Loads particle fractions for all species from the desired file and returns kTRUE on success.
1745 // The maps are assumed to be of Type TH3D, to sit in the main directory and to have names
1746 // Form("hFraction_%e", AliPID::ParticleName(i)) for sysError = kFALSE and
1747 // Form("hFractionSysError_%e", AliPID::ParticleName(i)) for sysError = kTRUE.
1748
1749 TFile* f = TFile::Open(filePathName.Data());
1750 if (!f) {
1751 std::cout << "Failed to open file with particle fractions \"" << filePathName.Data() << "\"!" << std::endl;
1752 return kFALSE;
1753 }
1754
1755 TH3D* hist = 0x0;
1756 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1757 TString histName = Form("hFraction%s_%s", sysError ? "SysError" : "", AliPID::ParticleName(i));
1758 hist = dynamic_cast<TH3D*>(f->Get(histName.Data()));
1759 if (!hist) {
1760 std::cout << "Failed to load particle fractions for " << histName.Data() << "!";
1761 std::cout << std::endl << "Cleaning up particle fraction histos!" << std::endl;
1762 CleanupParticleFractionHistos();
1763 return kFALSE;
1764 }
1765
1766 if (!SetParticleFractionHisto(hist, i, sysError)) {
1767 std::cout << "Failed to load particle fractions for " << histName.Data() << "!";
1768 std::cout << std::endl << "Cleaning up particle fraction histos!" << std::endl;
1769 CleanupParticleFractionHistos();
1770 return kFALSE;
1771 }
1772 }
1773
1774 delete hist;
1775
1776 return kTRUE;
1777
1778}
1779
1780
1781//_____________________________________________________________________________
9d7ad2e4 1782Int_t AliAnalysisTaskPID::GetRandomParticleTypeAccordingToParticleFractions(Double_t trackPt, Double_t jetPt,
1783 Double_t centralityPercentile,
1784 Bool_t smearByError,
1785 Bool_t takeIntoAccountSysError) const
e131b05f 1786{
1787 // Uses the stored histograms with the particle fractions to generate a random particle type according to these fractions.
1788 // In case of problems (e.g. histo missing), AliPID::kUnknown is returned.
1789 // If smearByError is kTRUE, the used fractions will be random numbers distributed with a gauss with mean
1790 // being the corresponding particle fraction and sigma it's error.
1791 // Note that in this case only the fraction of a random species is varied in this way. The other fractions
1792 // will be re-normalised according their statistical errors.
1793 // The same holds for the systematic error of species "takeIntoAccountSpeciesSysError", but the random number will be
1794 // uniformly distributed within [mean - sys, mean + sys] and the re-normalisation will be weighted with the systematic errors.
1795 // Note that the fractions will be calculated first with only the systematic error taken into account (if desired), including
1796 // re-normalisation. Then, the resulting fractions will be used to calculate the final fractions - either with statistical error
1797 // or without. The species, for which the error will be used for smearing, is the same for sys and stat error.
1798
1799 Double_t prob[AliPID::kSPECIES];
1800 Int_t randomSpecies = (smearByError || takeIntoAccountSysError) ? (Int_t)(fRandom->Rndm() * AliPID::kSPECIES) : -1;
1801 Bool_t success = GetParticleFractions(trackPt, jetPt, centralityPercentile, prob, randomSpecies, randomSpecies);
1802
1803 if (!success)
1804 return AliPID::kUnknown;
1805
1806 Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
1807
1808 if (rnd <= prob[AliPID::kPion])
1809 return AliPID::kPion;
1810 else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon])
1811 return AliPID::kKaon;
1812 else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon] + prob[AliPID::kProton])
1813 return AliPID::kProton;
1814 else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon] + prob[AliPID::kProton] + prob[AliPID::kElectron])
1815 return AliPID::kElectron;
1816
1817 return AliPID::kMuon; //else it must be a muon (only species left)
1818}
1819
1820
1821//_____________________________________________________________________________
9d7ad2e4 1822AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::GenerateDetectorResponse(AliAnalysisTaskPID::ErrorCode errCode,
1823 Double_t mean, Double_t sigma,
1824 Double_t* responses, Int_t nResponses,
1825 Bool_t usePureGaus)
e131b05f 1826{
1827 // Generate detector response. If a previous generation was not successful or there is something wrong with this signal generation,
1828 // the function will return kFALSE
1829 if (!responses)
1830 return kError;
1831
1832 // Reset response array
1833 for (Int_t i = 0; i < nResponses; i++)
1834 responses[i] = -999;
1835
1836 if (errCode == kError)
1837 return kError;
1838
1839 ErrorCode ownErrCode = kNoErrors;
1840
1841 if (fUseConvolutedGaus && !usePureGaus) {
1842 // In case of convoluted gauss, calculate the probability density only once to save a lot of time!
1843
1844 TH1* hProbDensity = 0x0;
1845 ownErrCode = SetParamsForConvolutedGaus(mean, sigma);
1846 if (ownErrCode == kError)
1847 return kError;
1848
1849 hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
1850
1851 for (Int_t i = 0; i < nResponses; i++) {
1852 responses[i] = hProbDensity->GetRandom();
1853 //responses[i] fConvolutedGausDeltaPrime->GetRandom(); // MUCH slower than using the binned version via the histogram
1854 }
1855 }
1856 else {
1857 for (Int_t i = 0; i < nResponses; i++) {
1858 responses[i] = fRandom->Gaus(mean, sigma);
1859 }
1860 }
1861
1862 // If forwarded error code was a warning (error case has been handled before), return a warning
1863 if (errCode == kWarning)
1864 return kWarning;
1865
1866 return ownErrCode; // Forward success/warning
1867}
1868
1869
1870//_____________________________________________________________________________
1871void AliAnalysisTaskPID::PrintSettings(Bool_t printSystematicsSettings) const
1872{
1873 // Print current settings.
1874
1875 printf("\n\nSettings for task %s:\n", GetName());
1876 printf("Is pPb/Pbp: %d -> %s\n", GetIsPbpOrpPb(), GetIsPbpOrpPb() ? "Adapting vertex cuts" : "Using standard vertex cuts");
1877 printf("Track cuts: %s\n", fTrackFilter ? fTrackFilter->GetTitle() : "-");
1878 printf("Eta cut: %.2f <= |eta| <= %.2f\n", GetEtaAbsCutLow(), GetEtaAbsCutUp());
1879 printf("Phi' cut: %d\n", GetUsePhiCut());
a6852ea8 1880 printf("TPCCutMIGeo: %d\n", GetUseTPCCutMIGeo());
493982d9
ML		 1881 if (GetUseTPCCutMIGeo()) {
1882 printf("GetCutGeo: %f\n", GetCutGeo());
1883 printf("GetCutNcr: %f\n", GetCutNcr());
1884 printf("GetCutNcl: %f\n", GetCutNcl());
1885 }
1886 printf("TPCnclCut: %d\n", GetUseTPCnclCut());
1887 if (GetUseTPCnclCut()) {
1888 printf("GetCutPureNcl: %d\n", GetCutPureNcl());
1889 }
e131b05f 1890
1891 printf("\n");
1892
1893 printf("Centrality estimator: %s\n", GetCentralityEstimator().Data());
1894
1895 printf("\n");
1896
1897 printf("Use MC-ID for signal generation: %d\n", GetUseMCidForGeneration());
1898 printf("Use ITS: %d\n", GetUseITS());
1899 printf("Use TOF: %d\n", GetUseTOF());
1900 printf("Use priors: %d\n", GetUsePriors());
1901 printf("Use TPC default priors: %d\n", GetUseTPCDefaultPriors());
1902 printf("Use convoluted Gauss: %d\n", GetUseConvolutedGaus());
1903 printf("Accuracy of non-Gaussian tail: %e\n", GetAccuracyNonGaussianTail());
1904 printf("Take into account muons: %d\n", GetTakeIntoAccountMuons());
1905 printf("\nParams for transition from gauss to asymmetric shape:\n");
1906 printf("[0]: %e\n", GetConvolutedGaussTransitionPar(0));
1907 printf("[1]: %e\n", GetConvolutedGaussTransitionPar(1));
1908 printf("[2]: %e\n", GetConvolutedGaussTransitionPar(2));
1909
9e95a906 1910 printf("\n");
1911
1912 printf("Do PID: %d\n", fDoPID);
1913 printf("Do Efficiency: %d\n", fDoEfficiency);
1914 printf("Do PtResolution: %d\n", fDoPtResolution);
1915
e131b05f 1916 printf("\n");
1917
1918 printf("Input from other task: %d\n", GetInputFromOtherTask());
1919 printf("Store additional jet information: %d\n", GetStoreAdditionalJetInformation());
1920 printf("Store centrality percentile: %d", GetStoreCentralityPercentile());
1921
1922 if (printSystematicsSettings)
1923 PrintSystematicsSettings();
1924 else
1925 printf("\n\n\n");
1926}
1927
1928
1929//_____________________________________________________________________________
1930void AliAnalysisTaskPID::PrintSystematicsSettings() const
1931{
1932 // Print current settings for systematic studies.
1933
1934 printf("\n\nSettings for systematics for task %s:\n", GetName());
1935 printf("Splines:\t%f\n", GetSystematicScalingSplines());
1936 printf("EtaCorrMomThr:\t%f\n", GetSystematicScalingEtaCorrectionMomentumThr());
1937 printf("EtaCorrLowP:\t%f\n", GetSystematicScalingEtaCorrectionLowMomenta());
1938 printf("EtaCorrHighP:\t%f\n", GetSystematicScalingEtaCorrectionHighMomenta());
1939 printf("SigmaPara:\t%f\n", GetSystematicScalingEtaSigmaPara());
1940 printf("MultCorr:\t%f\n", GetSystematicScalingMultCorrection());
1941
1942 printf("\n\n");
1943}
1944
1945
1946//_____________________________________________________________________________
1947Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePDGcode, Double_t centralityPercentile,
1948 Double_t jetPt)
1949{
1950 // Process the track (generate expected response, fill histos, etc.).
1951 // particlePDGcode == 0 means data. Otherwise, the corresponding MC ID will be assumed.
1952
1953 //Printf("Debug: Task %s is starting to process track: dEdx %f, pTPC %f, eta %f, ncl %d\n", GetName(), track->GetTPCsignal(), track->GetTPCmomentum(),
1954 // track->Eta(), track->GetTPCsignalN());
1955
9e95a906 1956 if(fDebug > 1)
1957 printf("File: %s, Line: %d: ProcessTrack\n", (char*)__FILE__, __LINE__);
1958
1959 if (!fDoPID)
1960 return kFALSE;
1961
1962 if(fDebug > 2)
1963 printf("File: %s, Line: %d: ProcessTrack -> Processing started\n", (char*)__FILE__, __LINE__);
1964
e131b05f 1965 const Bool_t isMC = (particlePDGcode == 0) ? kFALSE : kTRUE;
1966
1967 Int_t binMC = -1;
1968
1969 if (isMC) {
1970 if (TMath::Abs(particlePDGcode) == 211) {//Pion
1971 binMC = 3;
1972 }
1973 else if (TMath::Abs(particlePDGcode) == 321) {//Kaon
1974 binMC = 1;
1975 }
1976 else if (TMath::Abs(particlePDGcode) == 2212) {//Proton
1977 binMC = 4;
1978 }
1979 else if (TMath::Abs(particlePDGcode) == 11) {//Electron
1980 binMC = 0;
1981 }
1982 else if (TMath::Abs(particlePDGcode) == 13) {//Muon
1983 binMC = 2;
1984 }
1985 else // In MC-ID case, set to underflow bin such that the response from this track is only used for unidentified signal generation
1986 // or signal generation with PID response and the track is still there (as in data) - e.g. relevant w.r.t. deuterons.
1987 // This is important to be as much as possible consistent with data. And the tracks can still be removed by disabling the
1988 // underflow bin for the projections
1989 binMC = -1;
1990 }
1991
1992 // Momenta
1993 //Double_t p = track->GetP();
1994 //Double_t pTPC = track->GetTPCmomentum();
1995 Double_t pT = track->Pt();
1996
1997 Double_t z = -1, xi = -1;
1998 GetJetTrackObservables(pT, jetPt, z, xi);
1999
2000
2001 Double_t trackCharge = track->Charge();
2002
2003 // TPC signal
2004 Double_t dEdxTPC = fPIDResponse->IsTunedOnData() ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
2005
2006 if (dEdxTPC <= 0) {
2007 Printf("Skipping track with strange dEdx value: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), track->GetTPCmomentum(),
2008 track->Eta(), track->GetTPCsignalN());
2009 return kFALSE;
2010 }
2011
2012
2013
2014
2015 Double_t dEdxEl, dEdxKa, dEdxPi, dEdxMu, dEdxPr;
2016 Double_t sigmaEl, sigmaKa, sigmaPi, sigmaMu, sigmaPr;
2017
2018 if (fDoAnySystematicStudiesOnTheExpectedSignal) {
2019 // Get the uncorrected signal first and the corresponding correction factors.
2020 // Then modify the correction factors and properly recalculate the corrected dEdx
2021
2022 // Get pure spline values for dEdx_expected, without any correction
2023 dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2024 dEdxKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2025 dEdxPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2026 dEdxMu = !fTakeIntoAccountMuons ? -1 :
2027 fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2028 dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2029
2030 // Scale splines, if desired
2031 if (TMath::Abs(fSystematicScalingSplines - 1.0) > fgkEpsilon) {
2032 dEdxEl *= fSystematicScalingSplines;
2033 dEdxKa *= fSystematicScalingSplines;
2034 dEdxPi *= fSystematicScalingSplines;
2035 dEdxMu *= fTakeIntoAccountMuons ? fSystematicScalingSplines : 1.;
2036 dEdxPr *= fSystematicScalingSplines;
2037 }
2038
2039 // Get the eta correction factors for the (modified) expected dEdx
2040 Double_t etaCorrEl = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxEl) : 1.;
2041 Double_t etaCorrKa = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxKa) : 1.;
2042 Double_t etaCorrPi = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPi) : 1.;
2043 Double_t etaCorrMu = fTakeIntoAccountMuons && !fPIDResponse->UseTPCEtaCorrection() ?
2044 fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxMu) : 1.;
2045 Double_t etaCorrPr = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPr) : 1.;
2046
2047 // Scale eta correction factors, if desired (and eta correction maps are to be used, otherwise it is not possible!)
2048 if (fPIDResponse->UseTPCEtaCorrection() &&
2049 (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon ||
2050 TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon)) {
2051 // Since we do not want to scale the splines with this, but only the eta variation, only scale the deviation of the correction factor!
2052 // 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!
2053
2054
2055 // Due to additional azimuthal effects, there is an additional eta dependence for low momenta which is not corrected successfully so far.
2056 // One can assign a different (higher) systematic scale factor for this low-p region and a threshold which separates low- and high-p.
2057 // 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
2058 Double_t usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionHighMomenta;
2059
2060 if (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - fSystematicScalingEtaCorrectionLowMomenta) > fgkEpsilon) {
2061 const Double_t pTPC = track->GetTPCmomentum();
2062 const Double_t fractionHighMomentumScaleFactor = 0.5 * (1. + TMath::Erf((pTPC - fSystematicScalingEtaCorrectionMomentumThr) / 0.1));
2063 usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionLowMomenta * (1 - fractionHighMomentumScaleFactor)
2064 + fSystematicScalingEtaCorrectionHighMomenta * fractionHighMomentumScaleFactor;
2065 }
2066
2067 etaCorrEl = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrEl - 1.0);
2068 etaCorrKa = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrKa - 1.0);
2069 etaCorrPi = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPi - 1.0);
2070 etaCorrMu = fTakeIntoAccountMuons ? (1.0 + usedSystematicScalingEtaCorrection * (etaCorrMu - 1.0)) : 1.0;
2071 etaCorrPr = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPr - 1.0);
2072 }
2073
2074 // Get the multiplicity correction factors for the (modified) expected dEdx
2075 const Int_t currEvtMultiplicity = fPIDResponse->GetTPCResponse().GetCurrentEventMultiplicity();
2076
2077 Double_t multiplicityCorrEl = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2078 dEdxEl * etaCorrEl, currEvtMultiplicity) : 1.;
2079 Double_t multiplicityCorrKa = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2080 dEdxKa * etaCorrKa, currEvtMultiplicity) : 1.;
2081 Double_t multiplicityCorrPi = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2082 dEdxPi * etaCorrPi, currEvtMultiplicity) : 1.;
2083 Double_t multiplicityCorrMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCMultiplicityCorrection() ?
2084 fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track, dEdxMu * etaCorrMu, currEvtMultiplicity) : 1.;
2085 Double_t multiplicityCorrPr = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2086 dEdxPr * etaCorrPr, currEvtMultiplicity) : 1.;
2087
2088 Double_t multiplicityCorrSigmaEl = fPIDResponse->UseTPCMultiplicityCorrection() ?
2089 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxEl * etaCorrEl, currEvtMultiplicity) : 1.;
2090 Double_t multiplicityCorrSigmaKa = fPIDResponse->UseTPCMultiplicityCorrection() ?
2091 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxKa * etaCorrKa, currEvtMultiplicity) : 1.;
2092 Double_t multiplicityCorrSigmaPi = fPIDResponse->UseTPCMultiplicityCorrection() ?
2093 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxPi * etaCorrPi, currEvtMultiplicity) : 1.;
2094 Double_t multiplicityCorrSigmaMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCMultiplicityCorrection() ?
2095 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxMu * etaCorrMu, currEvtMultiplicity) : 1.;
2096 Double_t multiplicityCorrSigmaPr = fPIDResponse->UseTPCMultiplicityCorrection() ?
2097 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxPr * etaCorrPr, currEvtMultiplicity) : 1.;
2098
2099 // Scale multiplicity correction factors, if desired (and multiplicity correction functions are to be used, otherwise it is not possible!)
2100 if (fPIDResponse->UseTPCMultiplicityCorrection() && TMath::Abs(fSystematicScalingMultCorrection - 1.0) > fgkEpsilon) {
2101 // Since we do not want to scale the splines with this, but only the multiplicity variation, only scale the deviation of the correction factor!
2102 // 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!
2103
2104 multiplicityCorrEl = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrEl - 1.0);
2105 multiplicityCorrKa = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrKa - 1.0);
2106 multiplicityCorrPi = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrPi - 1.0);
2107 multiplicityCorrMu = fTakeIntoAccountMuons ? (1.0 + fSystematicScalingMultCorrection * (multiplicityCorrMu - 1.0)) : 1.0;
2108 multiplicityCorrPr = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrPr - 1.0);
2109 }
2110
2111 // eta correction must be enabled in order to use the new sigma parametrisation maps. Since this is the absolute sigma
2112 // for a track calculated with the unscaled paramaters, we have to devide out dEdxExpectedEtaCorrected and then need
2113 // to scale with the multiplicitySigmaCorrFactor * fSystematicScalingEtaSigmaPara. In the end, one has to scale with the
2114 // (modified) dEdx to get the absolute sigma
2115 // This means there is no extra parameter for the multiplicitySigmaCorrFactor, but only for the sigma map itself.
2116 // This is valid, since it appears only as a product. One has to assume a larger systematic shift in case of additional
2117 // multiplicity dependence....
2118 Double_t sigmaRelEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2119 / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2120 * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaEl;
2121
2122 Double_t sigmaRelKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2123 / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2124 * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaKa;
2125
2126 Double_t sigmaRelPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2127 / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2128 * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaPi;
2129
2130 Double_t sigmaRelMu = fTakeIntoAccountMuons ?
2131 fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2132 / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2133 * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaMu
2134 : 999.;
2135
2136 Double_t sigmaRelPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2137 / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2138 * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaPr;
2139
2140 // Now scale the (possibly modified) spline values with the (possibly modified) correction factors
2141 dEdxEl *= etaCorrEl * multiplicityCorrEl;
2142 dEdxKa *= etaCorrKa * multiplicityCorrKa;
2143 dEdxPi *= etaCorrPi * multiplicityCorrPi;
2144 dEdxMu *= etaCorrMu * multiplicityCorrMu;
2145 dEdxPr *= etaCorrPr * multiplicityCorrPr;
2146
2147 // Finally, get the absolute sigma
2148 sigmaEl = sigmaRelEl * dEdxEl;
2149 sigmaKa = sigmaRelKa * dEdxKa;
2150 sigmaPi = sigmaRelPi * dEdxPi;
2151 sigmaMu = sigmaRelMu * dEdxMu;
2152 sigmaPr = sigmaRelPr * dEdxPr;
2153
2154 }
2155 else {
2156 // No systematic studies on expected signal - just take it as it comve from the TPCPIDResponse
2157 dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2158 fPIDResponse->UseTPCEtaCorrection(),
2159 fPIDResponse->UseTPCMultiplicityCorrection());
2160 dEdxKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2161 fPIDResponse->UseTPCEtaCorrection(),
2162 fPIDResponse->UseTPCMultiplicityCorrection());
2163 dEdxPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2164 fPIDResponse->UseTPCEtaCorrection(),
2165 fPIDResponse->UseTPCMultiplicityCorrection());
2166 dEdxMu = !fTakeIntoAccountMuons ? -1 :
2167 fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2168 fPIDResponse->UseTPCEtaCorrection(),
2169 fPIDResponse->UseTPCMultiplicityCorrection());
2170 dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2171 fPIDResponse->UseTPCEtaCorrection(),
2172 fPIDResponse->UseTPCMultiplicityCorrection());
2173
2174 sigmaEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2175 fPIDResponse->UseTPCEtaCorrection(),
2176 fPIDResponse->UseTPCMultiplicityCorrection());
2177 sigmaKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2178 fPIDResponse->UseTPCEtaCorrection(),
2179 fPIDResponse->UseTPCMultiplicityCorrection());
2180 sigmaPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2181 fPIDResponse->UseTPCEtaCorrection(),
2182 fPIDResponse->UseTPCMultiplicityCorrection());
2183 sigmaMu = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2184 fPIDResponse->UseTPCEtaCorrection(),
2185 fPIDResponse->UseTPCMultiplicityCorrection());
2186 sigmaPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2187 fPIDResponse->UseTPCEtaCorrection(),
2188 fPIDResponse->UseTPCMultiplicityCorrection());
2189 }
2190 /*OLD with deltaSpecies
2191 Double_t deltaElectron = dEdxTPC - dEdxEl;
2192 Double_t deltaKaon = dEdxTPC - dEdxKa;
2193 Double_t deltaPion = dEdxTPC - dEdxPi;
2194 Double_t deltaProton = dEdxTPC - dEdxPr;
2195 */
2196
2197 Double_t deltaPrimeElectron = (dEdxEl > 0) ? dEdxTPC / dEdxEl : -1;
2198 if (dEdxEl <= 0) {
2199 Printf("Error: Expected TPC signal <= 0 for electron hypothesis");
2200 return kFALSE;
2201 }
2202
2203 Double_t deltaPrimeKaon = (dEdxKa > 0) ? dEdxTPC / dEdxKa : -1;
2204 if (dEdxKa <= 0) {
2205 Printf("Error: Expected TPC signal <= 0 for kaon hypothesis");
2206 return kFALSE;
2207 }
2208
2209 Double_t deltaPrimePion = (dEdxPi > 0) ? dEdxTPC / dEdxPi : -1;
2210 if (dEdxPi <= 0) {
2211 Printf("Error: Expected TPC signal <= 0 for pion hypothesis");
2212 return kFALSE;
2213 }
2214
2215 Double_t deltaPrimeProton = (dEdxPr > 0) ? dEdxTPC / dEdxPr : -1;
2216 if (dEdxPr <= 0) {
2217 Printf("Error: Expected TPC signal <= 0 for proton hypothesis");
2218 return kFALSE;
2219 }
2220
2221 /*TODO for TOF
2222 // TOF signal
2223 Double_t times[AliPID::kSPECIES];
2224 track->GetIntegratedTimes(times);
2225 AliTOFPIDResponse tofPIDResponse = fPIDResponse->GetTOFResponse();
2226 Float_t electronDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kElectron);
2227 Float_t pionDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kPion);
2228 Float_t kaonDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kKaon);
2229 Float_t protonDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kProton);
2230 */
2231
9e95a906 2232 if(fDebug > 2)
2233 printf("File: %s, Line: %d: ProcessTrack -> Compute probabilities\n", (char*)__FILE__, __LINE__);
e131b05f 2234
2235 // Use probabilities to weigh the response generation for the different species.
2236 // Also determine the most probable particle type.
2237 Double_t prob[AliPID::kSPECIESC];
2238 for (Int_t i = 0; i < AliPID::kSPECIESC; i++)
2239 prob[i] = 0;
2240
2241 fPIDcombined->ComputeProbabilities(track, fPIDResponse, prob);
2242
2243 // Bug: One can set the number of species for PIDcombined, but PIDcombined will call PIDresponse, which writes without testing
2244 // the probs for kSPECIESC (including light nuclei) into the array.
2245 // In this case, when only kSPECIES are considered, the probabilities have to be rescaled!
2246 for (Int_t i = AliPID::kSPECIES; i < AliPID::kSPECIESC; i++)
2247 prob[i] = 0;
2248
2249 // If muons are not to be taken into account, just set their probability to zero and normalise the remaining probabilities
2250 if (!fTakeIntoAccountMuons)
2251 prob[AliPID::kMuon] = 0;
2252
2253 Double_t probSum = 0.;
2254 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2255 probSum += prob[i];
2256
2257 if (probSum > 0) {
2258 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2259 prob[i] /= probSum;
2260 }
2261
2262 if (!isMC) {
2263 // If there is no MC information, take the most probable species for the ID
2264 Float_t max = 0.;
2265 Int_t maxIndex = -1;
2266 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2267 if (prob[i] > max) {
2268 max = prob[i];
2269 maxIndex = i;
2270 }
2271 }
2272
2273 // Translate from AliPID numbering to numbering of this class
2274 if (max > 0) {
2275 if (maxIndex == AliPID::kElectron)
2276 binMC = 0;
2277 else if (maxIndex == AliPID::kKaon)
2278 binMC = 1;
2279 else if (maxIndex == AliPID::kMuon)
2280 binMC = 2;
2281 else if (maxIndex == AliPID::kPion)
2282 binMC = 3;
2283 else if (maxIndex == AliPID::kProton)
2284 binMC = 4;
2285 else
2286 binMC = -1;
2287 }
2288 else {
2289 // Only take track into account for expectation values, if valid pid response is available.. Otherwise: Set to underflow bin.
2290 binMC = -1;
2291 }
2292 }
2293
2294 /*
2295 //For testing: Swap p<->pT to analyse pure p-dependence => Needs to be removed later
2296 Double_t temp = pT;
2297 pT = pTPC;
2298 pTPC = temp;
2299 */
2300
2301
2302 Double_t entry[fStoreAdditionalJetInformation ? kDataNumAxes : kDataNumAxes - fgkNumJetAxes];
2303 entry[kDataMCID] = binMC;
2304 entry[kDataSelectSpecies] = 0;
2305 entry[kDataPt] = pT;
2306 entry[kDataDeltaPrimeSpecies] = deltaPrimeElectron;
2307 entry[kDataCentrality] = centralityPercentile;
2308
2309 if (fStoreAdditionalJetInformation) {
2310 entry[kDataJetPt] = jetPt;
2311 entry[kDataZ] = z;
2312 entry[kDataXi] = xi;
2313 }
2314
2315 entry[GetIndexOfChargeAxisData()] = trackCharge;
2316
2317 /*OLD with TOF, p_TPC_Inner and p_vertex and deltaSpecies
2318 // MC PID, SelectSpecies, P(TPC_inner), pT, p(Vertex), deltaSpecies, deltaPrimeSpecies, deltaTOFspecies
2319 Double_t entry[8] = { binMC, 0, pTPC, pT, p, deltaElectron, deltaPrimeElectron, electronDeltaTOF };
2320 */
2321
2322 fhPIDdataAll->Fill(entry);
2323
2324 entry[kDataSelectSpecies] = 1;
2325 //OLD with deltaSpecies entry[5] = deltaKaon;
2326 entry[kDataDeltaPrimeSpecies] = deltaPrimeKaon;
2327 //TODO for TOF entry[7] = kaonDeltaTOF;
2328 fhPIDdataAll->Fill(entry);
2329
2330 entry[kDataSelectSpecies] = 2;
2331 //OLD with deltaSpecies entry[5] = deltaPion;
2332 entry[kDataDeltaPrimeSpecies] = deltaPrimePion;
2333 //TODO for TOF entry[7] = pionDeltaTOF;
2334 fhPIDdataAll->Fill(entry);
2335
2336 entry[kDataSelectSpecies] = 3;
2337 //OLD with deltaSpecies entry[5] = deltaProton;
2338 entry[kDataDeltaPrimeSpecies] = deltaPrimeProton;
2339 //TODO for TOF entry[7] = protonDeltaTOF;
2340 fhPIDdataAll->Fill(entry);
2341
2342
2343 // Construct the expected shape for the transition p -> pT
2344
2345 Double_t genEntry[fStoreAdditionalJetInformation ? kGenNumAxes : kGenNumAxes - fgkNumJetAxes];
2346 genEntry[kGenMCID] = binMC;
2347 genEntry[kGenSelectSpecies] = 0;
2348 genEntry[kGenPt] = pT;
2349 genEntry[kGenDeltaPrimeSpecies] = -999;
2350 genEntry[kGenCentrality] = centralityPercentile;
2351
2352 if (fStoreAdditionalJetInformation) {
2353 genEntry[kGenJetPt] = jetPt;
2354 genEntry[kGenZ] = z;
2355 genEntry[kGenXi] = xi;
2356 }
2357
2358 genEntry[GetIndexOfChargeAxisGen()] = trackCharge;
2359
2360 //OLD with deltaSpecies Double_t genEntry[5] = { binMC, 0, pT, -999, -999 }; // MC PID, SelectSpecies, pT, deltaSpecies, deltaPrimeSpecies
2361
2362 // Generate numGenEntries "responses" with fluctuations around the expected values.
2363 // The higher the (transverse) momentum, the more "responses" will be generated in order not to run out of statistics too fast.
2364 Int_t numGenEntries = 10;
2365 if (pT >= 5)
2366 numGenEntries *= 5;
2367 else if (pT >= 2)
2368 numGenEntries *= 2;
2369
2370 // Jets have even worse statistics, therefore, scale numGenEntries further
2371 if (jetPt >= 40)
2372 numGenEntries *= 20;
2373 else if (jetPt >= 20)
2374 numGenEntries *= 10;
2375 else if (jetPt >= 10)
2376 numGenEntries *= 2;
2377
2378
2379 // Do not generate more entries than available in memory!
2380 if (numGenEntries > fgkMaxNumGenEntries)// fgkMaxNumGenEntries = 1000
2381 numGenEntries = fgkMaxNumGenEntries;
2382
2383 ErrorCode errCode = kNoErrors;
2384
9e95a906 2385 if(fDebug > 2)
2386 printf("File: %s, Line: %d: ProcessTrack -> Generate Responses\n", (char*)__FILE__, __LINE__);
2387
e131b05f 2388 // Electrons
2389 errCode = GenerateDetectorResponse(errCode, 1., sigmaEl / dEdxEl, fGenRespElDeltaPrimeEl, numGenEntries);
2390 errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxKa, sigmaEl / dEdxKa, fGenRespElDeltaPrimeKa, numGenEntries);
2391 errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxPi, sigmaEl / dEdxPi, fGenRespElDeltaPrimePi, numGenEntries);
2392 errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxPr, sigmaEl / dEdxPr, fGenRespElDeltaPrimePr, numGenEntries);
2393
2394 // Kaons
2395 errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxEl, sigmaKa / dEdxEl, fGenRespKaDeltaPrimeEl, numGenEntries);
2396 errCode = GenerateDetectorResponse(errCode, 1., sigmaKa / dEdxKa, fGenRespKaDeltaPrimeKa, numGenEntries);
2397 errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxPi, sigmaKa / dEdxPi, fGenRespKaDeltaPrimePi, numGenEntries);
2398 errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxPr, sigmaKa / dEdxPr, fGenRespKaDeltaPrimePr, numGenEntries);
2399
2400 // Pions
2401 errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxEl, sigmaPi / dEdxEl, fGenRespPiDeltaPrimeEl, numGenEntries);
2402 errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxKa, sigmaPi / dEdxKa, fGenRespPiDeltaPrimeKa, numGenEntries);
2403 errCode = GenerateDetectorResponse(errCode, 1., sigmaPi / dEdxPi, fGenRespPiDeltaPrimePi, numGenEntries);
2404 errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxPr, sigmaPi / dEdxPr, fGenRespPiDeltaPrimePr, numGenEntries);
2405
2406 // Muons, if desired
2407 if (fTakeIntoAccountMuons) {
2408 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxEl, sigmaMu / dEdxEl, fGenRespMuDeltaPrimeEl, numGenEntries);
2409 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxKa, sigmaMu / dEdxKa, fGenRespMuDeltaPrimeKa, numGenEntries);
2410 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxPi, sigmaMu / dEdxPi, fGenRespMuDeltaPrimePi, numGenEntries);
2411 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxPr, sigmaMu / dEdxPr, fGenRespMuDeltaPrimePr, numGenEntries);
2412 }
2413
2414 // Protons
2415 errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxEl, sigmaPr / dEdxEl, fGenRespPrDeltaPrimeEl, numGenEntries);
2416 errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxKa, sigmaPr / dEdxKa, fGenRespPrDeltaPrimeKa, numGenEntries);
2417 errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxPi, sigmaPr / dEdxPi, fGenRespPrDeltaPrimePi, numGenEntries);
2418 errCode = GenerateDetectorResponse(errCode, 1., sigmaPr / dEdxPr, fGenRespPrDeltaPrimePr, numGenEntries);
2419
2420
2421 /*OLD with deltaSpecies
2422 // Delta
2423
2424 // Electrons
2425 errCode = GenerateDetectorResponse(errCode, 0., sigmaEl, fGenRespElDeltaEl, numGenEntries, usePureGausForDelta);
2426 errCode = GenerateDetectorResponse(errCode, dEdxEl - dEdxKa, sigmaEl, fGenRespElDeltaKa, numGenEntries, usePureGausForDelta);
2427 errCode = GenerateDetectorResponse(errCode, dEdxEl - dEdxPi, sigmaEl, fGenRespElDeltaPi, numGenEntries, usePureGausForDelta);
2428 errCode = GenerateDetectorResponse(errCode, dEdxEl - dEdxPr, sigmaEl, fGenRespElDeltaPr, numGenEntries, usePureGausForDelta);
2429
2430 // Kaons
2431 errCode = GenerateDetectorResponse(errCode, dEdxKa - dEdxEl, sigmaKa, fGenRespKaDeltaEl, numGenEntries, usePureGausForDelta);
2432 errCode = GenerateDetectorResponse(errCode, 0., sigmaKa, fGenRespKaDeltaKa, numGenEntries, usePureGausForDelta);
2433 errCode = GenerateDetectorResponse(errCode, dEdxKa - dEdxPi, sigmaKa, fGenRespKaDeltaPi, numGenEntries, usePureGausForDelta);
2434 errCode = GenerateDetectorResponse(errCode, dEdxKa - dEdxPr, sigmaKa, fGenRespKaDeltaPr, numGenEntries, usePureGausForDelta);
2435
2436 // Pions
2437 errCode = GenerateDetectorResponse(errCode, dEdxPi - dEdxEl, sigmaPi, fGenRespPiDeltaEl, numGenEntries, usePureGausForDelta);
2438 errCode = GenerateDetectorResponse(errCode, dEdxPi - dEdxKa, sigmaPi, fGenRespPiDeltaKa, numGenEntries, usePureGausForDelta);
2439 errCode = GenerateDetectorResponse(errCode, 0., sigmaPi, fGenRespPiDeltaPi, numGenEntries, usePureGausForDelta);
2440 errCode = GenerateDetectorResponse(errCode, dEdxPi - dEdxPr, sigmaPi, fGenRespPiDeltaPr, numGenEntries, usePureGausForDelta);
2441
2442 // Muons
2443 errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxEl, sigmaMu, fGenRespMuDeltaEl, numGenEntries, usePureGausForDelta);
2444 errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxKa, sigmaMu, fGenRespMuDeltaKa, numGenEntries, usePureGausForDelta);
2445 errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxPi, sigmaMu, fGenRespMuDeltaPi, numGenEntries, usePureGausForDelta);
2446 errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxPr, sigmaMu, fGenRespMuDeltaPr, numGenEntries, usePureGausForDelta);
2447
2448 // Protons
2449 errCode = GenerateDetectorResponse(errCode, dEdxPr - dEdxEl, sigmaPr, fGenRespPrDeltaEl, numGenEntries, usePureGausForDelta);
2450 errCode = GenerateDetectorResponse(errCode, dEdxPr - dEdxKa, sigmaPr, fGenRespPrDeltaKa, numGenEntries, usePureGausForDelta);
2451 errCode = GenerateDetectorResponse(errCode, dEdxPr - dEdxPi, sigmaPr, fGenRespPrDeltaPi, numGenEntries, usePureGausForDelta);
2452 errCode = GenerateDetectorResponse(errCode, 0., sigmaPr, fGenRespPrDeltaPr, numGenEntries, usePureGausForDelta);
2453 */
2454
2455 if (errCode != kNoErrors) {
2456 if (errCode == kWarning) {
2457 //Printf("Warning: Questionable detector response for track, most likely due to very low number of PID clusters! Debug output (dEdx_expected, sigma_expected):");
2458 }
2459 else
2460 Printf("Error: Failed to generate detector response for track - skipped! Debug output (dEdx_expected, sigma_expected):");
2461
2462 /*
2463 Printf("Pr: %e, %e", dEdxPr, sigmaPr);
2464 Printf("Pi: %e, %e", dEdxPi, sigmaPi);
2465 Printf("El: %e, %e", dEdxEl, sigmaEl);
2466 Printf("Mu: %e, %e", dEdxMu, sigmaMu);
2467 Printf("Ka: %e, %e", dEdxKa, sigmaKa);
2468 Printf("track: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), track->GetTPCmomentum(), track->Eta(),
2469 track->GetTPCsignalN());
2470 */
2471
2472 if (errCode != kWarning) {
2473 fhSkippedTracksForSignalGeneration->Fill(track->GetTPCmomentum(), track->GetTPCsignalN());
2474 return kFALSE;// Skip generated response in case of error
2475 }
2476 }
2477
2478 for (Int_t n = 0; n < numGenEntries; n++) {
2479 if (!isMC || !fUseMCidForGeneration) {
2480 // If no MC info is available or shall not be used, use weighting with priors to generate entries for the different species
2481 Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
2482
2483 // Consider generated response as originating from...
2484 if (rnd <= prob[AliPID::kElectron])
2485 genEntry[kGenMCID] = 0; // ... an electron
2486 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon])
2487 genEntry[kGenMCID] = 1; // ... a kaon
2488 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon])
2489 genEntry[kGenMCID] = 2; // ... a muon -> NOTE: prob[AliPID::kMuon] = 0 in case of fTakeIntoAccountMuons = kFALSE
2490 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon] + prob[AliPID::kPion])
2491 genEntry[kGenMCID] = 3; // ... a pion
2492 else
2493 genEntry[kGenMCID] = 4; // ... a proton
2494 }
2495
2496 // Electrons
2497 genEntry[kGenSelectSpecies] = 0;
2498 //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaEl[n];
2499 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeEl[n];
2500 fhGenEl->Fill(genEntry);
2501
2502 genEntry[kGenSelectSpecies] = 1;
2503 //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaKa[n];
2504 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeKa[n];
2505 fhGenEl->Fill(genEntry);
2506
2507 genEntry[kGenSelectSpecies] = 2;
2508 //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaPi[n];
2509 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePi[n];
2510 fhGenEl->Fill(genEntry);
2511
2512 genEntry[kGenSelectSpecies] = 3;
2513 //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaPr[n];
2514 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePr[n];
2515 fhGenEl->Fill(genEntry);
2516
2517 // Kaons
2518 genEntry[kGenSelectSpecies] = 0;
2519 //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaEl[n];
2520 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeEl[n];
2521 fhGenKa->Fill(genEntry);
2522
2523 genEntry[kGenSelectSpecies] = 1;
2524 //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaKa[n];
2525 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeKa[n];
2526 fhGenKa->Fill(genEntry);
2527
2528 genEntry[kGenSelectSpecies] = 2;
2529 //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaPi[n];
2530 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePi[n];
2531 fhGenKa->Fill(genEntry);
2532
2533 genEntry[kGenSelectSpecies] = 3;
2534 //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaPr[n];
2535 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePr[n];
2536 fhGenKa->Fill(genEntry);
2537
2538 // Pions
2539 genEntry[kGenSelectSpecies] = 0;
2540 //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaEl[n];
2541 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeEl[n];
2542 fhGenPi->Fill(genEntry);
2543
2544 genEntry[kGenSelectSpecies] = 1;
2545 //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaKa[n];
2546 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeKa[n];
2547 fhGenPi->Fill(genEntry);
2548
2549 genEntry[kGenSelectSpecies] = 2;
2550 //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaPi[n];
2551 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePi[n];
2552 fhGenPi->Fill(genEntry);
2553
2554 genEntry[kGenSelectSpecies] = 3;
2555 //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaPr[n];
2556 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePr[n];
2557 fhGenPi->Fill(genEntry);
2558
2559 if (fTakeIntoAccountMuons) {
2560 // Muons, if desired
2561 genEntry[kGenSelectSpecies] = 0;
2562 //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaEl[n];
2563 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeEl[n];
2564 fhGenMu->Fill(genEntry);
2565
2566 genEntry[kGenSelectSpecies] = 1;
2567 //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaKa[n];
2568 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeKa[n];
2569 fhGenMu->Fill(genEntry);
2570
2571 genEntry[kGenSelectSpecies] = 2;
2572 //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaPi[n];
2573 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePi[n];
2574 fhGenMu->Fill(genEntry);
2575
2576 genEntry[kGenSelectSpecies] = 3;
2577 //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaPr[n];
2578 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePr[n];
2579 fhGenMu->Fill(genEntry);
2580 }
2581
2582 // Protons
2583 genEntry[kGenSelectSpecies] = 0;
2584 //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaEl[n];
2585 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeEl[n];
2586 fhGenPr->Fill(genEntry);
2587
2588 genEntry[kGenSelectSpecies] = 1;
2589 //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaKa[n];
2590 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeKa[n];
2591 fhGenPr->Fill(genEntry);
2592
2593 genEntry[kGenSelectSpecies] = 2;
2594 //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaPi[n];
2595 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePi[n];
2596 fhGenPr->Fill(genEntry);
2597
2598 genEntry[kGenSelectSpecies] = 3;
2599 //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaPr[n];
2600 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePr[n];
2601 fhGenPr->Fill(genEntry);
2602 }
2603
9e95a906 2604 if(fDebug > 2)
2605 printf("File: %s, Line: %d: ProcessTrack -> Done\n", (char*)__FILE__, __LINE__);
2606
e131b05f 2607 return kTRUE;
2608}
2609
2610
2611//_____________________________________________________________________________
2612Bool_t AliAnalysisTaskPID::SetConvolutedGaussLambdaParameter(Double_t lambda)
2613{
2614 // Set the lambda parameter of the convolution to the desired value. Automatically
2615 // calculates the parameters for the transition (restricted) gauss -> convoluted gauss.
2616 fConvolutedGaussTransitionPars[2] = lambda;
2617
2618 // Save old parameters and settings of function to restore them later:
2619 Double_t* oldFuncParams = new Double_t[fkConvolutedGausNPar];
2620 for (Int_t i = 0; i < fkConvolutedGausNPar; i++)
2621 oldFuncParams[i] = fConvolutedGausDeltaPrime->GetParameter(i);
2622 Int_t oldFuncNpx = fConvolutedGausDeltaPrime->GetNpx();
2623 Double_t oldFuncRangeLow = 0, oldFuncRangeUp = 100;
2624 fConvolutedGausDeltaPrime->GetRange(oldFuncRangeLow, oldFuncRangeUp);
2625
2626 // Choose some sufficiently large range
2627 const Double_t rangeStart = 0.5;
2628 const Double_t rangeEnd = 2.0;
2629
2630 // To get the parameters for the transition, just choose arbitrarily input parameters for mu and sigma
2631 // (it makes sense to choose typical values). The ratio sigma_gauss / sigma_convolution is independent
2632 // of mu and as well the difference mu_gauss - mu_convolution.
2633 Double_t muInput = 1.0;
2634 Double_t sigmaInput = fgkSigmaReferenceForTransitionPars;
2635
2636
2637 // Step 1: Generate distribution with input parameters
2638 const Int_t nPar = 3;
2639 Double_t inputPar[nPar] = { muInput, sigmaInput, lambda };
2640
2641 TH1D* hInput = new TH1D("hInput", "Input distribution", 2000, rangeStart, rangeEnd);
2642
2643 fConvolutedGausDeltaPrime->SetParameters(inputPar);
2644 fConvolutedGausDeltaPrime->SetRange(rangeStart, rangeEnd);
2645 fConvolutedGausDeltaPrime->SetNpx(2000);
2646
2647 /*OLD
2648 // The resolution and mean of the AliPIDResponse are determined in finite intervals
2649 // of ncl (also second order effects due to finite dEdx and tanTheta).
2650 // Take this into account for the transition parameters via assuming an approximately flat
2651 // distritubtion in ncl in this interval.
2652 // NOTE: The ncl interval should be the same as the one used for the sigma map creation!
2653 const Int_t nclStart = 151;
2654 const Int_t nclEnd = 160;
2655 const Int_t nclSteps = (nclEnd - nclStart) + 1;
2656 for (Int_t ncl = nclStart; ncl <= nclEnd; ncl++) {
2657 // Resolution scales with 1/sqrt(ncl)
2658 fConvolutedGausDeltaPrime->SetParameter(1, inputPar[1] * sqrt(nclEnd) / sqrt(ncl));
2659 TH1* hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
2660
2661 for (Int_t i = 0; i < 50000000 / nclSteps; i++)
2662 hInput->Fill(hProbDensity->GetRandom());
2663 }
2664 */
2665
2666 TH1* hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
2667
2668 for (Int_t i = 0; i < 50000000; i++)
2669 hInput->Fill(hProbDensity->GetRandom());
2670
2671 // Step 2: Fit generated distribution with restricted gaussian
2672 Int_t maxBin = hInput->GetMaximumBin();
2673 Double_t max = hInput->GetBinContent(maxBin);
2674
2675 UChar_t usedBins = 1;
2676 if (maxBin > 1) {
2677 max += hInput->GetBinContent(maxBin - 1);
2678 usedBins++;
2679 }
2680 if (maxBin < hInput->GetNbinsX()) {
2681 max += hInput->GetBinContent(maxBin + 1);
2682 usedBins++;
2683 }
2684 max /= usedBins;
2685
2686 // NOTE: The range (<-> fraction of maximum) should be the same
2687 // as for the spline and eta maps creation
2688 const Double_t lowThreshold = hInput->GetXaxis()->GetBinLowEdge(hInput->FindFirstBinAbove(0.1 * max));
2689 const Double_t highThreshold = hInput->GetXaxis()->GetBinUpEdge(hInput->FindLastBinAbove(0.1 * max));
2690
2691 TFitResultPtr fitResGaussFirstStep = hInput->Fit("gaus", "QNRS", "", lowThreshold, highThreshold);
2692
2693 TFitResultPtr fitResGauss;
2694
2695 if ((Int_t)fitResGaussFirstStep == 0) {
2696 TF1 fGauss("fGauss", "[0]*TMath::Gaus(x, [1], [2], 1)", rangeStart, rangeEnd);
2697 fGauss.SetParameter(0, fitResGaussFirstStep->GetParams()[0]);
2698 fGauss.SetParError(0, fitResGaussFirstStep->GetErrors()[0]);
2699 fGauss.SetParameter(2, fitResGaussFirstStep->GetParams()[2]);
2700 fGauss.SetParError(2, fitResGaussFirstStep->GetErrors()[2]);
2701
2702 fGauss.FixParameter(1, fitResGaussFirstStep->GetParams()[1]);
2703 fitResGauss = hInput->Fit(&fGauss, "QNS", "s", rangeStart, rangeEnd);
2704 }
2705 else {
2706 fitResGauss = hInput->Fit("gaus", "QNRS", "same", rangeStart, rangeEnd);
2707 }
2708 //OLD TFitResultPtr fitResGauss = hInput->Fit("gaus", "QNRS", "", hInput->GetXaxis()->GetBinLowEdge(hInput->FindFirstBinAbove(0.1 * max)),
2709 // hInput->GetXaxis()->GetBinUpEdge(hInput->FindLastBinAbove(0.1 * max)));
2710
2711
2712 // Step 3: Use parameters from gaussian fit to obtain parameters for the transition "restricted gauss" -> "convoluted gauss"
2713
2714 // 3.1 The ratio sigmaInput / sigma_gaussFit ONLY depends on lambda (which is fixed per period) -> Calculate this first
2715 // for an arbitrary (here: typical) sigma. The ratio is then ~the same for ALL sigma for given lambda!
2716 if ((Int_t)fitResGauss != 0) {
2717 AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Gauss Fit failed!\n");
2718 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2719 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2720 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2721
2722 delete hInput;
c4856fb1 2723 delete [] oldFuncParams;
e131b05f 2724
2725 return kFALSE;
2726 }
2727
2728 if (fitResGauss->GetParams()[2] <= 0) {
2729 AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Sigma of gauss fit <= 0!\n");
2730 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2731 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2732 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2733
2734 delete hInput;
c4856fb1 2735 delete [] oldFuncParams;
e131b05f 2736
2737 return kFALSE;
2738 }
2739
2740 // sigma correction factor
2741 fConvolutedGaussTransitionPars[1] = sigmaInput / fitResGauss->GetParams()[2];
2742
2743 // 3.2 Now that sigma und lambda are determined, one can calculate mu by shifting the maximum to the desired position,
2744 // i.e. the maximum of the input distribution should coincide with that of the re-calculated distribution,
2745 // which is achieved by getting the same mu for the same sigma.
2746 // NOTE: For fixed lambda, the shift is proportional to sigma and independent of mu!
2747 // So, one can calculate the shift for an arbitrary fixed (here: typical)
2748 // sigma and then simply use this shift for any other sigma by scaling it correspondingly!!!
2749
2750 // Mu shift correction:
2751 // Shift in mu (difference between mean of gauss and mean of convolution) is proportional to sigma!
2752 // Thus, choose a reference sigma (typical -> 0.05), such that for arbitrary sigma one can simple scale
2753 // this shift correction with sigma / referenceSigma.
2754 fConvolutedGaussTransitionPars[0] = (fitResGauss->GetParams()[1] - muInput);
2755
2756
2757 /*Changed on 03.07.2013
2758
2759 // Maximum of fConvolutedGausDeltaPrime should agree with maximum of input
2760 Double_t par[nPar] = { fitResGauss->GetParams()[1], // just as a guess of the maximum position
2761 sigmaInput,
2762 lambda };
2763
2764 fConvolutedGausDeltaPrime->SetParameters(par);
2765
2766 Double_t maxXInput = fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001, muInput - 3. * sigmaInput),
2767 muInput + 10. * sigmaInput,
2768 0.0001);
2769
2770 // Maximum shifts proportional to sigma and is linear in mu (~mean of gauss)
2771 // Maximum should be typically located within [gaussMean, gaussMean + 3 gaussSigma].
2772 // -> Larger search range for safety reasons (also: sigma and/or mean might be not 100% accurate).
2773 Double_t maxXconvoluted = fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001,
2774 fitResGauss->GetParams()[1] - 3. * fitResGauss->GetParams()[2]),
2775 fitResGauss->GetParams()[1] + 10. * fitResGauss->GetParams()[2],
2776 0.0001);
2777 if (maxXconvoluted <= 0) {
2778 AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Maximum of fConvolutedGausDeltaPrime <= 0!\n");
2779
2780 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2781 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2782 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2783
2784 delete hInput;
c4856fb1 2785 delete [] oldFuncParams;
e131b05f 2786
2787 return kFALSE;
2788 }
2789
2790 // maxX perfectly shifts as par[0] (scaled by sigma) -> Can shift maxX to input value.
2791 // Mu shift correction:
2792 fConvolutedGaussTransitionPars[0] = maxXconvoluted - maxXInput;
2793 */
2794
2795
2796
2797 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2798 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2799 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2800
2801 delete hInput;
c4856fb1 2802 delete [] oldFuncParams;
e131b05f 2803
2804 return kTRUE;
2805}
2806
2807
2808//_____________________________________________________________________________
9d7ad2e4 2809AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::SetParamsForConvolutedGaus(Double_t gausMean, Double_t gausSigma)
e131b05f 2810{
2811 // Set parameters for convoluted gauss using parameters for a pure gaussian.
2812 // If SetConvolutedGaussLambdaParameter has not been called before to initialise the translation parameters,
2813 // some default parameters will be used and an error will show up.
2814
9e95a906 2815 if(fDebug > 1)
2816 printf("File: %s, Line: %d: SetParamsForConvolutedGaus: mean %e, sigma %e\n", (char*)__FILE__, __LINE__, gausMean, gausSigma);
2817
e131b05f 2818 if (fConvolutedGaussTransitionPars[1] < -998) {
2819 AliError("Transition parameters not initialised! Default parameters will be used. Please call SetConvolutedGaussLambdaParameter(...) before any calculations!");
2820 SetConvolutedGaussLambdaParameter(2.0);
2821 AliError(Form("Parameters set to:\n[0]: %f\n[1]: %f\n[2]: %f\n", fConvolutedGaussTransitionPars[0],
2822 fConvolutedGaussTransitionPars[1], fConvolutedGaussTransitionPars[2]));
2823 }
2824
2825 Double_t par[fkConvolutedGausNPar];
2826 par[2] = fConvolutedGaussTransitionPars[2];
2827 par[1] = fConvolutedGaussTransitionPars[1] * gausSigma;
2828 // maxX perfectly shifts as par[0] (scaled by sigma) -> Can shift maxX so that it sits at the right place.
2829 par[0] = gausMean - fConvolutedGaussTransitionPars[0] * par[1] / fgkSigmaReferenceForTransitionPars;
2830
2831 ErrorCode errCode = kNoErrors;
2832 fConvolutedGausDeltaPrime->SetParameters(par);
2833
9e95a906 2834 if(fDebug > 3)
2835 printf("File: %s, Line: %d: SetParamsForConvolutedGaus -> Parameters set to: %e, %e, %e (transition pars: %e, %e, %e, %e)\n",
2836 (char*)__FILE__, __LINE__, par[0], par[1], par[2], fConvolutedGaussTransitionPars[0], fConvolutedGaussTransitionPars[1],
2837 fConvolutedGaussTransitionPars[2], fgkSigmaReferenceForTransitionPars);
2838
e131b05f 2839 fConvolutedGausDeltaPrime->SetNpx(20); // Small value speeds up following algorithm (valid, since extrema far apart)
2840
2841 // Accuracy of 10^-5 is enough to get 0.1% precise peak for MIPS w.r.t. to dEdx = 2000 of protons
2842 // (should boost up the algorithm, because 10^-10 is the default value!)
2843 Double_t maxX= fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001, gausMean - 2. * gausSigma),
2844 gausMean + 6. * gausSigma, 1.0E-5);
2845
2846 const Double_t maximum = fConvolutedGausDeltaPrime->Eval(maxX);
2847 const Double_t maximumFraction = maximum * fAccuracyNonGaussianTail;
2848
2849 // Estimate lower boundary for subsequent search:
2850 Double_t lowBoundSearchBoundLow = TMath::Max(1e-4, maxX - 5. * gausSigma);
2851 Double_t lowBoundSearchBoundUp = maxX;
2852
2853 Bool_t lowerBoundaryFixedAtZero = kFALSE;
2854
2855 while (fConvolutedGausDeltaPrime->Eval(lowBoundSearchBoundLow) >= maximumFraction) {
2856 if (lowBoundSearchBoundLow <= 0) {
2857 // This should only happen to low dEdx particles with very few clusters and therefore large sigma, such that the gauss goes below zero deltaPrime
2858 if (maximum <= 0) { // Something is weired
2859 printf("Error generating signal: maximum is <= 0!\n");
2860 return kError;
2861 }
2862 else {
2863 const Double_t valueAtZero = fConvolutedGausDeltaPrime->Eval(0);
2864 if (valueAtZero / maximum > 0.05) {
2865 // Too large fraction below zero deltaPrime. Signal generation cannot be reliable in this case
2866 printf("Error generating signal: Too large fraction below zero deltaPrime: convGauss(0) / convGauss(max) = %e / %e = %e!\n",
2867 valueAtZero, maximum, valueAtZero / maximum);
2868 return kError;
2869 }
2870 }
2871
2872 /*
2873 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",
2874 fConvolutedGausDeltaPrime->Eval(maxX), fAccuracyNonGaussianTail, maximumFraction, maxX, par[0], par[1], par[2], gausMean, gausSigma);
2875 */
2876
2877 lowerBoundaryFixedAtZero = kTRUE;
2878
2879 if (errCode != kError)
2880 errCode = kWarning;
2881
2882 break;
2883 }
2884
2885 lowBoundSearchBoundUp -= gausSigma;
2886 lowBoundSearchBoundLow -= gausSigma;
2887
2888 if (lowBoundSearchBoundLow < 0) {
2889 lowBoundSearchBoundLow = 0;
2890 lowBoundSearchBoundUp += gausSigma;
2891 }
2892 }
2893
2894 // Determine lower boundary inside estimated range. For small values of the maximum: Need more precision, since finer binning!
2895 Double_t rangeStart = lowerBoundaryFixedAtZero ? 0 :
2896 fConvolutedGausDeltaPrime->GetX(maximumFraction, lowBoundSearchBoundLow, lowBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
2897
2898 // .. and the same for the upper boundary
2899 Double_t rangeEnd = 0;
2900 // If distribution starts beyond upper boundary, everything ends up in the overflow bin. So, just reduce range and Npx to minimum
2901 if (rangeStart > fkDeltaPrimeUpLimit) {
2902 rangeEnd = rangeStart + 0.00001;
2903 fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
2904 fConvolutedGausDeltaPrime->SetNpx(4);
2905 }
2906 else {
2907 // Estimate upper boundary for subsequent search:
2908 Double_t upBoundSearchBoundUp = maxX + 5 * gausSigma;
2909 Double_t upBoundSearchBoundLow = maxX;
2910 while (fConvolutedGausDeltaPrime->Eval(upBoundSearchBoundUp) >= maximumFraction) {
2911 upBoundSearchBoundUp += gausSigma;
2912 upBoundSearchBoundLow += gausSigma;
2913 }
2914
2915 // For small values of the maximum: Need more precision, since finer binning!
2916 rangeEnd = fConvolutedGausDeltaPrime->GetX(maximumFraction, upBoundSearchBoundLow, upBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
2917
2918 fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
2919 fConvolutedGausDeltaPrime->SetNpx(fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindBin(rangeEnd)
2920 - fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindBin(rangeStart) + 1);
2921 //fConvolutedGausDeltaPrime->SetNpx((rangeEnd - rangeStart) / fDeltaPrimeBinWidth + 1);
2922 }
2923
9e95a906 2924 if(fDebug > 3)
2925 printf("File: %s, Line: %d: SetParamsForConvolutedGaus -> range %f - %f, error code %d\n", (char*)__FILE__, __LINE__,
2926 rangeStart, rangeEnd, errCode);
2927
e131b05f 2928 return errCode;
2929}
2930
2931
2932//________________________________________________________________________
2933void AliAnalysisTaskPID::SetUpGenHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const
2934{
2935 // Sets bin limits for axes which are not standard binned and the axes titles.
2936
2937 hist->SetBinEdges(kGenPt, binsPt);
2938 hist->SetBinEdges(kGenDeltaPrimeSpecies, binsDeltaPrime);
2939 hist->SetBinEdges(kGenCentrality, binsCent);
2940
2941 if (fStoreAdditionalJetInformation)
2942 hist->SetBinEdges(kGenJetPt, binsJetPt);
2943
2944 // Set axes titles
2945 hist->GetAxis(kGenMCID)->SetTitle("MC PID");
2946 hist->GetAxis(kGenMCID)->SetBinLabel(1, "e");
2947 hist->GetAxis(kGenMCID)->SetBinLabel(2, "K");
2948 hist->GetAxis(kGenMCID)->SetBinLabel(3, "#mu");
2949 hist->GetAxis(kGenMCID)->SetBinLabel(4, "#pi");
2950 hist->GetAxis(kGenMCID)->SetBinLabel(5, "p");
2951
2952 hist->GetAxis(kGenSelectSpecies)->SetTitle("Select Species");
2953 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(1, "e");
2954 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(2, "K");
2955 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(3, "#pi");
2956 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(4, "p");
2957
2958 hist->GetAxis(kGenPt)->SetTitle("P_{T} (GeV/c)");
2959
2960 hist->GetAxis(kGenDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
2961
2962 hist->GetAxis(kGenCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
2963
2964 if (fStoreAdditionalJetInformation) {
2965 hist->GetAxis(kGenJetPt)->SetTitle("P_{T}^{jet} (GeV/c)");
2966
2967 hist->GetAxis(kGenZ)->SetTitle("z = P_{T}^{track} / P_{T}^{jet}");
2968
2969 hist->GetAxis(kGenXi)->SetTitle("#xi = ln(P_{T}^{jet} / P_{T}^{track})");
2970 }
2971
2972 hist->GetAxis(GetIndexOfChargeAxisGen())->SetTitle("Charge (e_{0})");
2973}
2974
2975
2976//________________________________________________________________________
2977void AliAnalysisTaskPID::SetUpGenYieldHist(THnSparse* hist, Double_t* binsPt, Double_t* binsCent, Double_t* binsJetPt) const
2978{
2979 // Sets bin limits for axes which are not standard binned and the axes titles.
2980
2981 hist->SetBinEdges(kGenYieldPt, binsPt);
2982 hist->SetBinEdges(kGenYieldCentrality, binsCent);
2983 if (fStoreAdditionalJetInformation)
2984 hist->SetBinEdges(kGenYieldJetPt, binsJetPt);
2985
2986 for (Int_t i = 0; i < 5; i++)
2987 hist->GetAxis(kGenYieldMCID)->SetBinLabel(i + 1, AliPID::ParticleLatexName(i));
2988
2989 // Set axes titles
2990 hist->GetAxis(kGenYieldMCID)->SetTitle("MC PID");
2991 hist->GetAxis(kGenYieldPt)->SetTitle("P_{T}^{gen} (GeV/c)");
2992 hist->GetAxis(kGenYieldCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
2993
2994 if (fStoreAdditionalJetInformation) {
2995 hist->GetAxis(kGenYieldJetPt)->SetTitle("P_{T}^{jet, gen} (GeV/c)");
2996
2997 hist->GetAxis(kGenYieldZ)->SetTitle("z = P_{T}^{track} / P_{T}^{jet}");
2998
2999 hist->GetAxis(kGenYieldXi)->SetTitle("#xi = ln(P_{T}^{jet} / P_{T}^{track})");
3000 }
3001
3002 hist->GetAxis(GetIndexOfChargeAxisGenYield())->SetTitle("Charge (e_{0})");
3003}
3004
3005
3006//________________________________________________________________________
3007void AliAnalysisTaskPID::SetUpHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const
3008{
3009 // Sets bin limits for axes which are not standard binned and the axes titles.
3010
3011 hist->SetBinEdges(kDataPt, binsPt);
3012 hist->SetBinEdges(kDataDeltaPrimeSpecies, binsDeltaPrime);
3013 hist->SetBinEdges(kDataCentrality, binsCent);
3014
3015 if (fStoreAdditionalJetInformation)
3016 hist->SetBinEdges(kDataJetPt, binsJetPt);
3017
3018 // Set axes titles
3019 hist->GetAxis(kDataMCID)->SetTitle("MC PID");
3020 hist->GetAxis(kDataMCID)->SetBinLabel(1, "e");
3021 hist->GetAxis(kDataMCID)->SetBinLabel(2, "K");
3022 hist->GetAxis(kDataMCID)->SetBinLabel(3, "#mu");
3023 hist->GetAxis(kDataMCID)->SetBinLabel(4, "#pi");
3024 hist->GetAxis(kDataMCID)->SetBinLabel(5, "p");
3025
3026 hist->GetAxis(kDataSelectSpecies)->SetTitle("Select Species");
3027 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(1, "e");
3028 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(2, "K");
3029 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(3, "#pi");
3030 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(4, "p");
3031
3032 hist->GetAxis(kDataPt)->SetTitle("P_{T} (GeV/c)");
3033
3034 hist->GetAxis(kDataDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
3035
3036 hist->GetAxis(kDataCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
3037
3038 if (fStoreAdditionalJetInformation) {
3039 hist->GetAxis(kDataJetPt)->SetTitle("P_{T}^{jet} (GeV/c)");
3040
3041 hist->GetAxis(kDataZ)->SetTitle("z = P_{T}^{track} / P_{T}^{jet}");
3042
3043 hist->GetAxis(kDataXi)->SetTitle("#xi = ln(P_{T}^{jet} / P_{T}^{track})");
3044 }
3045
3046 hist->GetAxis(GetIndexOfChargeAxisData())->SetTitle("Charge (e_{0})");
3047
3048 /*OLD with TOF, p_TPC_Inner and p_vertex
3049 // MC PID, SelectSpecies, P(TPC_inner), pT, p(Vertex), deltaSpecies, deltaPrimeSpecies, deltaTOFspecies
3050 hist->SetBinEdges(2, binsPt);
3051 hist->SetBinEdges(3, binsPt);
3052 hist->SetBinEdges(4, binsPt);
3053
3054 // Set axes titles
3055 hist->GetAxis(0)->SetTitle("MC PID");
3056 hist->GetAxis(0)->SetBinLabel(1, "e");
3057 hist->GetAxis(0)->SetBinLabel(2, "K");
3058 hist->GetAxis(0)->SetBinLabel(3, "#mu");
3059 hist->GetAxis(0)->SetBinLabel(4, "#pi");
3060 hist->GetAxis(0)->SetBinLabel(5, "p");
3061
3062 hist->GetAxis(1)->SetTitle("Select Species");
3063 hist->GetAxis(1)->SetBinLabel(1, "e");
3064 hist->GetAxis(1)->SetBinLabel(2, "K");
3065 hist->GetAxis(1)->SetBinLabel(3, "#pi");
3066 hist->GetAxis(1)->SetBinLabel(4, "p");
3067
3068 hist->GetAxis(2)->SetTitle("P_{TPC_inner} (GeV/c)");
3069 hist->GetAxis(3)->SetTitle("P_{T} (GeV/c)");
3070 hist->GetAxis(4)->SetTitle("P_{vertex} (GeV/c)");
3071
3072 hist->GetAxis(5)->SetTitle("TPC #Delta_{species} (arb. unit)");
3073
3074 hist->GetAxis(6)->SetTitle("TPC #Delta'_{species} (arb. unit)");
3075
3076 hist->GetAxis(7)->SetTitle("#Delta TOF_{species} (ps)");
3077 */
3078}
9e95a906 3079
3080
3081//________________________________________________________________________
e4351829 3082void AliAnalysisTaskPID::SetUpPtResHist(THnSparse* hist, Double_t* binsPt, Double_t* binsJetPt, Double_t* binsCent) const
9e95a906 3083{
3084 // Sets bin limits for axes which are not standard binned and the axes titles.
3085
3086 hist->SetBinEdges(kPtResJetPt, binsJetPt);
3087 hist->SetBinEdges(kPtResGenPt, binsPt);
3088 hist->SetBinEdges(kPtResRecPt, binsPt);
e4351829 3089 hist->SetBinEdges(kPtResCentrality, binsCent);
9e95a906 3090
3091 // Set axes titles
3092 hist->GetAxis(kPtResJetPt)->SetTitle("P_{T}^{jet, rec} (GeV/c)");
3093 hist->GetAxis(kPtResGenPt)->SetTitle("P_{T}^{gen} (GeV/c)");
3094 hist->GetAxis(kPtResRecPt)->SetTitle("P_{T}^{rec} (GeV/c)");
e4351829 3095
3096 hist->GetAxis(kPtResCharge)->SetTitle("Charge (e_{0})");
3097 hist->GetAxis(kPtResCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
9e95a906 3098}