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