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