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