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