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