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