]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/UserTasks/AliAnalysisTaskPID.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
add control histograms: reduced chi2 vs centrality, p-value versus centrality, kolgor...
[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
a6852ea8 981 Double_t dEdxTPC = fPIDResponse->IsTunedOnData() ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
982 if (dEdxTPC <= 0)
983 continue;
e131b05f 984
985 if(fTrackFilter && !fTrackFilter->IsSelected(track))
986 continue;
987
493982d9 988 if (GetUseTPCCutMIGeo()) {
a6852ea8 989 if (!TPCCutMIGeo(track, fEvent))
990 continue;
991 }
493982d9
ML		 992 else if (GetUseTPCnclCut()) {
993 if (!TPCnclCut(track))
994 continue;
995 }
e131b05f 996
997 if(fUsePhiCut) {
998 if (!PhiPrimeCut(track, magField))
999 continue; // reject track
1000 }
1001
e131b05f 1002 Int_t pdg = 0; // = 0 indicates data for the moment
1003 AliMCParticle* mcTrack = 0x0;
1004 Int_t mcID = AliPID::kUnknown;
1005 Int_t label = 0;
1006
1007 if (fMC) {
1008 label = track->GetLabel();
1009
1010 //if (label < 0)
1011 // continue;
1012
1013 mcTrack = dynamic_cast<AliMCParticle*>(fMC->GetTrack(TMath::Abs(label)));
1014 if (!mcTrack) {
1015 Printf("ERROR: Could not retrieve mcTrack with label %d for track %d", label, iTracks);
1016 continue;
1017 }
1018
1019 pdg = mcTrack->PdgCode();
1020 mcID = PDGtoMCID(mcTrack->PdgCode());
1021
9e95a906 1022 if (fDoEfficiency) {
1023 // For efficiency: Reconstructed track has survived all cuts on the detector level (excluding acceptance)
1024 // and has an associated MC track which is a physical primary and was generated inside the acceptance
1025 if (fMC->IsPhysicalPrimary(TMath::Abs(label)) &&
77324970 1026 IsInAcceptedEtaRange(TMath::Abs(mcTrack->Eta()))) {
e131b05f 1027
9e95a906 1028 // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
1029 Double_t value[kEffNumAxes] = { mcID, mcTrack->Pt(), mcTrack->Eta(), mcTrack->Charge() / 3., centralityPercentile,
1030 -1, -1, -1 };
1031 fContainerEff->Fill(value, kStepRecWithGenCuts);
1032
1033 Double_t valueMeas[kEffNumAxes] = { mcID, track->Pt(), track->Eta(), track->Charge(), centralityPercentile,
1034 -1, -1, -1 };
1035 fContainerEff->Fill(valueMeas, kStepRecWithGenCutsMeasuredObs);
1036 }
e131b05f 1037 }
1038 }
1039
1040 // Only process tracks inside the desired eta window
77324970 1041 if (!IsInAcceptedEtaRange(TMath::Abs(track->Eta()))) continue;
e131b05f 1042
9e95a906 1043 if (fDoPID)
1044 ProcessTrack(track, pdg, centralityPercentile, -1); // No jet information in this case -> Set jet pT to -1
e131b05f 1045
9e95a906 1046 if (fDoPtResolution) {
1047 if (mcTrack && fMC->IsPhysicalPrimary(TMath::Abs(label))) {
e4351829 1048 // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
1049 Double_t valuePtRes[kPtResNumAxes] = { -1, mcTrack->Pt(), track->Pt(), mcTrack->Charge() / 3., centralityPercentile };
44ed76dd 1050 FillPtResolution(mcID, valuePtRes);
9e95a906 1051 }
1052 }
1053
1054 if (fDoEfficiency) {
1055 if (mcTrack) {
1056 Double_t valueRecAllCuts[kEffNumAxes] = { mcID, track->Pt(), track->Eta(), track->Charge(), centralityPercentile,
1057 -1, -1, -1 };
1058 fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObs);
1059
1060 Double_t weight = IsSecondaryWithStrangeMotherMC(fMC, TMath::Abs(label)) ?
1061 GetMCStrangenessFactorCMS(fMC, mcTrack) : 1.0;
1062 fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObsStrangenessScaled, weight);
1063
1064 // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
1065 Double_t valueGenAllCuts[kEffNumAxes] = { mcID, mcTrack->Pt(), mcTrack->Eta(), mcTrack->Charge() / 3.,
1066 centralityPercentile, -1, -1, -1 };
1067 if (fMC->IsPhysicalPrimary(TMath::Abs(label))) {
1068 fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObsPrimaries);
1069 fContainerEff->Fill(valueGenAllCuts, kStepRecWithRecCutsPrimaries);
1070 }
1071 }
e131b05f 1072 }
1073 } //track loop
1074
9e95a906 1075 if(fDebug > 2)
1076 printf("File: %s, Line: %d: UserExec -> Processing done\n", (char*)__FILE__, __LINE__);
1077
e131b05f 1078 PostOutputData();
9e95a906 1079
1080 if(fDebug > 2)
1081 printf("File: %s, Line: %d: UserExec -> Done\n", (char*)__FILE__, __LINE__);
e131b05f 1082}
1083
1084//________________________________________________________________________
1085void AliAnalysisTaskPID::Terminate(const Option_t *)
1086{
1087 // Draw result to the screen
1088 // Called once at the end of the query
1089}
1090
1091
1092//_____________________________________________________________________________
1093void AliAnalysisTaskPID::CheckDoAnyStematicStudiesOnTheExpectedSignal()
1094{
1095 // Check whether at least one scale factor indicates the ussage of systematic studies
1096 // and set internal flag accordingly.
1097
1098 fDoAnySystematicStudiesOnTheExpectedSignal = kFALSE;
1099
21b0adb1 1100
1101 if ((TMath::Abs(fSystematicScalingSplinesBelowThreshold - 1.0) > fgkEpsilon) ||
1102 (TMath::Abs(fSystematicScalingSplinesAboveThreshold - 1.0) > fgkEpsilon)) {
e131b05f 1103 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1104 return;
1105 }
1106
1107 if ((TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon) ||
1108 (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon)) {
1109 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1110 return;
1111 }
1112
5709c40a 1113 if ((TMath::Abs(fSystematicScalingEtaSigmaParaBelowThreshold - 1.0) > fgkEpsilon) ||
1114 (TMath::Abs(fSystematicScalingEtaSigmaParaAboveThreshold - 1.0) > fgkEpsilon)) {
e131b05f 1115 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1116 return;
1117 }
1118
1119 if (TMath::Abs(fSystematicScalingMultCorrection - 1.0) > fgkEpsilon) {
1120 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1121 return;
1122 }
1123}
1124
1125
1126//_____________________________________________________________________________
1127Int_t AliAnalysisTaskPID::PDGtoMCID(Int_t pdg)
1128{
1129 // Returns the corresponding AliPID index to the given pdg code.
1130 // Returns AliPID::kUnkown if pdg belongs to a not considered species.
1131
1132 Int_t absPDGcode = TMath::Abs(pdg);
1133 if (absPDGcode == 211) {//Pion
1134 return AliPID::kPion;
1135 }
1136 else if (absPDGcode == 321) {//Kaon
1137 return AliPID::kKaon;
1138 }
1139 else if (absPDGcode == 2212) {//Proton
1140 return AliPID::kProton;
1141 }
1142 else if (absPDGcode == 11) {//Electron
1143 return AliPID::kElectron;
1144 }
1145 else if (absPDGcode == 13) {//Muon
1146 return AliPID::kMuon;
1147 }
1148
1149 return AliPID::kUnknown;
1150}
1151
1152
1153//_____________________________________________________________________________
9d7ad2e4 1154void AliAnalysisTaskPID::GetJetTrackObservables(Double_t trackPt, Double_t jetPt, Double_t& z, Double_t& xi)
e131b05f 1155{
1156 // Uses trackPt and jetPt to obtain z and xi.
1157
1158 z = (jetPt > 0 && trackPt >= 0) ? (trackPt / jetPt) : -1;
1159 xi = (z > 0) ? TMath::Log(1. / z) : -1;
1160
1161 if(trackPt > (1. - 1e-06) * jetPt && trackPt < (1. + 1e-06) * jetPt) { // case z=1 : move entry to last histo bin <1
1162 z = 1. - 1e-06;
1163 xi = 1e-06;
1164 }
1165}
1166
1167
1168//_____________________________________________________________________________
1169void AliAnalysisTaskPID::CleanupParticleFractionHistos()
1170{
1171 // Delete histos with particle fractions
1172
1173 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1174 delete fFractionHists[i];
1175 fFractionHists[i] = 0x0;
1176
1177 delete fFractionSysErrorHists[i];
1178 fFractionSysErrorHists[i] = 0x0;
1179 }
1180}
1181
1182
1183//_____________________________________________________________________________
1184Double_t AliAnalysisTaskPID::ConvolutedGaus(const Double_t* xx, const Double_t* par) const
1185{
1186 // Convolutes gauss with an exponential tail which describes dEdx-response better than pure gaussian
1187
1188 const Double_t mean = par[0];
1189 const Double_t sigma = par[1];
1190 const Double_t lambda = par[2];
1191
9e95a906 1192 if(fDebug > 5)
1193 printf("File: %s, Line: %d: ConvolutedGaus: mean %e, sigma %e, lambda %e\n", (char*)__FILE__, __LINE__, mean, sigma, lambda);
1194
e131b05f 1195 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);
1196}
1197
1198
1199//_____________________________________________________________________________
9d7ad2e4 1200inline Double_t AliAnalysisTaskPID::FastGaus(Double_t x, Double_t mean, Double_t sigma) const
e131b05f 1201{
1202 // Calculate an unnormalised gaussian function with mean and sigma.
1203
1204 if (sigma < fgkEpsilon)
1205 return 1.e30;
1206
1207 const Double_t arg = (x - mean) / sigma;
1208 return exp(-0.5 * arg * arg);
1209}
1210
1211
1212//_____________________________________________________________________________
9d7ad2e4 1213inline Double_t AliAnalysisTaskPID::FastNormalisedGaus(Double_t x, Double_t mean, Double_t sigma) const
e131b05f 1214{
1215 // Calculate a normalised (divided by sqrt(2*Pi)*sigma) gaussian function with mean and sigma.
1216
1217 if (sigma < fgkEpsilon)
1218 return 1.e30;
1219
1220 const Double_t arg = (x - mean) / sigma;
1221 const Double_t res = exp(-0.5 * arg * arg);
1222 return res / (2.50662827463100024 * sigma); //sqrt(2*Pi)=2.50662827463100024
1223}
1224
1225
1226//_____________________________________________________________________________
1227Int_t AliAnalysisTaskPID::FindBinWithinRange(TAxis* axis, Double_t value) const
1228{
1229 // Find the corresponding bin of the axis. Values outside the range (also under and overflow) will be set to the first/last
1230 // available bin
1231
1232 Int_t bin = axis->FindFixBin(value);
1233
1234 if (bin <= 0)
1235 bin = 1;
1236 if (bin > axis->GetNbins())
1237 bin = axis->GetNbins();
1238
1239 return bin;
1240}
1241
1242
1243//_____________________________________________________________________________
9d7ad2e4 1244Int_t AliAnalysisTaskPID::FindFirstBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yBin,
1245 Int_t zBin) const
e131b05f 1246{
1247 // Kind of projects a TH3 to 1 bin combination in y and z
1248 // and looks for the first x bin above a threshold for this projection.
1249 // If no such bin is found, -1 is returned.
1250
1251 if (!hist)
1252 return -1;
1253
1254 Int_t nBinsX = hist->GetNbinsX();
1255 for (Int_t xBin = 1; xBin <= nBinsX; xBin++) {
1256 if (hist->GetBinContent(xBin, yBin, zBin) > threshold)
1257 return xBin;
1258 }
1259
1260 return -1;
1261}
1262
1263
1264//_____________________________________________________________________________
9d7ad2e4 1265Int_t AliAnalysisTaskPID::FindLastBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yBin,
1266 Int_t zBin) const
e131b05f 1267{
1268 // Kind of projects a TH3 to 1 bin combination in y and z
1269 // and looks for the last x bin above a threshold for this projection.
1270 // If no such bin is found, -1 is returned.
1271
1272 if (!hist)
1273 return -1;
1274
1275 Int_t nBinsX = hist->GetNbinsX();
1276 for (Int_t xBin = nBinsX; xBin >= 1; xBin--) {
1277 if (hist->GetBinContent(xBin, yBin, zBin) > threshold)
1278 return xBin;
1279 }
1280
1281 return -1;
1282}
1283
1284
1285//_____________________________________________________________________________
9d7ad2e4 1286Bool_t AliAnalysisTaskPID::GetParticleFraction(Double_t trackPt, Double_t jetPt, Double_t centralityPercentile,
1287 AliPID::EParticleType species,
e131b05f 1288 Double_t& fraction, Double_t& fractionErrorStat, Double_t& fractionErrorSys) const
1289{
1290 // Computes the particle fraction for the corresponding species for the given trackPt, jetPt and centrality.
1291 // Use jetPt = -1 for inclusive spectra and centralityPercentile = -1 for pp.
1292 // On success (return value kTRUE), fraction contains the particle fraction, fractionErrorStat(Sys) the sigma of its
1293 // statistical (systematic) error
1294
1295 fraction = -999.;
1296 fractionErrorStat = 999.;
1297 fractionErrorSys = 999.;
1298
1299 if (species > AliPID::kProton || species < AliPID::kElectron) {
1300 AliError(Form("Only fractions for species index %d to %d availabe, but not for the requested one: %d", 0, AliPID::kProton, species));
1301 return kFALSE;
1302 }
1303
1304 if (!fFractionHists[species]) {
1305 AliError(Form("Histo with particle fractions for species %d not loaded!", species));
1306
1307 return kFALSE;
1308 }
1309
1310 Int_t jetPtBin = FindBinWithinRange(fFractionHists[species]->GetYaxis(), jetPt);
1311 Int_t centBin = FindBinWithinRange(fFractionHists[species]->GetZaxis(), centralityPercentile);
1312
1313 // The following interpolation takes the bin content of the first/last available bin,
1314 // if requested point lies beyond bin center of first/last bin.
1315 // The interpolation is only done for the x-axis (track pT), i.e. jetPtBin and centBin are fix,
1316 // because the analysis will anyhow be run in bins of jetPt and centrality and
1317 // it is not desired to correlate different jetPt bins via interpolation.
1318
1319 // The same procedure is used for the error of the fraction
1320 TAxis* xAxis = fFractionHists[species]->GetXaxis();
1321
1322 // No interpolation to values beyond the centers of the first/last bins (we don't know exactly where the spectra start or stop,
1323 // thus, search for the first and last bin above 0.0 to constrain the range
1324 Int_t firstBin = TMath::Max(1, FindFirstBinAboveIn3dSubset(fFractionHists[species], 0.0, jetPtBin, centBin));
1325 Int_t lastBin = TMath::Min(fFractionHists[species]->GetNbinsX(),
1326 FindLastBinAboveIn3dSubset(fFractionHists[species], 0.0, jetPtBin, centBin));
1327
1328 if (trackPt <= xAxis->GetBinCenter(firstBin)) {
1329 fraction = fFractionHists[species]->GetBinContent(firstBin, jetPtBin, centBin);
1330 fractionErrorStat = fFractionHists[species]->GetBinError(firstBin, jetPtBin, centBin);
1331 fractionErrorSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(firstBin, jetPtBin, centBin) : 0.;
1332 }
1333 else if (trackPt >= xAxis->GetBinCenter(lastBin)) {
1334 fraction = fFractionHists[species]->GetBinContent(lastBin, jetPtBin, centBin);
1335 fractionErrorStat = fFractionHists[species]->GetBinError(lastBin, jetPtBin, centBin);
1336 fractionErrorSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(lastBin, jetPtBin, centBin) : 0.;
1337 }
1338 else {
1339 Double_t x0 = 0., x1 = 0., y0 = 0., y1 = 0.;
1340 Double_t y0errStat = 0., y1errStat = 0., y0errSys = 0., y1errSys = 0.;
1341 Int_t trackPtBin = xAxis->FindBin(trackPt);
1342
1343 // Linear interpolation between nearest neighbours in trackPt
1344 if (trackPt <= xAxis->GetBinCenter(trackPtBin)) {
1345 y0 = fFractionHists[species]->GetBinContent(trackPtBin - 1, jetPtBin, centBin);
1346 y0errStat = fFractionHists[species]->GetBinError(trackPtBin - 1, jetPtBin, centBin);
1347 y0errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin - 1, jetPtBin, centBin)
1348 : 0.;
1349 x0 = xAxis->GetBinCenter(trackPtBin - 1);
1350 y1 = fFractionHists[species]->GetBinContent(trackPtBin, jetPtBin, centBin);
1351 y1errStat = fFractionHists[species]->GetBinError(trackPtBin, jetPtBin, centBin);
1352 y1errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin, jetPtBin, centBin)
1353 : 0.;
1354 x1 = xAxis->GetBinCenter(trackPtBin);
1355 }
1356 else {
1357 y0 = fFractionHists[species]->GetBinContent(trackPtBin, jetPtBin, centBin);
1358 y0errStat = fFractionHists[species]->GetBinError(trackPtBin, jetPtBin, centBin);
1359 y0errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin, jetPtBin, centBin)
1360 : 0.;
1361 x0 = xAxis->GetBinCenter(trackPtBin);
1362 y1 = fFractionHists[species]->GetBinContent(trackPtBin + 1, jetPtBin, centBin);
1363 y1errStat = fFractionHists[species]->GetBinError(trackPtBin + 1, jetPtBin, centBin);
1364 y1errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin + 1, jetPtBin, centBin)
1365 : 0.;
1366 x1 = xAxis->GetBinCenter(trackPtBin + 1);
1367 }
1368
1369 // Per construction: x0 < trackPt < x1
1370 fraction = y0 + (trackPt - x0) * ((y1 - y0) / (x1 - x0));
1371 fractionErrorStat = y0errStat + (trackPt - x0) * ((y1errStat - y0errStat) / (x1 - x0));
1372 fractionErrorSys = fFractionSysErrorHists[species] ? (y0errSys + (trackPt - x0) * ((y1errSys - y0errSys) / (x1 - x0))) : 0.;
1373 }
1374
1375 return kTRUE;
1376}
1377
1378
1379//_____________________________________________________________________________
9d7ad2e4 1380Bool_t AliAnalysisTaskPID::GetParticleFractions(Double_t trackPt, Double_t jetPt, Double_t centralityPercentile,
1381 Double_t* prob, Int_t smearSpeciesByError,
1382 Int_t takeIntoAccountSpeciesSysError, Bool_t uniformSystematicError) const
e131b05f 1383{
1384 // Fills the particle fractions for the given trackPt, jetPt and centrality into "prob".
1385 // Use jetPt = -1 for inclusive spectra and centralityPercentile = -1 for pp.
1386 // If smearSpeciesByError is >= 0 && < AliPID::kSPECIES, the returned fractions will be a random number distributed
1387 // with a gauss with mean being the corresponding particle fraction and sigma it's error for the considered species
1388 // "smearSpeciesByError".
1389 // Note that in this case the fractions for all species will NOT sum up to 1!
1390 // Thus, all other species fractions will be re-scaled weighted with their corresponding statistical error.
1391 // A similar procedure is used for "takeIntoAccountSpeciesSysError": The systematic error of the corresponding species
1392 // is used to generate a random number with uniform distribution in [mean - sysError, mean + sysError] for the new mean
1393 // (in cace of uniformSystematicError = kTRUE, otherwise it will be a gaus(mean, sysError)),
1394 // then the other species will be re-scaled according to their systematic errors.
1395 // First, the systematic error uncertainty procedure will be performed (that is including re-scaling), then the statistical
1396 // uncertainty procedure.
1397 // On success, kTRUE is returned.
1398
1399 if (!prob || smearSpeciesByError >= AliPID::kSPECIES || takeIntoAccountSpeciesSysError >= AliPID::kSPECIES)
1400 return kFALSE;
1401
1402 Double_t probTemp[AliPID::kSPECIES];
1403 Double_t probErrorStat[AliPID::kSPECIES];
1404 Double_t probErrorSys[AliPID::kSPECIES];
1405
1406 Bool_t success = kTRUE;
1407 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kElectron,
1408 probTemp[AliPID::kElectron], probErrorStat[AliPID::kElectron],
1409 probErrorSys[AliPID::kElectron]);
1410 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kMuon,
1411 probTemp[AliPID::kMuon], probErrorStat[AliPID::kMuon], probErrorSys[AliPID::kMuon]);
1412 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kPion,
1413 probTemp[AliPID::kPion], probErrorStat[AliPID::kPion], probErrorSys[AliPID::kPion]);
1414 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kKaon,
1415 probTemp[AliPID::kKaon], probErrorStat[AliPID::kKaon], probErrorSys[AliPID::kKaon]);
1416 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kProton,
1417 probTemp[AliPID::kProton], probErrorStat[AliPID::kProton], probErrorSys[AliPID::kProton]);
1418
1419 if (!success)
1420 return kFALSE;
1421
1422 // If desired, take into account the systematic error of the corresponding species and re-generate probTemp accordingly
1423 if (takeIntoAccountSpeciesSysError >= 0) {
1424 // Generate random fraction of the considered species "smearSpeciesByError" according to mean and sigma
1425 Double_t generatedFraction = uniformSystematicError
1426 ? fRandom->Rndm() * 2. * probErrorSys[takeIntoAccountSpeciesSysError]
1427 - probErrorSys[takeIntoAccountSpeciesSysError]
1428 + probTemp[takeIntoAccountSpeciesSysError]
1429 : fRandom->Gaus(probTemp[takeIntoAccountSpeciesSysError],
1430 probErrorSys[takeIntoAccountSpeciesSysError]);
1431
1432 // Catch cases with invalid fraction (can happen for large errors), i.e. fraction < 0 or > 1
1433 if (generatedFraction < 0.)
1434 generatedFraction = 0.;
1435 else if (generatedFraction > 1.)
1436 generatedFraction = 1.;
1437
1438 // Calculate difference from original fraction (original fractions sum up to 1!)
1439 Double_t deltaFraction = generatedFraction - probTemp[takeIntoAccountSpeciesSysError];
1440
1441 // Fractions must (including errors) lie inside [0,1] -> Adapt weights accordingly by setting the errors
1442 if (deltaFraction > 0) {
1443 // Some part will be SUBTRACTED from the other fractions
1444 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1445 if (probTemp[i] - probErrorSys[i] < 0)
1446 probErrorSys[i] = probTemp[i];
1447 }
1448 }
1449 else {
1450 // Some part will be ADDED to the other fractions
1451 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1452 if (probTemp[i] + probErrorSys[i] > 1)
1453 probErrorSys[i] = 1. - probTemp[i];
1454 }
1455 }
1456
1457 // Compute summed weight of all fractions except for the considered one
1458 Double_t summedWeight = 0.;
1459 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1460 if (i != takeIntoAccountSpeciesSysError)
1461 summedWeight += probErrorSys[i];
1462 }
1463
1464 // Compute the weight for the other species
1465 /*
1466 if (summedWeight <= 1e-13) {
1467 // If this happens for some reason (it should not!), just assume flat weight
1468 printf("Error: summedWeight (sys error) ~ 0 for trackPt %f, jetPt %f, centralityPercentile %f. Setting flat weight!\n",
1469 trackPt, jetPt, centralityPercentile);
1470 }*/
1471
1472 Double_t weight[AliPID::kSPECIES];
1473 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1474 if (i != takeIntoAccountSpeciesSysError) {
1475 if (summedWeight > 1e-13)
1476 weight[i] = probErrorSys[i] / summedWeight;
1477 else
1478 weight[i] = probErrorSys[i] / (AliPID::kSPECIES - 1);
1479 }
1480 }
1481
1482 // For the final generated fractions, set the generated value for the considered species
1483 // and the generated value minus delta times statistical weight
1484 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1485 if (i != takeIntoAccountSpeciesSysError)
1486 probTemp[i] = probTemp[i] - weight[i] * deltaFraction;
1487 else
1488 probTemp[i] = generatedFraction;
1489 }
1490 }
1491
1492 // Using the values of probTemp (either the original ones or those after taking into account the systematic error),
1493 // calculate the final fractions - if the statistical error is to be taken into account, smear the corresponding
1494 // fraction. If not, just write probTemp to the final result array.
1495 if (smearSpeciesByError >= 0) {
1496 // Generate random fraction of the considered species "smearSpeciesByError" according to mean and sigma
1497 Double_t generatedFraction = fRandom->Gaus(probTemp[smearSpeciesByError], probErrorStat[smearSpeciesByError]);
1498
1499 // Catch cases with invalid fraction (can happen for large errors), i.e. fraction < 0 or > 1
1500 if (generatedFraction < 0.)
1501 generatedFraction = 0.;
1502 else if (generatedFraction > 1.)
1503 generatedFraction = 1.;
1504
1505 // Calculate difference from original fraction (original fractions sum up to 1!)
1506 Double_t deltaFraction = generatedFraction - probTemp[smearSpeciesByError];
1507
1508 // Fractions must (including errors) lie inside [0,1] -> Adapt weights accordingly by setting the errors
1509 if (deltaFraction > 0) {
1510 // Some part will be SUBTRACTED from the other fractions
1511 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1512 if (probTemp[i] - probErrorStat[i] < 0)
1513 probErrorStat[i] = probTemp[i];
1514 }
1515 }
1516 else {
1517 // Some part will be ADDED to the other fractions
1518 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1519 if (probTemp[i] + probErrorStat[i] > 1)
1520 probErrorStat[i] = 1. - probTemp[i];
1521 }
1522 }
1523
1524 // Compute summed weight of all fractions except for the considered one
1525 Double_t summedWeight = 0.;
1526 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1527 if (i != smearSpeciesByError)
1528 summedWeight += probErrorStat[i];
1529 }
1530
1531 // Compute the weight for the other species
1532 /*
1533 if (summedWeight <= 1e-13) {
1534 // If this happens for some reason (it should not!), just assume flat weight
1535 printf("Error: summedWeight (stat error) ~ 0 for trackPt %f, jetPt %f, centralityPercentile %f. Setting flat weight!\n",
1536 trackPt, jetPt, centralityPercentile);
1537 }*/
1538
1539 Double_t weight[AliPID::kSPECIES];
1540 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1541 if (i != smearSpeciesByError) {
1542 if (summedWeight > 1e-13)
1543 weight[i] = probErrorStat[i] / summedWeight;
1544 else
1545 weight[i] = probErrorStat[i] / (AliPID::kSPECIES - 1);
1546 }
1547 }
1548
1549 // For the final generated fractions, set the generated value for the considered species
1550 // and the generated value minus delta times statistical weight
1551 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1552 if (i != smearSpeciesByError)
1553 prob[i] = probTemp[i] - weight[i] * deltaFraction;
1554 else
1555 prob[i] = generatedFraction;
1556 }
1557 }
1558 else {
1559 // Just take the generated values
1560 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
1561 prob[i] = probTemp[i];
1562 }
1563
1564
1565 // Should already be normalised, but make sure that it really is:
1566 Double_t probSum = 0.;
1567 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1568 probSum += prob[i];
1569 }
1570
1571 if (probSum <= 0)
1572 return kFALSE;
1573
1574 if (TMath::Abs(probSum - 1.0) > 1e-4) {
1575 printf("Warning: Re-normalising sum of fractions: Sum is %e\n", probSum);
1576 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1577 prob[i] /= probSum;
1578 }
1579 }
1580
1581 return kTRUE;
1582}
1583
1584
1585//_____________________________________________________________________________
9d7ad2e4 1586const TH3D* AliAnalysisTaskPID::GetParticleFractionHisto(Int_t species, Bool_t sysError) const
e131b05f 1587{
1588 if (species < AliPID::kElectron || species > AliPID::kProton)
1589 return 0x0;
1590
1591 return sysError ? fFractionSysErrorHists[species] : fFractionHists[species];
1592}
1593
1594
1595//_____________________________________________________________________________
1596Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(Int_t motherPDG, Double_t motherGenPt)
1597{
1598 // Strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
1599 // -> Based on function in PWGJE/AliAnalysisTaskFragmentationFunction, which uses
1600 // the following data from CMS pp @ 7 TeV inclusive (JHEP 05 (2011) 064)
1601
1602 Double_t fac = 1;
1603
1604 const Int_t absMotherPDG = TMath::Abs(motherPDG);
1605
1606 if (absMotherPDG == 310 || absMotherPDG == 321) { // K0s / K+ / K-
1607 if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.768049;
1608 else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.732933;
1609 else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.650298;
1610 else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.571332;
1611 else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.518734;
1612 else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.492543;
1613 else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.482704;
1614 else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.488056;
1615 else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.488861;
1616 else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.492862;
1617 else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.504332;
1618 else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.501858;
1619 else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.512970;
1620 else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.524131;
1621 else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.539130;
1622 else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.554101;
1623 else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.560348;
1624 else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.568869;
1625 else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.583310;
1626 else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.604818;
1627 else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.632630;
1628 else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.710070;
1629 else if(6.00 <= motherGenPt && motherGenPt < 8.00) fac = 0.736365;
1630 else if(8.00 <= motherGenPt && motherGenPt < 10.00) fac = 0.835865;
1631 }
1632
1633 if (absMotherPDG == 3122) { // Lambda
1634 if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.645162;
1635 else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.627431;
1636 else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.457136;
1637 else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.384369;
1638 else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.330597;
1639 else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.309571;
1640 else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.293620;
1641 else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.283709;
1642 else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.282047;
1643 else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.277261;
1644 else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.275772;
1645 else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.280726;
1646 else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.288540;
1647 else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.288315;
1648 else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.296619;
1649 else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.302993;
1650 else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.338121;
1651 else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.349800;
1652 else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.356802;
1653 else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.391202;
1654 else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.422573;
1655 else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.573815;
1656 else if(6.00 <= motherGenPt && motherGenPt < 8.00) fac = 0.786984;
1657 else if(8.00 <= motherGenPt && motherGenPt < 10.00) fac = 1.020021;
1658 }
1659
1660 if (absMotherPDG == 3312 || absMotherPDG == 3322) { // xi
1661 if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.666620;
1662 else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.575908;
1663 else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.433198;
1664 else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.340901;
1665 else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.290896;
1666 else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.236074;
1667 else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.218681;
1668 else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.207763;
1669 else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.222848;
1670 else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.208806;
1671 else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.197275;
1672 else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.183645;
1673 else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.188788;
1674 else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.188282;
1675 else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.207442;
1676 else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.240388;
1677 else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.241916;
1678 else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.208276;
1679 else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.234550;
1680 else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.251689;
1681 else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.310204;
1682 else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.343492;
1683 }
1684
1685 const Double_t weight = 1. / fac;
1686
1687 return weight;
1688}
1689
1690
1691//_____________________________________________________________________________
1692Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(AliMCEvent* mcEvent, AliMCParticle* daughter)
1693{
1694 // Strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
1695 // -> Based on function in PWGJE/AliAnalysisTaskFragmentationFunction
1696
1697 if (!mcEvent)
1698 return 1.;
1699
1700 AliMCParticle* currentMother = daughter;
1701 AliMCParticle* currentDaughter = daughter;
1702
1703
1704 // find first primary mother K0s, Lambda or Xi
1705 while(1) {
1706 Int_t daughterPDG = currentDaughter->PdgCode();
1707
1708 Int_t motherLabel = currentDaughter->GetMother();
1709 if(motherLabel >= mcEvent->GetNumberOfTracks()){ // protection
1710 currentMother = currentDaughter;
1711 break;
1712 }
1713
1714 currentMother = (AliMCParticle*)mcEvent->GetTrack(motherLabel);
1715
1716 if (!currentMother) {
1717 currentMother = currentDaughter;
1718 break;
1719 }
1720
1721 Int_t motherPDG = currentMother->PdgCode();
1722
1723 // phys. primary found ?
1724 if (mcEvent->IsPhysicalPrimary(motherLabel))
1725 break;
1726
1727 if (TMath::Abs(daughterPDG) == 321) {
1728 // K+/K- e.g. from phi (ref data not feeddown corrected)
1729 currentMother = currentDaughter;
1730 break;
1731 }
1732 if (TMath::Abs(motherPDG) == 310) {
1733 // K0s e.g. from phi (ref data not feeddown corrected)
1734 break;
1735 }
1736 if (TMath::Abs(motherPDG) == 3212 && TMath::Abs(daughterPDG) == 3122) {
1737 // Mother Sigma0, daughter Lambda (this case not included in feeddown corr.)
1738 currentMother = currentDaughter;
1739 break;
1740 }
1741
1742 currentDaughter = currentMother;
1743 }
1744
1745
1746 Int_t motherPDG = currentMother->PdgCode();
1747 Double_t motherGenPt = currentMother->Pt();
1748
1749 return GetMCStrangenessFactorCMS(motherPDG, motherGenPt);
1750}
1751
1752
77324970
CKB		 1753// _________________________________________________________________________________
1754AliAnalysisTaskPID::TOFpidInfo AliAnalysisTaskPID::GetTOFType(const AliVTrack* track, Int_t tofMode) const
1755{
1756 // Get the (locally defined) particle type judged by TOF
cadb37bc 1757
77324970
CKB		 1758 if (!fPIDResponse) {
1759 Printf("ERROR: fEvent not available -> Cannot determine TOF type!");
1760 return kNoTOFinfo;
1761 }
1762
1763 // Check kTOFout, kTIME, mismatch
1764 const AliPIDResponse::EDetPidStatus tofStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, track);
1765 if (tofStatus != AliPIDResponse::kDetPidOk)
1766 return kNoTOFinfo;
cadb37bc 1767
1768 Double_t nsigma[kNumTOFspecies + 1] = { -999., -999., -999., -999. };
1769 const Int_t kTOFelectron = kTOFproton + 1;
77324970
CKB		 1770 nsigma[kTOFpion] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion);
1771 nsigma[kTOFkaon] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon);
1772 nsigma[kTOFproton] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton);
cadb37bc 1773 nsigma[kTOFelectron] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kElectron);
1774
77324970
CKB		 1775 Double_t inclusion = -999;
1776 Double_t exclusion = -999;
1777
1778 if (tofMode == 0) {
275a4969 1779 inclusion = 1.;
1780 exclusion = 2.;
77324970 1781 }
275a4969 1782 else if (tofMode == 1) { // default
1783 inclusion = 1.;
ea34c37d 1784 exclusion = 2.5;
77324970
CKB		 1785 }
1786 else if (tofMode == 2) {
275a4969 1787 inclusion = 1.5;
1788 exclusion = 2.;
77324970
CKB		 1789 }
1790 else {
1791 Printf("ERROR: Bad TOF mode: %d!", tofMode);
1792 return kNoTOFinfo;
1793 }
cadb37bc 1794
1795 // Smaller exclusion cut for electron band in order not to sacrifise too much TOF pions,
1796 // but still have a reasonably small electron contamination
ea34c37d 1797 Double_t exclusionForEl = 1.5;
cadb37bc 1798
1799 // Exclusion cut on electrons for pions because the precision of pions is good and
1800 // the contamination of electron can not be ignored (although effect on pions is small
1801 // due to overall small electron fraction, the contamination would completely bias the
1802 // electron fraction).
1803 // The electron exclsuion cut is also applied to kaons and protons for consistency, but
1804 // there should be no effect. This is because there is already the exclusion cut on pions
1805 // and pions and electrons completely overlap in the region, where electrons and pions
1806 // fall inside the inclusion cut of kaons/protons.
1807 if (TMath::Abs(nsigma[kTOFpion]) < inclusion && TMath::Abs(nsigma[kTOFkaon]) > exclusion && TMath::Abs(nsigma[kTOFproton]) > exclusion
1808 && TMath::Abs(nsigma[kTOFelectron]) > exclusionForEl)
77324970 1809 return kTOFpion;
cadb37bc 1810 if (TMath::Abs(nsigma[kTOFpion]) > exclusion && TMath::Abs(nsigma[kTOFkaon]) < inclusion && TMath::Abs(nsigma[kTOFproton]) > exclusion
1811 && TMath::Abs(nsigma[kTOFelectron]) > exclusionForEl)
77324970 1812 return kTOFkaon;
cadb37bc 1813 if (TMath::Abs(nsigma[kTOFpion]) > exclusion && TMath::Abs(nsigma[kTOFkaon]) > exclusion && TMath::Abs(nsigma[kTOFproton]) < inclusion
1814 && TMath::Abs(nsigma[kTOFelectron]) > exclusionForEl)
77324970 1815 return kTOFproton;
cadb37bc 1816
1817 // There are no TOF electrons selected because the purity is rather bad, even for small momenta
1818 // (also a small mismatch probability significantly affects electrons because their fraction
1819 // is small). There is no need for TOF electrons anyway, since the dEdx distribution of kaons and
1820 // protons in a given pT bin (also at the dEdx crossings) is very different from that of electrons.
1821 // This is due to the steeply falling dEdx of p and K with momentum, whereas the electron dEdx stays
1822 // rather constant.
1823 // As a result, the TPC fit yields a more accurate electron fraction than the TOF selection can do.
1824
77324970
CKB		 1825 return kNoTOFpid;
1826}
1827
1828
e131b05f 1829// _________________________________________________________________________________
1830Bool_t AliAnalysisTaskPID::IsSecondaryWithStrangeMotherMC(AliMCEvent* mcEvent, Int_t partLabel)
1831{
1832 // Check whether particle is a secondary with strange mother, i.e. returns kTRUE if a strange mother is found
1833 // and the particle is NOT a physical primary. In all other cases kFALSE is returned
1834
1835 if (!mcEvent || partLabel < 0)
1836 return kFALSE;
1837
1838 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(partLabel);
1839
1840 if (!part)
1841 return kFALSE;
1842
1843 if (mcEvent->IsPhysicalPrimary(partLabel))
1844 return kFALSE;
1845
1846 Int_t iMother = part->GetMother();
1847 if (iMother < 0)
1848 return kFALSE;
1849
1850
1851 AliMCParticle* partM = (AliMCParticle*)mcEvent->GetTrack(iMother);
1852 if (!partM)
1853 return kFALSE;
1854
1855 Int_t codeM = TMath::Abs(partM->PdgCode());
1856 Int_t mfl = Int_t(codeM / TMath::Power(10, Int_t(TMath::Log10(codeM))));
1857 if (mfl == 3 && codeM != 3) // codeM = 3 is for s quark
1858 return kTRUE;
1859
1860 return kFALSE;
1861}
1862
1863
1864//_____________________________________________________________________________
9d7ad2e4 1865Bool_t AliAnalysisTaskPID::SetParticleFractionHisto(const TH3D* hist, Int_t species, Bool_t sysError)
e131b05f 1866{
1867 // Store a clone of hist (containing the particle fractions of the corresponding species with statistical error (sysError = kFALSE)
1868 // or systematic error (sysError = kTRUE), respectively), internally
1869
1870 if (species < AliPID::kElectron || species > AliPID::kProton) {
1871 AliError(Form("Only fractions for species index %d to %d can be set, but not for the requested one: %d", 0,
1872 AliPID::kProton, species));
1873 return kFALSE;
1874 }
1875
1876 if (sysError) {
1877 delete fFractionSysErrorHists[species];
1878
1879 fFractionSysErrorHists[species] = new TH3D(*hist);
1880 }
1881 else {
1882 delete fFractionHists[species];
1883
1884 fFractionHists[species] = new TH3D(*hist);
1885 }
1886
1887 return kTRUE;
1888}
1889
1890
1891//_____________________________________________________________________________
9d7ad2e4 1892Bool_t AliAnalysisTaskPID::SetParticleFractionHistosFromFile(const TString filePathName, Bool_t sysError)
e131b05f 1893{
1894 // Loads particle fractions for all species from the desired file and returns kTRUE on success.
1895 // The maps are assumed to be of Type TH3D, to sit in the main directory and to have names
1896 // Form("hFraction_%e", AliPID::ParticleName(i)) for sysError = kFALSE and
1897 // Form("hFractionSysError_%e", AliPID::ParticleName(i)) for sysError = kTRUE.
1898
1899 TFile* f = TFile::Open(filePathName.Data());
1900 if (!f) {
1901 std::cout << "Failed to open file with particle fractions \"" << filePathName.Data() << "\"!" << std::endl;
1902 return kFALSE;
1903 }
1904
1905 TH3D* hist = 0x0;
1906 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1907 TString histName = Form("hFraction%s_%s", sysError ? "SysError" : "", AliPID::ParticleName(i));
1908 hist = dynamic_cast<TH3D*>(f->Get(histName.Data()));
1909 if (!hist) {
1910 std::cout << "Failed to load particle fractions for " << histName.Data() << "!";
1911 std::cout << std::endl << "Cleaning up particle fraction histos!" << std::endl;
1912 CleanupParticleFractionHistos();
1913 return kFALSE;
1914 }
1915
1916 if (!SetParticleFractionHisto(hist, i, sysError)) {
1917 std::cout << "Failed to load particle fractions for " << histName.Data() << "!";
1918 std::cout << std::endl << "Cleaning up particle fraction histos!" << std::endl;
1919 CleanupParticleFractionHistos();
1920 return kFALSE;
1921 }
1922 }
1923
1924 delete hist;
1925
1926 return kTRUE;
1927
1928}
1929
1930
1931//_____________________________________________________________________________
9d7ad2e4 1932Int_t AliAnalysisTaskPID::GetRandomParticleTypeAccordingToParticleFractions(Double_t trackPt, Double_t jetPt,
1933 Double_t centralityPercentile,
1934 Bool_t smearByError,
1935 Bool_t takeIntoAccountSysError) const
e131b05f 1936{
1937 // Uses the stored histograms with the particle fractions to generate a random particle type according to these fractions.
1938 // In case of problems (e.g. histo missing), AliPID::kUnknown is returned.
1939 // If smearByError is kTRUE, the used fractions will be random numbers distributed with a gauss with mean
1940 // being the corresponding particle fraction and sigma it's error.
1941 // Note that in this case only the fraction of a random species is varied in this way. The other fractions
1942 // will be re-normalised according their statistical errors.
1943 // The same holds for the systematic error of species "takeIntoAccountSpeciesSysError", but the random number will be
1944 // uniformly distributed within [mean - sys, mean + sys] and the re-normalisation will be weighted with the systematic errors.
1945 // Note that the fractions will be calculated first with only the systematic error taken into account (if desired), including
1946 // re-normalisation. Then, the resulting fractions will be used to calculate the final fractions - either with statistical error
1947 // or without. The species, for which the error will be used for smearing, is the same for sys and stat error.
1948
1949 Double_t prob[AliPID::kSPECIES];
1950 Int_t randomSpecies = (smearByError || takeIntoAccountSysError) ? (Int_t)(fRandom->Rndm() * AliPID::kSPECIES) : -1;
1951 Bool_t success = GetParticleFractions(trackPt, jetPt, centralityPercentile, prob, randomSpecies, randomSpecies);
1952
1953 if (!success)
1954 return AliPID::kUnknown;
1955
1956 Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
1957
1958 if (rnd <= prob[AliPID::kPion])
1959 return AliPID::kPion;
1960 else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon])
1961 return AliPID::kKaon;
1962 else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon] + prob[AliPID::kProton])
1963 return AliPID::kProton;
1964 else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon] + prob[AliPID::kProton] + prob[AliPID::kElectron])
1965 return AliPID::kElectron;
1966
1967 return AliPID::kMuon; //else it must be a muon (only species left)
1968}
1969
1970
1971//_____________________________________________________________________________
9d7ad2e4 1972AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::GenerateDetectorResponse(AliAnalysisTaskPID::ErrorCode errCode,
1973 Double_t mean, Double_t sigma,
1974 Double_t* responses, Int_t nResponses,
1975 Bool_t usePureGaus)
e131b05f 1976{
1977 // Generate detector response. If a previous generation was not successful or there is something wrong with this signal generation,
1978 // the function will return kFALSE
1979 if (!responses)
1980 return kError;
1981
1982 // Reset response array
1983 for (Int_t i = 0; i < nResponses; i++)
1984 responses[i] = -999;
1985
1986 if (errCode == kError)
1987 return kError;
1988
1989 ErrorCode ownErrCode = kNoErrors;
1990
1991 if (fUseConvolutedGaus && !usePureGaus) {
1992 // In case of convoluted gauss, calculate the probability density only once to save a lot of time!
1993
1994 TH1* hProbDensity = 0x0;
1995 ownErrCode = SetParamsForConvolutedGaus(mean, sigma);
1996 if (ownErrCode == kError)
1997 return kError;
1998
1999 hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
2000
2001 for (Int_t i = 0; i < nResponses; i++) {
2002 responses[i] = hProbDensity->GetRandom();
2003 //responses[i] fConvolutedGausDeltaPrime->GetRandom(); // MUCH slower than using the binned version via the histogram
2004 }
2005 }
2006 else {
2007 for (Int_t i = 0; i < nResponses; i++) {
2008 responses[i] = fRandom->Gaus(mean, sigma);
2009 }
2010 }
2011
2012 // If forwarded error code was a warning (error case has been handled before), return a warning
2013 if (errCode == kWarning)
2014 return kWarning;
2015
2016 return ownErrCode; // Forward success/warning
2017}
2018
2019
2020//_____________________________________________________________________________
2021void AliAnalysisTaskPID::PrintSettings(Bool_t printSystematicsSettings) const
2022{
2023 // Print current settings.
2024
2025 printf("\n\nSettings for task %s:\n", GetName());
2026 printf("Is pPb/Pbp: %d -> %s\n", GetIsPbpOrpPb(), GetIsPbpOrpPb() ? "Adapting vertex cuts" : "Using standard vertex cuts");
2027 printf("Track cuts: %s\n", fTrackFilter ? fTrackFilter->GetTitle() : "-");
2028 printf("Eta cut: %.2f <= |eta| <= %.2f\n", GetEtaAbsCutLow(), GetEtaAbsCutUp());
2029 printf("Phi' cut: %d\n", GetUsePhiCut());
a6852ea8 2030 printf("TPCCutMIGeo: %d\n", GetUseTPCCutMIGeo());
493982d9
ML		 2031 if (GetUseTPCCutMIGeo()) {
2032 printf("GetCutGeo: %f\n", GetCutGeo());
2033 printf("GetCutNcr: %f\n", GetCutNcr());
2034 printf("GetCutNcl: %f\n", GetCutNcl());
2035 }
2036 printf("TPCnclCut: %d\n", GetUseTPCnclCut());
2037 if (GetUseTPCnclCut()) {
2038 printf("GetCutPureNcl: %d\n", GetCutPureNcl());
2039 }
e131b05f 2040
2041 printf("\n");
2042
2043 printf("Centrality estimator: %s\n", GetCentralityEstimator().Data());
2044
2045 printf("\n");
2046
2047 printf("Use MC-ID for signal generation: %d\n", GetUseMCidForGeneration());
2048 printf("Use ITS: %d\n", GetUseITS());
2049 printf("Use TOF: %d\n", GetUseTOF());
2050 printf("Use priors: %d\n", GetUsePriors());
2051 printf("Use TPC default priors: %d\n", GetUseTPCDefaultPriors());
2052 printf("Use convoluted Gauss: %d\n", GetUseConvolutedGaus());
2053 printf("Accuracy of non-Gaussian tail: %e\n", GetAccuracyNonGaussianTail());
2054 printf("Take into account muons: %d\n", GetTakeIntoAccountMuons());
77324970 2055 printf("TOF mode: %d\n", GetTOFmode());
e131b05f 2056 printf("\nParams for transition from gauss to asymmetric shape:\n");
2057 printf("[0]: %e\n", GetConvolutedGaussTransitionPar(0));
2058 printf("[1]: %e\n", GetConvolutedGaussTransitionPar(1));
2059 printf("[2]: %e\n", GetConvolutedGaussTransitionPar(2));
2060
9e95a906 2061 printf("\n");
2062
2063 printf("Do PID: %d\n", fDoPID);
2064 printf("Do Efficiency: %d\n", fDoEfficiency);
2065 printf("Do PtResolution: %d\n", fDoPtResolution);
2066
e131b05f 2067 printf("\n");
2068
2069 printf("Input from other task: %d\n", GetInputFromOtherTask());
2070 printf("Store additional jet information: %d\n", GetStoreAdditionalJetInformation());
2071 printf("Store centrality percentile: %d", GetStoreCentralityPercentile());
2072
2073 if (printSystematicsSettings)
2074 PrintSystematicsSettings();
2075 else
2076 printf("\n\n\n");
2077}
2078
2079
2080//_____________________________________________________________________________
2081void AliAnalysisTaskPID::PrintSystematicsSettings() const
2082{
2083 // Print current settings for systematic studies.
2084
2085 printf("\n\nSettings for systematics for task %s:\n", GetName());
21b0adb1 2086 printf("SplinesThr:\t%f\n", GetSystematicScalingSplinesThreshold());
2087 printf("SplinesBelowThr:\t%f\n", GetSystematicScalingSplinesBelowThreshold());
2088 printf("SplinesAboveThr:\t%f\n", GetSystematicScalingSplinesAboveThreshold());
e131b05f 2089 printf("EtaCorrMomThr:\t%f\n", GetSystematicScalingEtaCorrectionMomentumThr());
2090 printf("EtaCorrLowP:\t%f\n", GetSystematicScalingEtaCorrectionLowMomenta());
2091 printf("EtaCorrHighP:\t%f\n", GetSystematicScalingEtaCorrectionHighMomenta());
5709c40a 2092 printf("SigmaParaThr:\t%f\n", GetSystematicScalingEtaSigmaParaThreshold());
2093 printf("SigmaParaBelowThr:\t%f\n", GetSystematicScalingEtaSigmaParaBelowThreshold());
2094 printf("SigmaParaAboveThr:\t%f\n", GetSystematicScalingEtaSigmaParaAboveThreshold());
e131b05f 2095 printf("MultCorr:\t%f\n", GetSystematicScalingMultCorrection());
71c60990 2096 printf("TOF mode: %d\n", GetTOFmode());
e131b05f 2097
2098 printf("\n\n");
2099}
2100
2101
2102//_____________________________________________________________________________
2103Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePDGcode, Double_t centralityPercentile,
2104 Double_t jetPt)
2105{
2106 // Process the track (generate expected response, fill histos, etc.).
2107 // particlePDGcode == 0 means data. Otherwise, the corresponding MC ID will be assumed.
2108
2109 //Printf("Debug: Task %s is starting to process track: dEdx %f, pTPC %f, eta %f, ncl %d\n", GetName(), track->GetTPCsignal(), track->GetTPCmomentum(),
2110 // track->Eta(), track->GetTPCsignalN());
2111
9e95a906 2112 if(fDebug > 1)
2113 printf("File: %s, Line: %d: ProcessTrack\n", (char*)__FILE__, __LINE__);
2114
2115 if (!fDoPID)
2116 return kFALSE;
2117
2118 if(fDebug > 2)
2119 printf("File: %s, Line: %d: ProcessTrack -> Processing started\n", (char*)__FILE__, __LINE__);
2120
e131b05f 2121 const Bool_t isMC = (particlePDGcode == 0) ? kFALSE : kTRUE;
2122
2123 Int_t binMC = -1;
2124
2125 if (isMC) {
2126 if (TMath::Abs(particlePDGcode) == 211) {//Pion
2127 binMC = 3;
2128 }
2129 else if (TMath::Abs(particlePDGcode) == 321) {//Kaon
2130 binMC = 1;
2131 }
2132 else if (TMath::Abs(particlePDGcode) == 2212) {//Proton
2133 binMC = 4;
2134 }
2135 else if (TMath::Abs(particlePDGcode) == 11) {//Electron
2136 binMC = 0;
2137 }
2138 else if (TMath::Abs(particlePDGcode) == 13) {//Muon
2139 binMC = 2;
2140 }
2141 else // In MC-ID case, set to underflow bin such that the response from this track is only used for unidentified signal generation
2142 // or signal generation with PID response and the track is still there (as in data) - e.g. relevant w.r.t. deuterons.
2143 // This is important to be as much as possible consistent with data. And the tracks can still be removed by disabling the
2144 // underflow bin for the projections
2145 binMC = -1;
2146 }
2147
2148 // Momenta
2149 //Double_t p = track->GetP();
44ed76dd 2150 const Double_t pTPC = track->GetTPCmomentum();
e131b05f 2151 Double_t pT = track->Pt();
2152
2153 Double_t z = -1, xi = -1;
2154 GetJetTrackObservables(pT, jetPt, z, xi);
2155
2156
2157 Double_t trackCharge = track->Charge();
2158
2159 // TPC signal
2160 Double_t dEdxTPC = fPIDResponse->IsTunedOnData() ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
2161
2162 if (dEdxTPC <= 0) {
44ed76dd 2163 Printf("Skipping track with strange dEdx value: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), pTPC,
e131b05f 2164 track->Eta(), track->GetTPCsignalN());
2165 return kFALSE;
2166 }
2167
44ed76dd 2168 // Completely remove tracks from light nuclei defined by el and pr <dEdx> as function of pT.
2169 // A very loose cut to be sure not to throw away electrons and/or protons
2170 Double_t nSigmaPr = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton);
2171 Double_t nSigmaEl = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
e131b05f 2172
44ed76dd 2173 if ((pTPC >= 0.3 && (nSigmaPr > 10 && nSigmaEl > 10)) ||
2174 (pTPC < 0.3 && (nSigmaPr > 15 && nSigmaEl > 15))) {
2175 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",
2176 track->GetTPCsignal(), pTPC, pT, track->Eta(), track->GetTPCsignalN(), nSigmaPr, nSigmaEl);
2177 return kFALSE;
2178 }
e131b05f 2179
2180
2181 Double_t dEdxEl, dEdxKa, dEdxPi, dEdxMu, dEdxPr;
2182 Double_t sigmaEl, sigmaKa, sigmaPi, sigmaMu, sigmaPr;
2183
2184 if (fDoAnySystematicStudiesOnTheExpectedSignal) {
2185 // Get the uncorrected signal first and the corresponding correction factors.
2186 // Then modify the correction factors and properly recalculate the corrected dEdx
2187
2188 // Get pure spline values for dEdx_expected, without any correction
2189 dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2190 dEdxKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2191 dEdxPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2192 dEdxMu = !fTakeIntoAccountMuons ? -1 :
2193 fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2194 dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2195
2196 // Scale splines, if desired
21b0adb1 2197 if ((TMath::Abs(fSystematicScalingSplinesBelowThreshold - 1.0) > fgkEpsilon) ||
2198 (TMath::Abs(fSystematicScalingSplinesAboveThreshold - 1.0) > fgkEpsilon)) {
2199
2200 // Tune turn-on of correction for pions (not so relevant for the other species, since at very large momenta)
21b0adb1 2201 Double_t bg = 0;
2202 Double_t scaleFactor = 1.;
2203 Double_t usedSystematicScalingSplines = 1.;
2204
2205 bg = pTPC / AliPID::ParticleMass(AliPID::kElectron);
2206 scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2207 usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2208 + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2209 dEdxEl *= usedSystematicScalingSplines;
2210
2211 bg = pTPC / AliPID::ParticleMass(AliPID::kKaon);
2212 scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2213 usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2214 + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2215 dEdxKa *= usedSystematicScalingSplines;
2216
2217 bg = pTPC / AliPID::ParticleMass(AliPID::kPion);
2218 scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2219 usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2220 + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2221 dEdxPi *= usedSystematicScalingSplines;
2222
2223 if (fTakeIntoAccountMuons) {
2224 bg = pTPC / AliPID::ParticleMass(AliPID::kMuon);
2225 scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2226 usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2227 + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2228 dEdxMu *= usedSystematicScalingSplines;
2229 }
2230
2231
2232 bg = pTPC / AliPID::ParticleMass(AliPID::kProton);
2233 scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2234 usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2235 + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2236 dEdxPr *= usedSystematicScalingSplines;
e131b05f 2237 }
2238
2239 // Get the eta correction factors for the (modified) expected dEdx
2240 Double_t etaCorrEl = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxEl) : 1.;
2241 Double_t etaCorrKa = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxKa) : 1.;
2242 Double_t etaCorrPi = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPi) : 1.;
2243 Double_t etaCorrMu = fTakeIntoAccountMuons && !fPIDResponse->UseTPCEtaCorrection() ?
2244 fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxMu) : 1.;
2245 Double_t etaCorrPr = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPr) : 1.;
2246
2247 // Scale eta correction factors, if desired (and eta correction maps are to be used, otherwise it is not possible!)
2248 if (fPIDResponse->UseTPCEtaCorrection() &&
2249 (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon ||
2250 TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon)) {
2251 // Since we do not want to scale the splines with this, but only the eta variation, only scale the deviation of the correction factor!
2252 // 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!
2253
2254
2255 // Due to additional azimuthal effects, there is an additional eta dependence for low momenta which is not corrected successfully so far.
2256 // One can assign a different (higher) systematic scale factor for this low-p region and a threshold which separates low- and high-p.
5709c40a 2257 // 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 2258 Double_t usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionHighMomenta;
2259
2260 if (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - fSystematicScalingEtaCorrectionLowMomenta) > fgkEpsilon) {
e131b05f 2261 const Double_t fractionHighMomentumScaleFactor = 0.5 * (1. + TMath::Erf((pTPC - fSystematicScalingEtaCorrectionMomentumThr) / 0.1));
2262 usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionLowMomenta * (1 - fractionHighMomentumScaleFactor)
2263 + fSystematicScalingEtaCorrectionHighMomenta * fractionHighMomentumScaleFactor;
2264 }
2265
2266 etaCorrEl = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrEl - 1.0);
2267 etaCorrKa = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrKa - 1.0);
2268 etaCorrPi = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPi - 1.0);
2269 etaCorrMu = fTakeIntoAccountMuons ? (1.0 + usedSystematicScalingEtaCorrection * (etaCorrMu - 1.0)) : 1.0;
2270 etaCorrPr = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPr - 1.0);
2271 }
2272
2273 // Get the multiplicity correction factors for the (modified) expected dEdx
2274 const Int_t currEvtMultiplicity = fPIDResponse->GetTPCResponse().GetCurrentEventMultiplicity();
2275
2276 Double_t multiplicityCorrEl = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2277 dEdxEl * etaCorrEl, currEvtMultiplicity) : 1.;
2278 Double_t multiplicityCorrKa = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2279 dEdxKa * etaCorrKa, currEvtMultiplicity) : 1.;
2280 Double_t multiplicityCorrPi = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2281 dEdxPi * etaCorrPi, currEvtMultiplicity) : 1.;
2282 Double_t multiplicityCorrMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCMultiplicityCorrection() ?
2283 fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track, dEdxMu * etaCorrMu, currEvtMultiplicity) : 1.;
2284 Double_t multiplicityCorrPr = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2285 dEdxPr * etaCorrPr, currEvtMultiplicity) : 1.;
2286
2287 Double_t multiplicityCorrSigmaEl = fPIDResponse->UseTPCMultiplicityCorrection() ?
2288 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxEl * etaCorrEl, currEvtMultiplicity) : 1.;
2289 Double_t multiplicityCorrSigmaKa = fPIDResponse->UseTPCMultiplicityCorrection() ?
2290 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxKa * etaCorrKa, currEvtMultiplicity) : 1.;
2291 Double_t multiplicityCorrSigmaPi = fPIDResponse->UseTPCMultiplicityCorrection() ?
2292 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxPi * etaCorrPi, currEvtMultiplicity) : 1.;
2293 Double_t multiplicityCorrSigmaMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCMultiplicityCorrection() ?
2294 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxMu * etaCorrMu, currEvtMultiplicity) : 1.;
2295 Double_t multiplicityCorrSigmaPr = fPIDResponse->UseTPCMultiplicityCorrection() ?
2296 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxPr * etaCorrPr, currEvtMultiplicity) : 1.;
2297
2298 // Scale multiplicity correction factors, if desired (and multiplicity correction functions are to be used, otherwise it is not possible!)
2299 if (fPIDResponse->UseTPCMultiplicityCorrection() && TMath::Abs(fSystematicScalingMultCorrection - 1.0) > fgkEpsilon) {
2300 // Since we do not want to scale the splines with this, but only the multiplicity variation, only scale the deviation of the correction factor!
2301 // 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!
2302
2303 multiplicityCorrEl = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrEl - 1.0);
2304 multiplicityCorrKa = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrKa - 1.0);
2305 multiplicityCorrPi = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrPi - 1.0);
2306 multiplicityCorrMu = fTakeIntoAccountMuons ? (1.0 + fSystematicScalingMultCorrection * (multiplicityCorrMu - 1.0)) : 1.0;
2307 multiplicityCorrPr = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrPr - 1.0);
2308 }
2309
2310 // eta correction must be enabled in order to use the new sigma parametrisation maps. Since this is the absolute sigma
2311 // for a track calculated with the unscaled paramaters, we have to devide out dEdxExpectedEtaCorrected and then need
2312 // to scale with the multiplicitySigmaCorrFactor * fSystematicScalingEtaSigmaPara. In the end, one has to scale with the
2313 // (modified) dEdx to get the absolute sigma
2314 // This means there is no extra parameter for the multiplicitySigmaCorrFactor, but only for the sigma map itself.
2315 // This is valid, since it appears only as a product. One has to assume a larger systematic shift in case of additional
2316 // multiplicity dependence....
e131b05f 2317
5709c40a 2318 Bool_t doSigmaSystematics = (TMath::Abs(fSystematicScalingEtaSigmaParaBelowThreshold - 1.0) > fgkEpsilon) ||
2319 (TMath::Abs(fSystematicScalingEtaSigmaParaAboveThreshold - 1.0) > fgkEpsilon);
e131b05f 2320
e131b05f 2321
5709c40a 2322 Double_t dEdxElExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2323 kTRUE, kFALSE);
2324 Double_t systematicScalingEtaSigmaParaEl = 1.;
2325 if (doSigmaSystematics) {
2326 Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxElExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2327 systematicScalingEtaSigmaParaEl = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2328 + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2329 }
2330 Double_t sigmaRelEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2331 kTRUE, kFALSE)
2332 / dEdxElExpected * systematicScalingEtaSigmaParaEl * multiplicityCorrSigmaEl;
2333
2334
2335 Double_t dEdxKaExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2336 kTRUE, kFALSE);
2337 Double_t systematicScalingEtaSigmaParaKa = 1.;
2338 if (doSigmaSystematics) {
2339 Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxKaExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2340 systematicScalingEtaSigmaParaKa = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2341 + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2342 }
2343 Double_t sigmaRelKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2344 kTRUE, kFALSE)
2345 / dEdxKaExpected * systematicScalingEtaSigmaParaKa * multiplicityCorrSigmaKa;
2346
2347
2348 Double_t dEdxPiExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2349 kTRUE, kFALSE);
2350 Double_t systematicScalingEtaSigmaParaPi = 1.;
2351 if (doSigmaSystematics) {
2352 Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxPiExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2353 systematicScalingEtaSigmaParaPi = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2354 + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2355 }
2356 Double_t sigmaRelPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2357 kTRUE, kFALSE)
2358 / dEdxPiExpected * systematicScalingEtaSigmaParaPi * multiplicityCorrSigmaPi;
2359
2360
2361 Double_t sigmaRelMu = 999.;
2362 if (fTakeIntoAccountMuons) {
2363 Double_t dEdxMuExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2364 kTRUE, kFALSE);
2365 Double_t systematicScalingEtaSigmaParaMu = 1.;
2366 if (doSigmaSystematics) {
2367 Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxMuExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2368 systematicScalingEtaSigmaParaMu = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2369 + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2370 }
2371 sigmaRelMu = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2372 / dEdxMuExpected * systematicScalingEtaSigmaParaMu * multiplicityCorrSigmaMu;
2373 }
e131b05f 2374
5709c40a 2375
2376 Double_t dEdxPrExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2377 kTRUE, kFALSE);
2378 Double_t systematicScalingEtaSigmaParaPr = 1.;
2379 if (doSigmaSystematics) {
2380 Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxPrExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2381 systematicScalingEtaSigmaParaPr = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2382 + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2383 }
2384 Double_t sigmaRelPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2385 kTRUE, kFALSE)
2386 / dEdxPrExpected * systematicScalingEtaSigmaParaPr * multiplicityCorrSigmaPr;
e131b05f 2387
2388 // Now scale the (possibly modified) spline values with the (possibly modified) correction factors
2389 dEdxEl *= etaCorrEl * multiplicityCorrEl;
2390 dEdxKa *= etaCorrKa * multiplicityCorrKa;
2391 dEdxPi *= etaCorrPi * multiplicityCorrPi;
2392 dEdxMu *= etaCorrMu * multiplicityCorrMu;
2393 dEdxPr *= etaCorrPr * multiplicityCorrPr;
2394
2395 // Finally, get the absolute sigma
2396 sigmaEl = sigmaRelEl * dEdxEl;
2397 sigmaKa = sigmaRelKa * dEdxKa;
2398 sigmaPi = sigmaRelPi * dEdxPi;
2399 sigmaMu = sigmaRelMu * dEdxMu;
2400 sigmaPr = sigmaRelPr * dEdxPr;
2401
2402 }
2403 else {
2404 // No systematic studies on expected signal - just take it as it comve from the TPCPIDResponse
2405 dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2406 fPIDResponse->UseTPCEtaCorrection(),
2407 fPIDResponse->UseTPCMultiplicityCorrection());
2408 dEdxKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2409 fPIDResponse->UseTPCEtaCorrection(),
2410 fPIDResponse->UseTPCMultiplicityCorrection());
2411 dEdxPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2412 fPIDResponse->UseTPCEtaCorrection(),
2413 fPIDResponse->UseTPCMultiplicityCorrection());
2414 dEdxMu = !fTakeIntoAccountMuons ? -1 :
2415 fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2416 fPIDResponse->UseTPCEtaCorrection(),
2417 fPIDResponse->UseTPCMultiplicityCorrection());
2418 dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2419 fPIDResponse->UseTPCEtaCorrection(),
2420 fPIDResponse->UseTPCMultiplicityCorrection());
2421
2422 sigmaEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2423 fPIDResponse->UseTPCEtaCorrection(),
2424 fPIDResponse->UseTPCMultiplicityCorrection());
2425 sigmaKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2426 fPIDResponse->UseTPCEtaCorrection(),
2427 fPIDResponse->UseTPCMultiplicityCorrection());
2428 sigmaPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2429 fPIDResponse->UseTPCEtaCorrection(),
2430 fPIDResponse->UseTPCMultiplicityCorrection());
2431 sigmaMu = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2432 fPIDResponse->UseTPCEtaCorrection(),
2433 fPIDResponse->UseTPCMultiplicityCorrection());
2434 sigmaPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2435 fPIDResponse->UseTPCEtaCorrection(),
2436 fPIDResponse->UseTPCMultiplicityCorrection());
2437 }
e131b05f 2438
2439 Double_t deltaPrimeElectron = (dEdxEl > 0) ? dEdxTPC / dEdxEl : -1;
2440 if (dEdxEl <= 0) {
2441 Printf("Error: Expected TPC signal <= 0 for electron hypothesis");
2442 return kFALSE;
2443 }
2444
2445 Double_t deltaPrimeKaon = (dEdxKa > 0) ? dEdxTPC / dEdxKa : -1;
2446 if (dEdxKa <= 0) {
2447 Printf("Error: Expected TPC signal <= 0 for kaon hypothesis");
2448 return kFALSE;
2449 }
2450
2451 Double_t deltaPrimePion = (dEdxPi > 0) ? dEdxTPC / dEdxPi : -1;
2452 if (dEdxPi <= 0) {
2453 Printf("Error: Expected TPC signal <= 0 for pion hypothesis");
2454 return kFALSE;
2455 }
2456
2457 Double_t deltaPrimeProton = (dEdxPr > 0) ? dEdxTPC / dEdxPr : -1;
2458 if (dEdxPr <= 0) {
2459 Printf("Error: Expected TPC signal <= 0 for proton hypothesis");
2460 return kFALSE;
2461 }
e131b05f 2462
9e95a906 2463 if(fDebug > 2)
2464 printf("File: %s, Line: %d: ProcessTrack -> Compute probabilities\n", (char*)__FILE__, __LINE__);
e131b05f 2465
44ed76dd 2466
e131b05f 2467 // Use probabilities to weigh the response generation for the different species.
2468 // Also determine the most probable particle type.
2469 Double_t prob[AliPID::kSPECIESC];
2470 for (Int_t i = 0; i < AliPID::kSPECIESC; i++)
2471 prob[i] = 0;
2472
2473 fPIDcombined->ComputeProbabilities(track, fPIDResponse, prob);
2474
44ed76dd 2475 // If muons are not to be taken into account, just set their probability to zero and normalise the remaining probabilities later
e131b05f 2476 if (!fTakeIntoAccountMuons)
2477 prob[AliPID::kMuon] = 0;
2478
44ed76dd 2479 // Priors might cause problems in case of mismatches, i.e. mismatch particles will be identified as pions.
2480 // Can be easily caught for low p_TPC by just setting pion (and muon) probs to zero, if dEdx is larger than
2481 // expected dEdx of electrons and kaons
2482 if (pTPC < 1. && dEdxTPC > dEdxEl && dEdxTPC > dEdxKa) {
2483 prob[AliPID::kMuon] = 0;
2484 prob[AliPID::kPion] = 0;
2485 }
e131b05f 2486
44ed76dd 2487
2488 // Bug: One can set the number of species for PIDcombined, but PIDcombined will call PIDresponse, which writes without testing
2489 // the probs for kSPECIESC (including light nuclei) into the array.
2490 // In this case, when only kSPECIES are considered, the probabilities have to be rescaled!
2491
2492 // Exceptions:
2493 // 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
2494 // high-pT region (also contribution there completely negligable!)
2495 // 2. If dEdx is larger than the expected dEdx + 5 sigma of el and pr, it is most likely a light nucleus
2496 // (if dEdx is more than 5 sigma away from expectation, flat TPC probs are set... and the light nuclei splines are not
2497 // too accurate)
2498 // In these cases:
2499 // The highest expected dEdx is either for pr or for el. Assign 100% probability to the species with highest expected dEdx then.
2500 // In all other cases: Throw away light nuclei probs and rescale other probs to 100%.
2501
2502 // Find most probable species for the ID
2503 Int_t maxIndex = TMath::LocMax(AliPID::kSPECIESC, prob);
2504
2505 if ((prob[maxIndex] > 0 && maxIndex >= AliPID::kSPECIES && dEdxTPC > dEdxEl + 5. * sigmaEl) ||
2506 (dEdxTPC > dEdxEl + 5. * sigmaEl && dEdxTPC > dEdxPr + 5. * sigmaPr)) {
2507 for (Int_t i = 0; i < AliPID::kSPECIESC; i++)
2508 prob[i] = 0;
2509
2510 if (dEdxEl > dEdxPr) {
2511 prob[AliPID::kElectron] = 1.;
2512 maxIndex = AliPID::kElectron;
2513 }
2514 else {
2515 prob[AliPID::kProton] = 1.;
2516 maxIndex = AliPID::kProton;
2517 }
2518 }
2519 else {
2520 for (Int_t i = AliPID::kSPECIES; i < AliPID::kSPECIESC; i++)
2521 prob[i] = 0;
2522
2523 Double_t probSum = 0.;
e131b05f 2524 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
44ed76dd 2525 probSum += prob[i];
2526
2527 if (probSum > 0) {
2528 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2529 prob[i] /= probSum;
2530 }
2531
2532 // If all probabilities are equal (should only happen if no priors are used), set most probable PID to pion
2533 // because LocMax returns just the first maximal value (i.e. AliPID::kElectron)
2534 if (TMath::Abs(prob[maxIndex] - 1. / AliPID::kSPECIES) < 1e-5)
2535 maxIndex = AliPID::kPion;
2536
2537 // In all other cases, the maximum remains untouched from the scaling (and is < AliPID::kSPECIES by construction)
e131b05f 2538 }
2539
2540 if (!isMC) {
44ed76dd 2541 // Clean up the most probable PID for low momenta (not an issue for the probabilities, since fractions small,
2542 // but to get clean (and nice looking) templates (if most probable PID is used instead of generated responses), it should be done).
2543 // Problem: If more than 5 sigma away (and since priors also included), most probable PID for dEdx around protons could be pions.
2544 // Idea: Only clean selection for K and p possible. So, just take for the most probable PID the probs from TPC only calculated
2545 // by hand without restriction to 5 sigma and without priors. This will not affect the template generation, but afterwards one
2546 // can replace the K and p templates with the most probable PID distributions, so all other species templates remain the same.
2547 // I.e. it doesn't hurt if the most probable PID is distributed incorrectly between el, pi, mu.
2548 Bool_t maxIndexSetManually = kFALSE;
2549
2550 if (pTPC < 1.) {
2551 Double_t probManualTPC[AliPID::kSPECIES];
2552 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2553 probManualTPC[i] = 0;
2554
2555 probManualTPC[AliPID::kElectron] = TMath::Exp(-0.5*(dEdxTPC-dEdxEl)*(dEdxTPC-dEdxEl)/(sigmaEl*sigmaEl))/sigmaEl;
2556 // Muons are not important here, just ignore and save processing time
2557 probManualTPC[AliPID::kMuon] = 0.;//TMath::Exp(-0.5*(dEdxTPC-dEdxMu)*(dEdxTPC-dEdxMu)/(sigmaMu*sigmaMu))/sigmaMu;
2558 probManualTPC[AliPID::kPion] = TMath::Exp(-0.5*(dEdxTPC-dEdxPi)*(dEdxTPC-dEdxPi)/(sigmaPi*sigmaPi))/sigmaPi;
2559 probManualTPC[AliPID::kKaon] = TMath::Exp(-0.5*(dEdxTPC-dEdxKa)*(dEdxTPC-dEdxKa)/(sigmaKa*sigmaKa))/sigmaKa;
2560 probManualTPC[AliPID::kProton] = TMath::Exp(-0.5*(dEdxTPC-dEdxPr)*(dEdxTPC-dEdxPr)/(sigmaPr*sigmaPr))/sigmaPr;
2561
2562 const Int_t maxIndexManualTPC = TMath::LocMax(AliPID::kSPECIES, probManualTPC);
2563 // Should always be > 0, but in case the dEdx is far off such that the machine accuracy sets every prob to zero,
2564 // better take the "old" result
2565 if (probManualTPC[maxIndexManualTPC] > 0.)
2566 maxIndex = maxIndexManualTPC;
2567
2568 maxIndexSetManually = kTRUE;
e131b05f 2569 }
44ed76dd 2570
2571
e131b05f 2572 // Translate from AliPID numbering to numbering of this class
44ed76dd 2573 if (prob[maxIndex] > 0 || maxIndexSetManually) {
e131b05f 2574 if (maxIndex == AliPID::kElectron)
2575 binMC = 0;
2576 else if (maxIndex == AliPID::kKaon)
2577 binMC = 1;
2578 else if (maxIndex == AliPID::kMuon)
2579 binMC = 2;
2580 else if (maxIndex == AliPID::kPion)
2581 binMC = 3;
2582 else if (maxIndex == AliPID::kProton)
2583 binMC = 4;
2584 else
2585 binMC = -1;
2586 }
2587 else {
2588 // Only take track into account for expectation values, if valid pid response is available.. Otherwise: Set to underflow bin.
2589 binMC = -1;
2590 }
2591 }
2592
2593 /*
2594 //For testing: Swap p<->pT to analyse pure p-dependence => Needs to be removed later
2595 Double_t temp = pT;
2596 pT = pTPC;
2597 pTPC = temp;
2598 */
2599
77324970 2600 TOFpidInfo tofPIDinfo = GetTOFType(track, fTOFmode);
e131b05f 2601
2602 Double_t entry[fStoreAdditionalJetInformation ? kDataNumAxes : kDataNumAxes - fgkNumJetAxes];
2603 entry[kDataMCID] = binMC;
2604 entry[kDataSelectSpecies] = 0;
2605 entry[kDataPt] = pT;
2606 entry[kDataDeltaPrimeSpecies] = deltaPrimeElectron;
2607 entry[kDataCentrality] = centralityPercentile;
2608
2609 if (fStoreAdditionalJetInformation) {
2610 entry[kDataJetPt] = jetPt;
2611 entry[kDataZ] = z;
2612 entry[kDataXi] = xi;
2613 }
2614
2615 entry[GetIndexOfChargeAxisData()] = trackCharge;
77324970 2616 entry[GetIndexOfTOFpidInfoAxisData()] = tofPIDinfo;
e131b05f 2617
2618 fhPIDdataAll->Fill(entry);
2619
2620 entry[kDataSelectSpecies] = 1;
e131b05f 2621 entry[kDataDeltaPrimeSpecies] = deltaPrimeKaon;
e131b05f 2622 fhPIDdataAll->Fill(entry);
2623
2624 entry[kDataSelectSpecies] = 2;
e131b05f 2625 entry[kDataDeltaPrimeSpecies] = deltaPrimePion;
e131b05f 2626 fhPIDdataAll->Fill(entry);
2627
2628 entry[kDataSelectSpecies] = 3;
e131b05f 2629 entry[kDataDeltaPrimeSpecies] = deltaPrimeProton;
e131b05f 2630 fhPIDdataAll->Fill(entry);
2631
2632
2633 // Construct the expected shape for the transition p -> pT
2634
2635 Double_t genEntry[fStoreAdditionalJetInformation ? kGenNumAxes : kGenNumAxes - fgkNumJetAxes];
2636 genEntry[kGenMCID] = binMC;
2637 genEntry[kGenSelectSpecies] = 0;
2638 genEntry[kGenPt] = pT;
2639 genEntry[kGenDeltaPrimeSpecies] = -999;
2640 genEntry[kGenCentrality] = centralityPercentile;
2641
2642 if (fStoreAdditionalJetInformation) {
2643 genEntry[kGenJetPt] = jetPt;
2644 genEntry[kGenZ] = z;
2645 genEntry[kGenXi] = xi;
2646 }
2647
2648 genEntry[GetIndexOfChargeAxisGen()] = trackCharge;
77324970 2649 genEntry[GetIndexOfTOFpidInfoAxisGen()] = tofPIDinfo;
e131b05f 2650
77324970
CKB		 2651 // Generate numGenEntries "responses" with fluctuations around the expected values.
2652 // fgkMaxNumGenEntries = 500 turned out to give reasonable templates even for highest track pT in 15-20 GeV/c jet pT bin.
2653 Int_t numGenEntries = fgkMaxNumGenEntries; // fgkMaxNumGenEntries = 500
e131b05f 2654
77324970
CKB		 2655 /*OLD: Different number of responses depending on pT and jet pT for fgkMaxNumGenEntries = 1000
2656 * => Problem: If threshold to higher number of responses inside a bin (or after rebinning), then the template
2657 * is biased to the higher pT.
e131b05f 2658 // Generate numGenEntries "responses" with fluctuations around the expected values.
2659 // The higher the (transverse) momentum, the more "responses" will be generated in order not to run out of statistics too fast.
2660 Int_t numGenEntries = 10;
77324970 2661
e131b05f 2662 // Jets have even worse statistics, therefore, scale numGenEntries further
2663 if (jetPt >= 40)
2664 numGenEntries *= 20;
2665 else if (jetPt >= 20)
2666 numGenEntries *= 10;
2667 else if (jetPt >= 10)
2668 numGenEntries *= 2;
2669
2670
77324970 2671
e131b05f 2672 // Do not generate more entries than available in memory!
2673 if (numGenEntries > fgkMaxNumGenEntries)// fgkMaxNumGenEntries = 1000
2674 numGenEntries = fgkMaxNumGenEntries;
77324970
CKB		 2675 */
2676
2677
e131b05f 2678 ErrorCode errCode = kNoErrors;
2679
9e95a906 2680 if(fDebug > 2)
2681 printf("File: %s, Line: %d: ProcessTrack -> Generate Responses\n", (char*)__FILE__, __LINE__);
2682
e131b05f 2683 // Electrons
2684 errCode = GenerateDetectorResponse(errCode, 1., sigmaEl / dEdxEl, fGenRespElDeltaPrimeEl, numGenEntries);
2685 errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxKa, sigmaEl / dEdxKa, fGenRespElDeltaPrimeKa, numGenEntries);
2686 errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxPi, sigmaEl / dEdxPi, fGenRespElDeltaPrimePi, numGenEntries);
2687 errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxPr, sigmaEl / dEdxPr, fGenRespElDeltaPrimePr, numGenEntries);
2688
2689 // Kaons
2690 errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxEl, sigmaKa / dEdxEl, fGenRespKaDeltaPrimeEl, numGenEntries);
2691 errCode = GenerateDetectorResponse(errCode, 1., sigmaKa / dEdxKa, fGenRespKaDeltaPrimeKa, numGenEntries);
2692 errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxPi, sigmaKa / dEdxPi, fGenRespKaDeltaPrimePi, numGenEntries);
2693 errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxPr, sigmaKa / dEdxPr, fGenRespKaDeltaPrimePr, numGenEntries);
2694
2695 // Pions
2696 errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxEl, sigmaPi / dEdxEl, fGenRespPiDeltaPrimeEl, numGenEntries);
2697 errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxKa, sigmaPi / dEdxKa, fGenRespPiDeltaPrimeKa, numGenEntries);
2698 errCode = GenerateDetectorResponse(errCode, 1., sigmaPi / dEdxPi, fGenRespPiDeltaPrimePi, numGenEntries);
2699 errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxPr, sigmaPi / dEdxPr, fGenRespPiDeltaPrimePr, numGenEntries);
2700
2701 // Muons, if desired
2702 if (fTakeIntoAccountMuons) {
2703 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxEl, sigmaMu / dEdxEl, fGenRespMuDeltaPrimeEl, numGenEntries);
2704 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxKa, sigmaMu / dEdxKa, fGenRespMuDeltaPrimeKa, numGenEntries);
2705 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxPi, sigmaMu / dEdxPi, fGenRespMuDeltaPrimePi, numGenEntries);
2706 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxPr, sigmaMu / dEdxPr, fGenRespMuDeltaPrimePr, numGenEntries);
2707 }
2708
2709 // Protons
2710 errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxEl, sigmaPr / dEdxEl, fGenRespPrDeltaPrimeEl, numGenEntries);
2711 errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxKa, sigmaPr / dEdxKa, fGenRespPrDeltaPrimeKa, numGenEntries);
2712 errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxPi, sigmaPr / dEdxPi, fGenRespPrDeltaPrimePi, numGenEntries);
2713 errCode = GenerateDetectorResponse(errCode, 1., sigmaPr / dEdxPr, fGenRespPrDeltaPrimePr, numGenEntries);
2714
e131b05f 2715 if (errCode != kNoErrors) {
77324970
CKB		 2716 if (errCode == kWarning && fDebug > 1) {
2717 Printf("Warning: Questionable detector response for track, most likely due to very low number of PID clusters! Debug output (dEdx_expected, sigma_expected):");
e131b05f 2718 }
2719 else
2720 Printf("Error: Failed to generate detector response for track - skipped! Debug output (dEdx_expected, sigma_expected):");
2721
77324970
CKB		 2722 if (fDebug > 1) {
2723 Printf("Pr: %e, %e", dEdxPr, sigmaPr);
2724 Printf("Pi: %e, %e", dEdxPi, sigmaPi);
2725 Printf("El: %e, %e", dEdxEl, sigmaEl);
2726 Printf("Mu: %e, %e", dEdxMu, sigmaMu);
2727 Printf("Ka: %e, %e", dEdxKa, sigmaKa);
2728 Printf("track: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), track->GetTPCmomentum(), track->Eta(),
2729 track->GetTPCsignalN());
2730 }
e131b05f 2731
2732 if (errCode != kWarning) {
2733 fhSkippedTracksForSignalGeneration->Fill(track->GetTPCmomentum(), track->GetTPCsignalN());
2734 return kFALSE;// Skip generated response in case of error
2735 }
2736 }
2737
2738 for (Int_t n = 0; n < numGenEntries; n++) {
2739 if (!isMC || !fUseMCidForGeneration) {
2740 // If no MC info is available or shall not be used, use weighting with priors to generate entries for the different species
2741 Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
2742
2743 // Consider generated response as originating from...
2744 if (rnd <= prob[AliPID::kElectron])
2745 genEntry[kGenMCID] = 0; // ... an electron
2746 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon])
2747 genEntry[kGenMCID] = 1; // ... a kaon
2748 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon])
2749 genEntry[kGenMCID] = 2; // ... a muon -> NOTE: prob[AliPID::kMuon] = 0 in case of fTakeIntoAccountMuons = kFALSE
2750 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon] + prob[AliPID::kPion])
2751 genEntry[kGenMCID] = 3; // ... a pion
2752 else
2753 genEntry[kGenMCID] = 4; // ... a proton
2754 }
2755
2756 // Electrons
2757 genEntry[kGenSelectSpecies] = 0;
e131b05f 2758 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeEl[n];
2759 fhGenEl->Fill(genEntry);
2760
2761 genEntry[kGenSelectSpecies] = 1;
e131b05f 2762 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeKa[n];
2763 fhGenEl->Fill(genEntry);
2764
2765 genEntry[kGenSelectSpecies] = 2;
e131b05f 2766 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePi[n];
2767 fhGenEl->Fill(genEntry);
2768
2769 genEntry[kGenSelectSpecies] = 3;
e131b05f 2770 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePr[n];
2771 fhGenEl->Fill(genEntry);
2772
2773 // Kaons
2774 genEntry[kGenSelectSpecies] = 0;
e131b05f 2775 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeEl[n];
2776 fhGenKa->Fill(genEntry);
2777
2778 genEntry[kGenSelectSpecies] = 1;
e131b05f 2779 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeKa[n];
2780 fhGenKa->Fill(genEntry);
2781
2782 genEntry[kGenSelectSpecies] = 2;
e131b05f 2783 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePi[n];
2784 fhGenKa->Fill(genEntry);
2785
2786 genEntry[kGenSelectSpecies] = 3;
e131b05f 2787 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePr[n];
2788 fhGenKa->Fill(genEntry);
2789
2790 // Pions
2791 genEntry[kGenSelectSpecies] = 0;
e131b05f 2792 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeEl[n];
2793 fhGenPi->Fill(genEntry);
2794
2795 genEntry[kGenSelectSpecies] = 1;
e131b05f 2796 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeKa[n];
2797 fhGenPi->Fill(genEntry);
2798
2799 genEntry[kGenSelectSpecies] = 2;
e131b05f 2800 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePi[n];
2801 fhGenPi->Fill(genEntry);
2802
2803 genEntry[kGenSelectSpecies] = 3;
e131b05f 2804 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePr[n];
2805 fhGenPi->Fill(genEntry);
2806
2807 if (fTakeIntoAccountMuons) {
2808 // Muons, if desired
2809 genEntry[kGenSelectSpecies] = 0;
e131b05f 2810 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeEl[n];
2811 fhGenMu->Fill(genEntry);
2812
2813 genEntry[kGenSelectSpecies] = 1;
e131b05f 2814 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeKa[n];
2815 fhGenMu->Fill(genEntry);
2816
2817 genEntry[kGenSelectSpecies] = 2;
e131b05f 2818 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePi[n];
2819 fhGenMu->Fill(genEntry);
2820
2821 genEntry[kGenSelectSpecies] = 3;
e131b05f 2822 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePr[n];
2823 fhGenMu->Fill(genEntry);
2824 }
2825
2826 // Protons
2827 genEntry[kGenSelectSpecies] = 0;
e131b05f 2828 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeEl[n];
2829 fhGenPr->Fill(genEntry);
2830
2831 genEntry[kGenSelectSpecies] = 1;
e131b05f 2832 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeKa[n];
2833 fhGenPr->Fill(genEntry);
2834
2835 genEntry[kGenSelectSpecies] = 2;
e131b05f 2836 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePi[n];
2837 fhGenPr->Fill(genEntry);
2838
2839 genEntry[kGenSelectSpecies] = 3;
e131b05f 2840 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePr[n];
2841 fhGenPr->Fill(genEntry);
2842 }
2843
9e95a906 2844 if(fDebug > 2)
2845 printf("File: %s, Line: %d: ProcessTrack -> Done\n", (char*)__FILE__, __LINE__);
2846
e131b05f 2847 return kTRUE;
2848}
2849
2850
2851//_____________________________________________________________________________
2852Bool_t AliAnalysisTaskPID::SetConvolutedGaussLambdaParameter(Double_t lambda)
2853{
2854 // Set the lambda parameter of the convolution to the desired value. Automatically
2855 // calculates the parameters for the transition (restricted) gauss -> convoluted gauss.
2856 fConvolutedGaussTransitionPars[2] = lambda;
2857
2858 // Save old parameters and settings of function to restore them later:
2859 Double_t* oldFuncParams = new Double_t[fkConvolutedGausNPar];
2860 for (Int_t i = 0; i < fkConvolutedGausNPar; i++)
2861 oldFuncParams[i] = fConvolutedGausDeltaPrime->GetParameter(i);
2862 Int_t oldFuncNpx = fConvolutedGausDeltaPrime->GetNpx();
2863 Double_t oldFuncRangeLow = 0, oldFuncRangeUp = 100;
2864 fConvolutedGausDeltaPrime->GetRange(oldFuncRangeLow, oldFuncRangeUp);
2865
2866 // Choose some sufficiently large range
2867 const Double_t rangeStart = 0.5;
2868 const Double_t rangeEnd = 2.0;
2869
2870 // To get the parameters for the transition, just choose arbitrarily input parameters for mu and sigma
2871 // (it makes sense to choose typical values). The ratio sigma_gauss / sigma_convolution is independent
2872 // of mu and as well the difference mu_gauss - mu_convolution.
2873 Double_t muInput = 1.0;
2874 Double_t sigmaInput = fgkSigmaReferenceForTransitionPars;
2875
2876
2877 // Step 1: Generate distribution with input parameters
2878 const Int_t nPar = 3;
2879 Double_t inputPar[nPar] = { muInput, sigmaInput, lambda };
2880
2881 TH1D* hInput = new TH1D("hInput", "Input distribution", 2000, rangeStart, rangeEnd);
2882
2883 fConvolutedGausDeltaPrime->SetParameters(inputPar);
2884 fConvolutedGausDeltaPrime->SetRange(rangeStart, rangeEnd);
2885 fConvolutedGausDeltaPrime->SetNpx(2000);
2886
2887 /*OLD
2888 // The resolution and mean of the AliPIDResponse are determined in finite intervals
2889 // of ncl (also second order effects due to finite dEdx and tanTheta).
2890 // Take this into account for the transition parameters via assuming an approximately flat
2891 // distritubtion in ncl in this interval.
2892 // NOTE: The ncl interval should be the same as the one used for the sigma map creation!
2893 const Int_t nclStart = 151;
2894 const Int_t nclEnd = 160;
2895 const Int_t nclSteps = (nclEnd - nclStart) + 1;
2896 for (Int_t ncl = nclStart; ncl <= nclEnd; ncl++) {
2897 // Resolution scales with 1/sqrt(ncl)
2898 fConvolutedGausDeltaPrime->SetParameter(1, inputPar[1] * sqrt(nclEnd) / sqrt(ncl));
2899 TH1* hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
2900
2901 for (Int_t i = 0; i < 50000000 / nclSteps; i++)
2902 hInput->Fill(hProbDensity->GetRandom());
2903 }
2904 */
2905
2906 TH1* hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
2907
2908 for (Int_t i = 0; i < 50000000; i++)
2909 hInput->Fill(hProbDensity->GetRandom());
2910
2911 // Step 2: Fit generated distribution with restricted gaussian
2912 Int_t maxBin = hInput->GetMaximumBin();
2913 Double_t max = hInput->GetBinContent(maxBin);
2914
2915 UChar_t usedBins = 1;
2916 if (maxBin > 1) {
2917 max += hInput->GetBinContent(maxBin - 1);
2918 usedBins++;
2919 }
2920 if (maxBin < hInput->GetNbinsX()) {
2921 max += hInput->GetBinContent(maxBin + 1);
2922 usedBins++;
2923 }
2924 max /= usedBins;
2925
2926 // NOTE: The range (<-> fraction of maximum) should be the same
2927 // as for the spline and eta maps creation
2928 const Double_t lowThreshold = hInput->GetXaxis()->GetBinLowEdge(hInput->FindFirstBinAbove(0.1 * max));
2929 const Double_t highThreshold = hInput->GetXaxis()->GetBinUpEdge(hInput->FindLastBinAbove(0.1 * max));
2930
2931 TFitResultPtr fitResGaussFirstStep = hInput->Fit("gaus", "QNRS", "", lowThreshold, highThreshold);
2932
2933 TFitResultPtr fitResGauss;
2934
2935 if ((Int_t)fitResGaussFirstStep == 0) {
2936 TF1 fGauss("fGauss", "[0]*TMath::Gaus(x, [1], [2], 1)", rangeStart, rangeEnd);
2937 fGauss.SetParameter(0, fitResGaussFirstStep->GetParams()[0]);
2938 fGauss.SetParError(0, fitResGaussFirstStep->GetErrors()[0]);
2939 fGauss.SetParameter(2, fitResGaussFirstStep->GetParams()[2]);
2940 fGauss.SetParError(2, fitResGaussFirstStep->GetErrors()[2]);
2941
2942 fGauss.FixParameter(1, fitResGaussFirstStep->GetParams()[1]);
2943 fitResGauss = hInput->Fit(&fGauss, "QNS", "s", rangeStart, rangeEnd);
2944 }
2945 else {
2946 fitResGauss = hInput->Fit("gaus", "QNRS", "same", rangeStart, rangeEnd);
2947 }
2948 //OLD TFitResultPtr fitResGauss = hInput->Fit("gaus", "QNRS", "", hInput->GetXaxis()->GetBinLowEdge(hInput->FindFirstBinAbove(0.1 * max)),
2949 // hInput->GetXaxis()->GetBinUpEdge(hInput->FindLastBinAbove(0.1 * max)));
2950
2951
2952 // Step 3: Use parameters from gaussian fit to obtain parameters for the transition "restricted gauss" -> "convoluted gauss"
2953
2954 // 3.1 The ratio sigmaInput / sigma_gaussFit ONLY depends on lambda (which is fixed per period) -> Calculate this first
2955 // for an arbitrary (here: typical) sigma. The ratio is then ~the same for ALL sigma for given lambda!
2956 if ((Int_t)fitResGauss != 0) {
2957 AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Gauss Fit failed!\n");
2958 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2959 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2960 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2961
2962 delete hInput;
c4856fb1 2963 delete [] oldFuncParams;
e131b05f 2964
2965 return kFALSE;
2966 }
2967
2968 if (fitResGauss->GetParams()[2] <= 0) {
2969 AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Sigma of gauss fit <= 0!\n");
2970 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2971 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2972 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2973
2974 delete hInput;
c4856fb1 2975 delete [] oldFuncParams;
e131b05f 2976
2977 return kFALSE;
2978 }
2979
2980 // sigma correction factor
2981 fConvolutedGaussTransitionPars[1] = sigmaInput / fitResGauss->GetParams()[2];
2982
2983 // 3.2 Now that sigma und lambda are determined, one can calculate mu by shifting the maximum to the desired position,
2984 // i.e. the maximum of the input distribution should coincide with that of the re-calculated distribution,
2985 // which is achieved by getting the same mu for the same sigma.
2986 // NOTE: For fixed lambda, the shift is proportional to sigma and independent of mu!
2987 // So, one can calculate the shift for an arbitrary fixed (here: typical)
2988 // sigma and then simply use this shift for any other sigma by scaling it correspondingly!!!
2989
2990 // Mu shift correction:
2991 // Shift in mu (difference between mean of gauss and mean of convolution) is proportional to sigma!
2992 // Thus, choose a reference sigma (typical -> 0.05), such that for arbitrary sigma one can simple scale
2993 // this shift correction with sigma / referenceSigma.
2994 fConvolutedGaussTransitionPars[0] = (fitResGauss->GetParams()[1] - muInput);
2995
2996
2997 /*Changed on 03.07.2013
2998
2999 // Maximum of fConvolutedGausDeltaPrime should agree with maximum of input
3000 Double_t par[nPar] = { fitResGauss->GetParams()[1], // just as a guess of the maximum position
3001 sigmaInput,
3002 lambda };
3003
3004 fConvolutedGausDeltaPrime->SetParameters(par);
3005
3006 Double_t maxXInput = fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001, muInput - 3. * sigmaInput),
3007 muInput + 10. * sigmaInput,
3008 0.0001);
3009
3010 // Maximum shifts proportional to sigma and is linear in mu (~mean of gauss)
3011 // Maximum should be typically located within [gaussMean, gaussMean + 3 gaussSigma].
3012 // -> Larger search range for safety reasons (also: sigma and/or mean might be not 100% accurate).
3013 Double_t maxXconvoluted = fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001,
3014 fitResGauss->GetParams()[1] - 3. * fitResGauss->GetParams()[2]),
3015 fitResGauss->GetParams()[1] + 10. * fitResGauss->GetParams()[2],
3016 0.0001);
3017 if (maxXconvoluted <= 0) {
3018 AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Maximum of fConvolutedGausDeltaPrime <= 0!\n");
3019
3020 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
3021 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
3022 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
3023
3024 delete hInput;
c4856fb1 3025 delete [] oldFuncParams;
e131b05f 3026
3027 return kFALSE;
3028 }
3029
3030 // maxX perfectly shifts as par[0] (scaled by sigma) -> Can shift maxX to input value.
3031 // Mu shift correction:
3032 fConvolutedGaussTransitionPars[0] = maxXconvoluted - maxXInput;
3033 */
3034
3035
3036
3037 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
3038 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
3039 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
3040
3041 delete hInput;
c4856fb1 3042 delete [] oldFuncParams;
e131b05f 3043
3044 return kTRUE;
3045}
3046
3047
3048//_____________________________________________________________________________
9d7ad2e4 3049AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::SetParamsForConvolutedGaus(Double_t gausMean, Double_t gausSigma)
e131b05f 3050{
3051 // Set parameters for convoluted gauss using parameters for a pure gaussian.
3052 // If SetConvolutedGaussLambdaParameter has not been called before to initialise the translation parameters,
3053 // some default parameters will be used and an error will show up.
3054
9e95a906 3055 if(fDebug > 1)
3056 printf("File: %s, Line: %d: SetParamsForConvolutedGaus: mean %e, sigma %e\n", (char*)__FILE__, __LINE__, gausMean, gausSigma);
3057
e131b05f 3058 if (fConvolutedGaussTransitionPars[1] < -998) {
3059 AliError("Transition parameters not initialised! Default parameters will be used. Please call SetConvolutedGaussLambdaParameter(...) before any calculations!");
3060 SetConvolutedGaussLambdaParameter(2.0);
3061 AliError(Form("Parameters set to:\n[0]: %f\n[1]: %f\n[2]: %f\n", fConvolutedGaussTransitionPars[0],
3062 fConvolutedGaussTransitionPars[1], fConvolutedGaussTransitionPars[2]));
3063 }
3064
3065 Double_t par[fkConvolutedGausNPar];
3066 par[2] = fConvolutedGaussTransitionPars[2];
3067 par[1] = fConvolutedGaussTransitionPars[1] * gausSigma;
3068 // maxX perfectly shifts as par[0] (scaled by sigma) -> Can shift maxX so that it sits at the right place.
3069 par[0] = gausMean - fConvolutedGaussTransitionPars[0] * par[1] / fgkSigmaReferenceForTransitionPars;
3070
3071 ErrorCode errCode = kNoErrors;
3072 fConvolutedGausDeltaPrime->SetParameters(par);
3073
9e95a906 3074 if(fDebug > 3)
3075 printf("File: %s, Line: %d: SetParamsForConvolutedGaus -> Parameters set to: %e, %e, %e (transition pars: %e, %e, %e, %e)\n",
3076 (char*)__FILE__, __LINE__, par[0], par[1], par[2], fConvolutedGaussTransitionPars[0], fConvolutedGaussTransitionPars[1],
3077 fConvolutedGaussTransitionPars[2], fgkSigmaReferenceForTransitionPars);
3078
e131b05f 3079 fConvolutedGausDeltaPrime->SetNpx(20); // Small value speeds up following algorithm (valid, since extrema far apart)
3080
3081 // Accuracy of 10^-5 is enough to get 0.1% precise peak for MIPS w.r.t. to dEdx = 2000 of protons
3082 // (should boost up the algorithm, because 10^-10 is the default value!)
3083 Double_t maxX= fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001, gausMean - 2. * gausSigma),
3084 gausMean + 6. * gausSigma, 1.0E-5);
3085
3086 const Double_t maximum = fConvolutedGausDeltaPrime->Eval(maxX);
3087 const Double_t maximumFraction = maximum * fAccuracyNonGaussianTail;
3088
3089 // Estimate lower boundary for subsequent search:
3090 Double_t lowBoundSearchBoundLow = TMath::Max(1e-4, maxX - 5. * gausSigma);
3091 Double_t lowBoundSearchBoundUp = maxX;
3092
3093 Bool_t lowerBoundaryFixedAtZero = kFALSE;
3094
3095 while (fConvolutedGausDeltaPrime->Eval(lowBoundSearchBoundLow) >= maximumFraction) {
3096 if (lowBoundSearchBoundLow <= 0) {
3097 // This should only happen to low dEdx particles with very few clusters and therefore large sigma, such that the gauss goes below zero deltaPrime
3098 if (maximum <= 0) { // Something is weired
3099 printf("Error generating signal: maximum is <= 0!\n");
3100 return kError;
3101 }
3102 else {
3103 const Double_t valueAtZero = fConvolutedGausDeltaPrime->Eval(0);
3104 if (valueAtZero / maximum > 0.05) {
3105 // Too large fraction below zero deltaPrime. Signal generation cannot be reliable in this case
3106 printf("Error generating signal: Too large fraction below zero deltaPrime: convGauss(0) / convGauss(max) = %e / %e = %e!\n",
3107 valueAtZero, maximum, valueAtZero / maximum);
3108 return kError;
3109 }
3110 }
3111
3112 /*
3113 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",
3114 fConvolutedGausDeltaPrime->Eval(maxX), fAccuracyNonGaussianTail, maximumFraction, maxX, par[0], par[1], par[2], gausMean, gausSigma);
3115 */
3116
3117 lowerBoundaryFixedAtZero = kTRUE;
3118
3119 if (errCode != kError)
3120 errCode = kWarning;
3121
3122 break;
3123 }
3124
3125 lowBoundSearchBoundUp -= gausSigma;
3126 lowBoundSearchBoundLow -= gausSigma;
3127
3128 if (lowBoundSearchBoundLow < 0) {
3129 lowBoundSearchBoundLow = 0;
3130 lowBoundSearchBoundUp += gausSigma;
3131 }
3132 }
3133
3134 // Determine lower boundary inside estimated range. For small values of the maximum: Need more precision, since finer binning!
3135 Double_t rangeStart = lowerBoundaryFixedAtZero ? 0 :
3136 fConvolutedGausDeltaPrime->GetX(maximumFraction, lowBoundSearchBoundLow, lowBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
3137
3138 // .. and the same for the upper boundary
3139 Double_t rangeEnd = 0;
3140 // If distribution starts beyond upper boundary, everything ends up in the overflow bin. So, just reduce range and Npx to minimum
3141 if (rangeStart > fkDeltaPrimeUpLimit) {
3142 rangeEnd = rangeStart + 0.00001;
3143 fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
3144 fConvolutedGausDeltaPrime->SetNpx(4);
3145 }
3146 else {
3147 // Estimate upper boundary for subsequent search:
3148 Double_t upBoundSearchBoundUp = maxX + 5 * gausSigma;
3149 Double_t upBoundSearchBoundLow = maxX;
3150 while (fConvolutedGausDeltaPrime->Eval(upBoundSearchBoundUp) >= maximumFraction) {
3151 upBoundSearchBoundUp += gausSigma;
3152 upBoundSearchBoundLow += gausSigma;
3153 }
3154
3155 // For small values of the maximum: Need more precision, since finer binning!
3156 rangeEnd = fConvolutedGausDeltaPrime->GetX(maximumFraction, upBoundSearchBoundLow, upBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
3157
3158 fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
3159 fConvolutedGausDeltaPrime->SetNpx(fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindBin(rangeEnd)
3160 - fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindBin(rangeStart) + 1);
3161 //fConvolutedGausDeltaPrime->SetNpx((rangeEnd - rangeStart) / fDeltaPrimeBinWidth + 1);
3162 }
3163
9e95a906 3164 if(fDebug > 3)
3165 printf("File: %s, Line: %d: SetParamsForConvolutedGaus -> range %f - %f, error code %d\n", (char*)__FILE__, __LINE__,
3166 rangeStart, rangeEnd, errCode);
3167
e131b05f 3168 return errCode;
3169}
3170
3171
3172//________________________________________________________________________
3173void AliAnalysisTaskPID::SetUpGenHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const
3174{
3175 // Sets bin limits for axes which are not standard binned and the axes titles.
3176
3177 hist->SetBinEdges(kGenPt, binsPt);
3178 hist->SetBinEdges(kGenDeltaPrimeSpecies, binsDeltaPrime);
3179 hist->SetBinEdges(kGenCentrality, binsCent);
3180
3181 if (fStoreAdditionalJetInformation)
3182 hist->SetBinEdges(kGenJetPt, binsJetPt);
3183
3184 // Set axes titles
3185 hist->GetAxis(kGenMCID)->SetTitle("MC PID");
3186 hist->GetAxis(kGenMCID)->SetBinLabel(1, "e");
3187 hist->GetAxis(kGenMCID)->SetBinLabel(2, "K");
3188 hist->GetAxis(kGenMCID)->SetBinLabel(3, "#mu");
3189 hist->GetAxis(kGenMCID)->SetBinLabel(4, "#pi");
3190 hist->GetAxis(kGenMCID)->SetBinLabel(5, "p");
3191
3192 hist->GetAxis(kGenSelectSpecies)->SetTitle("Select Species");
3193 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(1, "e");
3194 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(2, "K");
3195 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(3, "#pi");
3196 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(4, "p");
3197
5709c40a 3198 hist->GetAxis(kGenPt)->SetTitle("p_{T} (GeV/c)");
e131b05f 3199
3200 hist->GetAxis(kGenDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
3201
3202 hist->GetAxis(kGenCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
3203
3204 if (fStoreAdditionalJetInformation) {
5709c40a 3205 hist->GetAxis(kGenJetPt)->SetTitle("p_{T}^{jet} (GeV/c)");
e131b05f 3206
5709c40a 3207 hist->GetAxis(kGenZ)->SetTitle("z = p_{T}^{track} / p_{T}^{jet}");
e131b05f 3208
5709c40a 3209 hist->GetAxis(kGenXi)->SetTitle("#xi = ln(p_{T}^{jet} / p_{T}^{track})");
e131b05f 3210 }
3211
3212 hist->GetAxis(GetIndexOfChargeAxisGen())->SetTitle("Charge (e_{0})");
77324970
CKB		 3213
3214 hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetTitle("TOF PID Info");
3215 // Offset is (TMath::Abs(kNoTOFinfo) + 1), such that bin 1 (= first label) corresponds to kNoTOFinfo (< 0)
3216 hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kNoTOFinfo + (TMath::Abs(kNoTOFinfo) + 1), "No TOF Info");
3217 hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kNoTOFpid + (TMath::Abs(kNoTOFinfo) + 1), "No TOF PID");
3218 hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kTOFpion + (TMath::Abs(kNoTOFinfo) + 1), "#pi");
3219 hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kTOFkaon + (TMath::Abs(kNoTOFinfo) + 1), "K");
3220 hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kTOFproton + (TMath::Abs(kNoTOFinfo) + 1), "p");
e131b05f 3221}
3222
3223
3224//________________________________________________________________________
3225void AliAnalysisTaskPID::SetUpGenYieldHist(THnSparse* hist, Double_t* binsPt, Double_t* binsCent, Double_t* binsJetPt) const
3226{
3227 // Sets bin limits for axes which are not standard binned and the axes titles.
3228
3229 hist->SetBinEdges(kGenYieldPt, binsPt);
3230 hist->SetBinEdges(kGenYieldCentrality, binsCent);
3231 if (fStoreAdditionalJetInformation)
3232 hist->SetBinEdges(kGenYieldJetPt, binsJetPt);
3233
3234 for (Int_t i = 0; i < 5; i++)
3235 hist->GetAxis(kGenYieldMCID)->SetBinLabel(i + 1, AliPID::ParticleLatexName(i));
3236
3237 // Set axes titles
3238 hist->GetAxis(kGenYieldMCID)->SetTitle("MC PID");
5709c40a 3239 hist->GetAxis(kGenYieldPt)->SetTitle("p_{T}^{gen} (GeV/c)");
e131b05f 3240 hist->GetAxis(kGenYieldCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
3241
3242 if (fStoreAdditionalJetInformation) {
5709c40a 3243 hist->GetAxis(kGenYieldJetPt)->SetTitle("p_{T}^{jet, gen} (GeV/c)");
e131b05f 3244
5709c40a 3245 hist->GetAxis(kGenYieldZ)->SetTitle("z = p_{T}^{track} / p_{T}^{jet}");
e131b05f 3246
5709c40a 3247 hist->GetAxis(kGenYieldXi)->SetTitle("#xi = ln(p_{T}^{jet} / p_{T}^{track})");
e131b05f 3248 }
3249
3250 hist->GetAxis(GetIndexOfChargeAxisGenYield())->SetTitle("Charge (e_{0})");
3251}
3252
3253
3254//________________________________________________________________________
3255void AliAnalysisTaskPID::SetUpHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const
3256{
3257 // Sets bin limits for axes which are not standard binned and the axes titles.
3258
3259 hist->SetBinEdges(kDataPt, binsPt);
3260 hist->SetBinEdges(kDataDeltaPrimeSpecies, binsDeltaPrime);
3261 hist->SetBinEdges(kDataCentrality, binsCent);
3262
3263 if (fStoreAdditionalJetInformation)
3264 hist->SetBinEdges(kDataJetPt, binsJetPt);
3265
3266 // Set axes titles
3267 hist->GetAxis(kDataMCID)->SetTitle("MC PID");
3268 hist->GetAxis(kDataMCID)->SetBinLabel(1, "e");
3269 hist->GetAxis(kDataMCID)->SetBinLabel(2, "K");
3270 hist->GetAxis(kDataMCID)->SetBinLabel(3, "#mu");
3271 hist->GetAxis(kDataMCID)->SetBinLabel(4, "#pi");
3272 hist->GetAxis(kDataMCID)->SetBinLabel(5, "p");
3273
3274 hist->GetAxis(kDataSelectSpecies)->SetTitle("Select Species");
3275 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(1, "e");
3276 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(2, "K");
3277 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(3, "#pi");
3278 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(4, "p");
3279
5709c40a 3280 hist->GetAxis(kDataPt)->SetTitle("p_{T} (GeV/c)");
e131b05f 3281
3282 hist->GetAxis(kDataDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
3283
3284 hist->GetAxis(kDataCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
3285
3286 if (fStoreAdditionalJetInformation) {
5709c40a 3287 hist->GetAxis(kDataJetPt)->SetTitle("p_{T}^{jet} (GeV/c)");
e131b05f 3288
5709c40a 3289 hist->GetAxis(kDataZ)->SetTitle("z = p_{T}^{track} / p_{T}^{jet}");
e131b05f 3290
5709c40a 3291 hist->GetAxis(kDataXi)->SetTitle("#xi = ln(p_{T}^{jet} / p_{T}^{track})");
e131b05f 3292 }
3293
3294 hist->GetAxis(GetIndexOfChargeAxisData())->SetTitle("Charge (e_{0})");
3295
77324970
CKB		 3296 hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetTitle("TOF PID Info");
3297 // Offset is (TMath::Abs(kNoTOFinfo) + 1), such that bin 1 (= first label) corresponds to kNoTOFinfo (< 0)
3298 hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kNoTOFinfo + (TMath::Abs(kNoTOFinfo) + 1), "No TOF Info");
3299 hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kNoTOFpid + (TMath::Abs(kNoTOFinfo) + 1), "No TOF PID");
3300 hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kTOFpion + (TMath::Abs(kNoTOFinfo) + 1), "#pi");
3301 hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kTOFkaon + (TMath::Abs(kNoTOFinfo) + 1), "K");
3302 hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kTOFproton + (TMath::Abs(kNoTOFinfo) + 1), "p");
e131b05f 3303}
9e95a906 3304
3305
3306//________________________________________________________________________
e4351829 3307void AliAnalysisTaskPID::SetUpPtResHist(THnSparse* hist, Double_t* binsPt, Double_t* binsJetPt, Double_t* binsCent) const
9e95a906 3308{
3309 // Sets bin limits for axes which are not standard binned and the axes titles.
3310
3311 hist->SetBinEdges(kPtResJetPt, binsJetPt);
3312 hist->SetBinEdges(kPtResGenPt, binsPt);
3313 hist->SetBinEdges(kPtResRecPt, binsPt);
e4351829 3314 hist->SetBinEdges(kPtResCentrality, binsCent);
9e95a906 3315
3316 // Set axes titles
5709c40a 3317 hist->GetAxis(kPtResJetPt)->SetTitle("p_{T}^{jet, rec} (GeV/c)");
3318 hist->GetAxis(kPtResGenPt)->SetTitle("p_{T}^{gen} (GeV/c)");
3319 hist->GetAxis(kPtResRecPt)->SetTitle("p_{T}^{rec} (GeV/c)");
e4351829 3320
3321 hist->GetAxis(kPtResCharge)->SetTitle("Charge (e_{0})");
3322 hist->GetAxis(kPtResCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
9e95a906 3323}