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