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