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