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