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