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