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