]> 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
Initial Tasks and macros for PIDed Fragmentation Function (Benjamin Hess)
[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"
17#include "AliESDInputHandler.h"
18#include "AliInputEventHandler.h"
19
20#include "AliVVertex.h"
21#include "AliAnalysisFilter.h"
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
40//________________________________________________________________________
41AliAnalysisTaskPID::AliAnalysisTaskPID()
42 : AliAnalysisTaskPIDV0base()
43 , fPIDcombined(new AliPIDCombined())
44 , fInputFromOtherTask(kFALSE)
45 , fStoreCentralityPercentile(kFALSE)
46 , fStoreAdditionalJetInformation(kFALSE)
47 , fTakeIntoAccountMuons(kFALSE)
48 , fUseITS(kFALSE)
49 , fUseTOF(kFALSE)
50 , fUsePriors(kFALSE)
51 , fTPCDefaultPriors(kFALSE)
52 , fUseMCidForGeneration(kTRUE)
53 , fUseConvolutedGaus(kFALSE)
54 , fkConvolutedGausNPar(3)
55 , fAccuracyNonGaussianTail(1e-8)
56 , fkDeltaPrimeLowLimit(0.02)
57 , fkDeltaPrimeUpLimit(40.0)
58 , fConvolutedGausDeltaPrime(0x0)
59 , fEtaAbsCutLow(0.0)
60 , fEtaAbsCutUp(0.9)
61 , fDoAnySystematicStudiesOnTheExpectedSignal(kFALSE)
62 , fSystematicScalingSplines(1.0)
63 , fSystematicScalingEtaCorrectionMomentumThr(0.35)
64 , fSystematicScalingEtaCorrectionLowMomenta(1.0)
65 , fSystematicScalingEtaCorrectionHighMomenta(1.0)
66 , fSystematicScalingEtaSigmaPara(1.0)
67 , fSystematicScalingMultCorrection(1.0)
68 , fCentralityEstimator("V0M")
69 , fhPIDdataAll(0x0)
70 , fhGenEl(0x0)
71 , fhGenKa(0x0)
72 , fhGenPi(0x0)
73 , fhGenMu(0x0)
74 , fhGenPr(0x0)
75 , fGenRespElDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
76 , fGenRespElDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
77 , fGenRespElDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
78 , fGenRespElDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
79 , fGenRespKaDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
80 , fGenRespKaDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
81 , fGenRespKaDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
82 , fGenRespKaDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
83 , fGenRespPiDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
84 , fGenRespPiDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
85 , fGenRespPiDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
86 , fGenRespPiDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
87 , fGenRespMuDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
88 , fGenRespMuDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
89 , fGenRespMuDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
90 , fGenRespMuDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
91 , fGenRespPrDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
92 , fGenRespPrDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
93 , fGenRespPrDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
94 , fGenRespPrDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
95 /*
96 , fGenRespElDeltaEl(new Double_t[fgkMaxNumGenEntries])
97 , fGenRespElDeltaKa(new Double_t[fgkMaxNumGenEntries])
98 , fGenRespElDeltaPi(new Double_t[fgkMaxNumGenEntries])
99 , fGenRespElDeltaPr(new Double_t[fgkMaxNumGenEntries])
100 , fGenRespKaDeltaEl(new Double_t[fgkMaxNumGenEntries])
101 , fGenRespKaDeltaKa(new Double_t[fgkMaxNumGenEntries])
102 , fGenRespKaDeltaPi(new Double_t[fgkMaxNumGenEntries])
103 , fGenRespKaDeltaPr(new Double_t[fgkMaxNumGenEntries])
104 , fGenRespPiDeltaEl(new Double_t[fgkMaxNumGenEntries])
105 , fGenRespPiDeltaKa(new Double_t[fgkMaxNumGenEntries])
106 , fGenRespPiDeltaPi(new Double_t[fgkMaxNumGenEntries])
107 , fGenRespPiDeltaPr(new Double_t[fgkMaxNumGenEntries])
108 , fGenRespMuDeltaEl(new Double_t[fgkMaxNumGenEntries])
109 , fGenRespMuDeltaKa(new Double_t[fgkMaxNumGenEntries])
110 , fGenRespMuDeltaPi(new Double_t[fgkMaxNumGenEntries])
111 , fGenRespMuDeltaPr(new Double_t[fgkMaxNumGenEntries])
112 , fGenRespPrDeltaEl(new Double_t[fgkMaxNumGenEntries])
113 , fGenRespPrDeltaKa(new Double_t[fgkMaxNumGenEntries])
114 , fGenRespPrDeltaPi(new Double_t[fgkMaxNumGenEntries])
115 , fGenRespPrDeltaPr(new Double_t[fgkMaxNumGenEntries])
116 */
117 , fhEventsProcessed(0x0)
118 , fhSkippedTracksForSignalGeneration(0x0)
119 , fhMCgeneratedYieldsPrimaries(0x0)
120 , fh2FFJetPtRec(0x0)
121 , fh2FFJetPtGen(0x0)
122 , fh1Xsec(0x0)
123 , fh1Trials(0x0)
124 , fContainerEff(0x0)
125 , fOutputContainer(0x0)
126{
127 // default Constructor
128
129 AliLog::SetClassDebugLevel("AliAnalysisTaskPID", AliLog::kInfo);
130
131 fConvolutedGausDeltaPrime = new TF1("convolutedGausDeltaPrime", this, &AliAnalysisTaskPID::ConvolutedGaus,
132 fkDeltaPrimeLowLimit, fkDeltaPrimeUpLimit,
133 fkConvolutedGausNPar, "AliAnalysisTaskPID", "ConvolutedGaus");
134
135 // Initialisation of translation parameters is time consuming.
136 // Therefore, default values will only be initialised if they are really needed.
137 // Otherwise, it is left to the user to set the parameter properly.
138 fConvolutedGaussTransitionPars[0] = -999;
139 fConvolutedGaussTransitionPars[1] = -999;
140 fConvolutedGaussTransitionPars[2] = -999;
141
142 // Fraction histos
143 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
144 fFractionHists[i] = 0x0;
145 fFractionSysErrorHists[i] = 0x0;
146 }
147}
148
149//________________________________________________________________________
150AliAnalysisTaskPID::AliAnalysisTaskPID(const char *name)
151 : AliAnalysisTaskPIDV0base(name)
152 , fPIDcombined(new AliPIDCombined())
153 , fInputFromOtherTask(kFALSE)
154 , fStoreCentralityPercentile(kFALSE)
155 , fStoreAdditionalJetInformation(kFALSE)
156 , fTakeIntoAccountMuons(kFALSE)
157 , fUseITS(kFALSE)
158 , fUseTOF(kFALSE)
159 , fUsePriors(kFALSE)
160 , fTPCDefaultPriors(kFALSE)
161 , fUseMCidForGeneration(kTRUE)
162 , fUseConvolutedGaus(kFALSE)
163 , fkConvolutedGausNPar(3)
164 , fAccuracyNonGaussianTail(1e-8)
165 , fkDeltaPrimeLowLimit(0.02)
166 , fkDeltaPrimeUpLimit(40.0)
167 , fConvolutedGausDeltaPrime(0x0)
168 , fEtaAbsCutLow(0.0)
169 , fEtaAbsCutUp(0.9)
170 , fDoAnySystematicStudiesOnTheExpectedSignal(kFALSE)
171 , fSystematicScalingSplines(1.0)
172 , fSystematicScalingEtaCorrectionMomentumThr(0.35)
173 , fSystematicScalingEtaCorrectionLowMomenta(1.0)
174 , fSystematicScalingEtaCorrectionHighMomenta(1.0)
175 , fSystematicScalingEtaSigmaPara(1.0)
176 , fSystematicScalingMultCorrection(1.0)
177 , fCentralityEstimator("V0M")
178 , fhPIDdataAll(0x0)
179 , fhGenEl(0x0)
180 , fhGenKa(0x0)
181 , fhGenPi(0x0)
182 , fhGenMu(0x0)
183 , fhGenPr(0x0)
184 , fGenRespElDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
185 , fGenRespElDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
186 , fGenRespElDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
187 , fGenRespElDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
188 , fGenRespKaDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
189 , fGenRespKaDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
190 , fGenRespKaDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
191 , fGenRespKaDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
192 , fGenRespPiDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
193 , fGenRespPiDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
194 , fGenRespPiDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
195 , fGenRespPiDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
196 , fGenRespMuDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
197 , fGenRespMuDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
198 , fGenRespMuDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
199 , fGenRespMuDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
200 , fGenRespPrDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
201 , fGenRespPrDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
202 , fGenRespPrDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
203 , fGenRespPrDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
204 /*
205 , fGenRespElDeltaEl(new Double_t[fgkMaxNumGenEntries])
206 , fGenRespElDeltaKa(new Double_t[fgkMaxNumGenEntries])
207 , fGenRespElDeltaPi(new Double_t[fgkMaxNumGenEntries])
208 , fGenRespElDeltaPr(new Double_t[fgkMaxNumGenEntries])
209 , fGenRespKaDeltaEl(new Double_t[fgkMaxNumGenEntries])
210 , fGenRespKaDeltaKa(new Double_t[fgkMaxNumGenEntries])
211 , fGenRespKaDeltaPi(new Double_t[fgkMaxNumGenEntries])
212 , fGenRespKaDeltaPr(new Double_t[fgkMaxNumGenEntries])
213 , fGenRespPiDeltaEl(new Double_t[fgkMaxNumGenEntries])
214 , fGenRespPiDeltaKa(new Double_t[fgkMaxNumGenEntries])
215 , fGenRespPiDeltaPi(new Double_t[fgkMaxNumGenEntries])
216 , fGenRespPiDeltaPr(new Double_t[fgkMaxNumGenEntries])
217 , fGenRespMuDeltaEl(new Double_t[fgkMaxNumGenEntries])
218 , fGenRespMuDeltaKa(new Double_t[fgkMaxNumGenEntries])
219 , fGenRespMuDeltaPi(new Double_t[fgkMaxNumGenEntries])
220 , fGenRespMuDeltaPr(new Double_t[fgkMaxNumGenEntries])
221 , fGenRespPrDeltaEl(new Double_t[fgkMaxNumGenEntries])
222 , fGenRespPrDeltaKa(new Double_t[fgkMaxNumGenEntries])
223 , fGenRespPrDeltaPi(new Double_t[fgkMaxNumGenEntries])
224 , fGenRespPrDeltaPr(new Double_t[fgkMaxNumGenEntries])
225 */
226 , fhEventsProcessed(0x0)
227 , fhSkippedTracksForSignalGeneration(0x0)
228 , fhMCgeneratedYieldsPrimaries(0x0)
229 , fh2FFJetPtRec(0x0)
230 , fh2FFJetPtGen(0x0)
231 , fh1Xsec(0x0)
232 , fh1Trials(0x0)
233 , fContainerEff(0x0)
234 , fOutputContainer(0x0)
235{
236 // Constructor
237
238 AliLog::SetClassDebugLevel("AliAnalysisTaskPID", AliLog::kInfo);
239
240 fConvolutedGausDeltaPrime = new TF1("convolutedGausDeltaPrime", this, &AliAnalysisTaskPID::ConvolutedGaus,
241 fkDeltaPrimeLowLimit, fkDeltaPrimeUpLimit,
242 fkConvolutedGausNPar, "AliAnalysisTaskPID", "ConvolutedGaus");
243
244 // Initialisation of translation parameters is time consuming.
245 // Therefore, default values will only be initialised if they are really needed.
246 // Otherwise, it is left to the user to set the parameter properly.
247 fConvolutedGaussTransitionPars[0] = -999;
248 fConvolutedGaussTransitionPars[1] = -999;
249 fConvolutedGaussTransitionPars[2] = -999;
250
251 // Fraction histos
252 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
253 fFractionHists[i] = 0x0;
254 fFractionSysErrorHists[i] = 0x0;
255 }
256
257 // Define input and output slots here
258 // Input slot #0 works with a TChain
259 DefineInput(0, TChain::Class());
260 // Output slot #0 writes into a TH1 container
261 DefineOutput(1, TObjArray::Class());
262
263 DefineOutput(2, AliCFContainer::Class());
264}
265
266
267//________________________________________________________________________
268AliAnalysisTaskPID::~AliAnalysisTaskPID()
269{
270 // dtor
271
272 CleanupParticleFractionHistos();
273
274 delete fOutputContainer;
275 fOutputContainer = 0;
276
277 delete fConvolutedGausDeltaPrime;
278 fConvolutedGausDeltaPrime = 0;
279
280 delete [] fGenRespElDeltaPrimeEl;
281 delete [] fGenRespElDeltaPrimeKa;
282 delete [] fGenRespElDeltaPrimePi;
283 delete [] fGenRespElDeltaPrimePr;
284
285 fGenRespElDeltaPrimeEl = 0x0;
286 fGenRespElDeltaPrimeKa = 0x0;
287 fGenRespElDeltaPrimePi = 0x0;
288 fGenRespElDeltaPrimePr = 0x0;
289
290 delete [] fGenRespKaDeltaPrimeEl;
291 delete [] fGenRespKaDeltaPrimeKa;
292 delete [] fGenRespKaDeltaPrimePi;
293 delete [] fGenRespKaDeltaPrimePr;
294
295 fGenRespKaDeltaPrimeEl = 0x0;
296 fGenRespKaDeltaPrimeKa = 0x0;
297 fGenRespKaDeltaPrimePi = 0x0;
298 fGenRespKaDeltaPrimePr = 0x0;
299
300 delete [] fGenRespPiDeltaPrimeEl;
301 delete [] fGenRespPiDeltaPrimeKa;
302 delete [] fGenRespPiDeltaPrimePi;
303 delete [] fGenRespPiDeltaPrimePr;
304
305 fGenRespPiDeltaPrimeEl = 0x0;
306 fGenRespPiDeltaPrimeKa = 0x0;
307 fGenRespPiDeltaPrimePi = 0x0;
308 fGenRespPiDeltaPrimePr = 0x0;
309
310 delete [] fGenRespMuDeltaPrimeEl;
311 delete [] fGenRespMuDeltaPrimeKa;
312 delete [] fGenRespMuDeltaPrimePi;
313 delete [] fGenRespMuDeltaPrimePr;
314
315 fGenRespMuDeltaPrimeEl = 0x0;
316 fGenRespMuDeltaPrimeKa = 0x0;
317 fGenRespMuDeltaPrimePi = 0x0;
318 fGenRespMuDeltaPrimePr = 0x0;
319
320 delete [] fGenRespPrDeltaPrimeEl;
321 delete [] fGenRespPrDeltaPrimeKa;
322 delete [] fGenRespPrDeltaPrimePi;
323 delete [] fGenRespPrDeltaPrimePr;
324
325 fGenRespPrDeltaPrimeEl = 0x0;
326 fGenRespPrDeltaPrimeKa = 0x0;
327 fGenRespPrDeltaPrimePi = 0x0;
328 fGenRespPrDeltaPrimePr = 0x0;
329
330 /*OLD with deltaSpecies
331 delete [] fGenRespElDeltaEl;
332 delete [] fGenRespElDeltaKa;
333 delete [] fGenRespElDeltaPi;
334 delete [] fGenRespElDeltaPr;
335
336 fGenRespElDeltaEl = 0x0;
337 fGenRespElDeltaKa = 0x0;
338 fGenRespElDeltaPi = 0x0;
339 fGenRespElDeltaPr = 0x0;
340
341 delete [] fGenRespKaDeltaEl;
342 delete [] fGenRespKaDeltaKa;
343 delete [] fGenRespKaDeltaPi;
344 delete [] fGenRespKaDeltaPr;
345
346 fGenRespKaDeltaEl = 0x0;
347 fGenRespKaDeltaKa = 0x0;
348 fGenRespKaDeltaPi = 0x0;
349 fGenRespKaDeltaPr = 0x0;
350
351 delete [] fGenRespPiDeltaEl;
352 delete [] fGenRespPiDeltaKa;
353 delete [] fGenRespPiDeltaPi;
354 delete [] fGenRespPiDeltaPr;
355
356 fGenRespPiDeltaEl = 0x0;
357 fGenRespPiDeltaKa = 0x0;
358 fGenRespPiDeltaPi = 0x0;
359 fGenRespPiDeltaPr = 0x0;
360
361 delete [] fGenRespMuDeltaEl;
362 delete [] fGenRespMuDeltaKa;
363 delete [] fGenRespMuDeltaPi;
364 delete [] fGenRespMuDeltaPr;
365
366 fGenRespMuDeltaEl = 0x0;
367 fGenRespMuDeltaKa = 0x0;
368 fGenRespMuDeltaPi = 0x0;
369 fGenRespMuDeltaPr = 0x0;
370
371 delete [] fGenRespPrDeltaEl;
372 delete [] fGenRespPrDeltaKa;
373 delete [] fGenRespPrDeltaPi;
374 delete [] fGenRespPrDeltaPr;
375
376 fGenRespPrDeltaEl = 0x0;
377 fGenRespPrDeltaKa = 0x0;
378 fGenRespPrDeltaPi = 0x0;
379 fGenRespPrDeltaPr = 0x0;
380 */
381}
382
383
384//________________________________________________________________________
385void AliAnalysisTaskPID::SetUpPIDcombined()
386{
387 // Initialise the PIDcombined object
388
389 if (!fPIDcombined) {
390 AliFatal("No PIDcombined object!\n");
391 return;
392 }
393
394 fPIDcombined->SetSelectedSpecies(AliPID::kSPECIESC);
395 fPIDcombined->SetEnablePriors(fUsePriors);
396
397 if (fTPCDefaultPriors)
398 fPIDcombined->SetDefaultTPCPriors();
399
400 //TODO use individual priors...
401
402 // Change detector mask (TPC,TOF,ITS)
403 Int_t detectorMask = AliPIDResponse::kDetTPC;
404
405 // Reject mismatch mask - mismatch only relevant for TOF at the moment - other detectors do not use it
406 Int_t rejectMismatchMask = AliPIDResponse::kDetTPC;
407
408
409 if (fUseITS) {
410 detectorMask = detectorMask | AliPIDResponse::kDetITS;
411 rejectMismatchMask = rejectMismatchMask | AliPIDResponse::kDetITS;
412 }
413 if (fUseTOF) {
414 detectorMask = detectorMask | AliPIDResponse::kDetTOF;
415 rejectMismatchMask = rejectMismatchMask | AliPIDResponse::kDetTOF;
416 }
417
418 fPIDcombined->SetDetectorMask(detectorMask);
419 fPIDcombined->SetRejectMismatchMask(rejectMismatchMask);
420}
421
422
423//________________________________________________________________________
424void AliAnalysisTaskPID::UserCreateOutputObjects()
425{
426 // Create histograms
427 // Called once
428
429 SetUpPIDcombined();
430
431 // Input handler
432 AliAnalysisManager* man = AliAnalysisManager::GetAnalysisManager();
433 AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
434
435 if (!inputHandler)
436 AliFatal("Input handler needed");
437 else {
438 // PID response object
439 fPIDResponse = inputHandler->GetPIDResponse();
440 if (!fPIDResponse)
441 AliFatal("PIDResponse object was not created");
442 }
443
444 OpenFile(1);
445
446 fOutputContainer = new TObjArray(1);
447 fOutputContainer->SetName(GetName()) ;
448 fOutputContainer->SetOwner(kTRUE);
449
450 const Int_t nPtBins = 68;
451 Double_t binsPt[nPtBins+1] = {0. , 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45,
452 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95,
453 1.0, 1.1 , 1.2, 1.3 , 1.4, 1.5 , 1.6, 1.7 , 1.8, 1.9 ,
454 2.0, 2.2 , 2.4, 2.6 , 2.8, 3.0 , 3.2, 3.4 , 3.6, 3.8 ,
455 4.0, 4.5 , 5.0, 5.5 , 6.0, 6.5 , 7.0, 8.0 , 9.0, 10.0,
456 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0,
457 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 };
458
459 const Int_t nCentBins = 12;
460 //-1 for pp; 90-100 has huge electro-magnetic impurities
461 Double_t binsCent[nCentBins+1] = {-1, 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 };
462
463 const Int_t nJetPtBins = 11;
464 Double_t binsJetPt[nJetPtBins+1] = {0, 2, 5, 10, 15, 20, 30, 40, 60, 80, 120, 200};
465
466 const Int_t nChargeBins = 2;
467 const Double_t binsCharge[nChargeBins+1] = { -1.0 - 1e-4, 0.0, 1.0 + 1e-4 };
468
469 const Int_t nBinsJets = kDataNumAxes;
470 const Int_t nBinsNoJets = nBinsJets - fgkNumJetAxes;
471
472 const Int_t nBins = fStoreAdditionalJetInformation ? nBinsJets : nBinsNoJets;
473
474 // deltaPrimeSpecies binning
475 const Int_t deltaPrimeNBins = 600;
476 Double_t deltaPrimeBins[deltaPrimeNBins + 1];
477
478 const Double_t fromLow = fkDeltaPrimeLowLimit;
479 const Double_t toHigh = fkDeltaPrimeUpLimit;
480 const Double_t factor = TMath::Power(toHigh/fromLow, 1./deltaPrimeNBins);
481
482 // Log binning for whole deltaPrime range
483 deltaPrimeBins[0] = fromLow;
484 for (Int_t i = 0 + 1; i <= deltaPrimeNBins; i++) {
485 deltaPrimeBins[i] = factor * deltaPrimeBins[i - 1];
486 }
487
488 const Int_t nMCPIDbins = 5;
489 const Double_t mcPIDmin = 0.;
490 const Double_t mcPIDmax = 5.;
491
492 const Int_t nSelSpeciesBins = 4;
493 const Double_t selSpeciesMin = 0.;
494 const Double_t selSpeciesMax = 4.;
495
496 const Int_t nZBins = 20;
497 const Double_t zMin = 0.;
498 const Double_t zMax = 1.;
499
500 const Int_t nXiBins = 70;
501 const Double_t xiMin = 0.;
502 const Double_t xiMax = 7.;
503
504 // MC PID, SelectSpecies, pT, deltaPrimeSpecies, centrality percentile, jet pT, z = track_pT/jet_pT, xi = log(1/z)
505 Int_t binsNoJets[nBinsNoJets] = { nMCPIDbins, nSelSpeciesBins, nPtBins, deltaPrimeNBins,
506 nCentBins, nChargeBins};
507 Int_t binsJets[nBinsJets] = { nMCPIDbins, nSelSpeciesBins, nPtBins, deltaPrimeNBins,
508 nCentBins, nJetPtBins, nZBins, nXiBins, nChargeBins };
509 Int_t *bins = fStoreAdditionalJetInformation ? &binsJets[0] : &binsNoJets[0];
510
511 Double_t xminNoJets[nBinsNoJets] = { mcPIDmin, selSpeciesMin, binsPt[0], deltaPrimeBins[0],
512 binsCent[0], binsCharge[0]};
513 Double_t xminJets[nBinsJets] = { mcPIDmin, selSpeciesMin, binsPt[0], deltaPrimeBins[0],
514 binsCent[0], binsJetPt[0], zMin, xiMin, binsCharge[0] };
515 Double_t *xmin = fStoreAdditionalJetInformation? &xminJets[0] : &xminNoJets[0];
516
517 Double_t xmaxNoJets[nBinsNoJets] = { mcPIDmax, selSpeciesMax, binsPt[nPtBins], deltaPrimeBins[deltaPrimeNBins],
518 binsCent[nCentBins], binsCharge[nChargeBins] };
519 Double_t xmaxJets[nBinsJets] = { mcPIDmax, selSpeciesMax, binsPt[nPtBins], deltaPrimeBins[deltaPrimeNBins],
520 binsCent[nCentBins], binsJetPt[nJetPtBins], zMax, xiMax, binsCharge[nChargeBins] };
521 Double_t *xmax = fStoreAdditionalJetInformation? &xmaxJets[0] : &xmaxNoJets[0];
522
523 /*OLD with TOF, p_TPC_Inner and p_vertex and deltaSpecies
524 const Int_t nBins = 8;
525 //TODO In case of memory trouble: Remove deltaTOFspecies and p(Vertex) (can be added later, if needed)?
526 //TODO IF everything is working fine: p(TPC_inner) and p(Vertex) can be removed, since everything in the analysis is only a
527 // function of pT -> Reduces memory consumption significantly!!!
528
529 // MC PID, SelectSpecies, P(TPC_inner), pT, p(Vertex), deltaSpecies, deltaPrimeSpecies, deltaTOFspecies
530 const Int_t deltaPrimeNBins = 201;
531
532 const Int_t deltaNBins = 801;
533 const Double_t deltaLowLimit = -200.;
534 const Double_t deltaUpLimit = 200.;
535
536 Int_t bins[nBins] =
537 { 5, 4, nPtBins, nPtBins, nPtBins, deltaNBins, deltaPrimeNBins, 501 };
538 Double_t xmin[nBins] =
539 { 0., 0., 0., 0., 0., deltaLowLimit, fkDeltaPrimeLowLimit, -5000.};
540 Double_t xmax[nBins] =
541 { 5., 4., 50.0, 50.0, 50.0, deltaUpLimit, fkDeltaPrimeUpLimit, 5000.};
542 */
543
544 fConvolutedGausDeltaPrime->SetNpx(deltaPrimeNBins);
545
546 fhPIDdataAll = new THnSparseD("hPIDdataAll","", nBins, bins, xmin, xmax);
547 SetUpHist(fhPIDdataAll, binsPt, deltaPrimeBins, binsCent, binsJetPt);
548 fOutputContainer->Add(fhPIDdataAll);
549
550 // Generated histograms (so far, bins are the same as for primary THnSparse)
551 const Int_t nGenBins = fStoreAdditionalJetInformation ? nBinsJets : nBinsNoJets;
552 // MC PID, SelectSpecies, Pt, deltaPrimeSpecies, jet pT, z = track_pT/jet_pT, xi = log(1/z)
553
554 Int_t *genBins = fStoreAdditionalJetInformation ? &binsJets[0] : &binsNoJets[0];
555 Double_t *genXmin = fStoreAdditionalJetInformation? &xminJets[0] : &xminNoJets[0];
556 Double_t *genXmax = fStoreAdditionalJetInformation? &xmaxJets[0] : &xmaxNoJets[0];
557
558 fhGenEl = new THnSparseD("hGenEl", "", nGenBins, genBins, genXmin, genXmax);
559 SetUpGenHist(fhGenEl, binsPt, deltaPrimeBins, binsCent, binsJetPt);
560 fOutputContainer->Add(fhGenEl);
561
562 fhGenKa = new THnSparseD("hGenKa", "", nGenBins, genBins, genXmin, genXmax);
563 SetUpGenHist(fhGenKa, binsPt, deltaPrimeBins, binsCent, binsJetPt);
564 fOutputContainer->Add(fhGenKa);
565
566 fhGenPi = new THnSparseD("hGenPi", "", nGenBins, genBins, genXmin, genXmax);
567 SetUpGenHist(fhGenPi, binsPt, deltaPrimeBins, binsCent, binsJetPt);
568 fOutputContainer->Add(fhGenPi);
569
570 if (fTakeIntoAccountMuons) {
571 fhGenMu = new THnSparseD("hGenMu", "", nGenBins, genBins, genXmin, genXmax);
572 SetUpGenHist(fhGenMu, binsPt, deltaPrimeBins, binsCent, binsJetPt);
573 fOutputContainer->Add(fhGenMu);
574 }
575
576 fhGenPr = new THnSparseD("hGenPr", "", nGenBins, genBins, genXmin, genXmax);
577 SetUpGenHist(fhGenPr, binsPt, deltaPrimeBins, binsCent, binsJetPt);
578 fOutputContainer->Add(fhGenPr);
579
580
581 fhEventsProcessed = new TH1D("fhEventsProcessed","Number of processed events;Centrality percentile", nCentBins, binsCent);
582 fhEventsProcessed->Sumw2();
583 fOutputContainer->Add(fhEventsProcessed);
584
585 fhSkippedTracksForSignalGeneration = new TH2D("fhSkippedTracksForSignalGeneration",
586 "Number of tracks skipped for the signal generation;P_{T}^{gen} (GeV/c);TPC signal N",
587 nPtBins, binsPt, 161, -0.5, 160.5);
588 fhSkippedTracksForSignalGeneration->Sumw2();
589 fOutputContainer->Add(fhSkippedTracksForSignalGeneration);
590
591
592 // Generated yields within acceptance
593 const Int_t nBinsGenYields = fStoreAdditionalJetInformation ? kGenYieldNumAxes : kGenYieldNumAxes - 3;
594 Int_t genYieldsBins[kGenYieldNumAxes] = { nMCPIDbins, nPtBins, nCentBins, nJetPtBins, nZBins, nXiBins,
595 nChargeBins };
596 genYieldsBins[GetIndexOfChargeAxisGenYield()] = nChargeBins;
597 Double_t genYieldsXmin[kGenYieldNumAxes] = { mcPIDmin, binsPt[0], binsCent[0], binsJetPt[0], zMin, xiMin,
598 binsCharge[0] };
599 genYieldsXmin[GetIndexOfChargeAxisGenYield()] = binsCharge[0];
600 Double_t genYieldsXmax[kGenYieldNumAxes] = { mcPIDmax, binsPt[nPtBins], binsCent[nCentBins], binsJetPt[nJetPtBins], zMax, xiMax,
601 binsCharge[nChargeBins] };
602 genYieldsXmax[GetIndexOfChargeAxisGenYield()] = binsCharge[nChargeBins];
603
604 fhMCgeneratedYieldsPrimaries = new THnSparseD("fhMCgeneratedYieldsPrimaries",
605 "Generated yields w/o reco and cuts inside acceptance (physical primaries)",
606 nBinsGenYields, genYieldsBins, genYieldsXmin, genYieldsXmax);
607 SetUpGenYieldHist(fhMCgeneratedYieldsPrimaries, binsPt, binsCent, binsJetPt);
608 fOutputContainer->Add(fhMCgeneratedYieldsPrimaries);
609
610 // Container with several process steps (generated and reconstructed level with some variations)
611 OpenFile(2);
612
613 // Array for the number of bins in each dimension
614 // Dimensions: MC-ID, trackPt, trackEta, trackCharge, cenrality percentile, jetPt, z, xi TODO phi???
615 const Int_t nEffDims = fStoreAdditionalJetInformation ? kEffNumAxes : kEffNumAxes - 3; // Number of dimensions for the efficiency
616
617 const Int_t nMCIDbins = AliPID::kSPECIES;
618 Double_t binsMCID[nMCIDbins];
619
620 for(Int_t i = 0; i <= nMCIDbins; i++) {
621 binsMCID[i]= i;
622 }
623
624 const Int_t nEtaBins = 18;
625 const Double_t binsEta[nEtaBins+1] = {-0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1,
626 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 };
627
628 const Int_t nEffBins[kEffNumAxes] = { nMCIDbins, nPtBins, nEtaBins, nChargeBins, nCentBins, nJetPtBins, nZBins, nXiBins };
629
630 fContainerEff = new AliCFContainer("containerEff", "Reconstruction Efficiency x Acceptance x Resolution and Secondary Correction",
631 kNumSteps, nEffDims, nEffBins);
632
633 // Setting the bin limits
634 fContainerEff->SetBinLimits(kEffMCID, binsMCID);
635 fContainerEff->SetBinLimits(kEffTrackPt, binsPt);
636 fContainerEff->SetBinLimits(kEffTrackEta, binsEta);
637 fContainerEff->SetBinLimits(kEffTrackCharge, binsCharge);
638 fContainerEff->SetBinLimits(kEffCentrality, binsCent);
639 if (fStoreAdditionalJetInformation) {
640 fContainerEff->SetBinLimits(kEffJetPt, binsJetPt);
641 fContainerEff->SetBinLimits(kEffZ, zMin, zMax);
642 fContainerEff->SetBinLimits(kEffXi, xiMin, xiMax);
643 }
644
645 fContainerEff->SetVarTitle(kEffMCID,"MC ID");
646 fContainerEff->SetVarTitle(kEffTrackPt,"P_{T} (GeV/c)");
647 fContainerEff->SetVarTitle(kEffTrackEta,"#eta");
648 fContainerEff->SetVarTitle(kEffTrackCharge,"Charge (e_{0})");
649 fContainerEff->SetVarTitle(kEffCentrality, "Centrality Percentile");
650 if (fStoreAdditionalJetInformation) {
651 fContainerEff->SetVarTitle(kEffJetPt, "P_{T}^{jet} (GeV/c)");
652 fContainerEff->SetVarTitle(kEffZ, "z = P_{T}^{track} / P_{T}^{jet}");
653 fContainerEff->SetVarTitle(kEffXi, "#xi = ln(P_{T}^{jet} / P_{T}^{track})");
654 }
655
656 // Define clean MC sample
657 fContainerEff->SetStepTitle(kStepGenWithGenCuts, "Particle level, cuts on particle level");
658 // For Acceptance x Efficiency correction of primaries
659 fContainerEff->SetStepTitle(kStepRecWithGenCuts, "Detector level (rec) with cuts on particle level");
660 // For (pT) resolution correction
661 fContainerEff->SetStepTitle(kStepRecWithGenCutsMeasuredObs,
662 "Detector level (rec) with cuts on particle level with measured observables");
663 // For secondary correction
664 fContainerEff->SetStepTitle(kStepRecWithRecCutsMeasuredObs, "Detector level, all cuts on detector level with measured observables");
665 fContainerEff->SetStepTitle(kStepRecWithRecCutsMeasuredObsPrimaries,
666 "Detector level, all cuts on detector level with measured observables, only MC primaries");
667 fContainerEff->SetStepTitle(kStepRecWithRecCutsMeasuredObsStrangenessScaled,
668 "Detector level (strangeness scaled), all cuts on detector level with measured observables");
669
670 // Generated jets
671 fh2FFJetPtRec = new TH2D("fh2FFJetPtRec", "Number of reconstructed jets;Centrality Percentile;P_{T}^{jet} (GeV/c)",
672 nCentBins, binsCent, nJetPtBins, binsJetPt);
673 fh2FFJetPtRec->Sumw2();
674 fOutputContainer->Add(fh2FFJetPtRec);
675 fh2FFJetPtGen = new TH2D("fh2FFJetPtGen", "Number of generated jets;Centrality Percentile;P_{T}^{jet} (GeV/c)",
676 nCentBins, binsCent, nJetPtBins, binsJetPt);
677 fh2FFJetPtGen->Sumw2();
678 fOutputContainer->Add(fh2FFJetPtGen);
679
680 // Pythia information
681 fh1Xsec = new TProfile("fh1Xsec", "xsec from pyxsec.root", 1, 0, 1);
682 fh1Xsec->Sumw2();
683 fh1Xsec->GetXaxis()->SetBinLabel(1, "<#sigma>");
684 fh1Trials = new TH1D("fh1Trials", "trials from pyxsec.root", 1, 0, 1);
685 fh1Trials->Sumw2();
686 fh1Trials->GetXaxis()->SetBinLabel(1, "#sum{ntrials}");
687
688 fOutputContainer->Add(fh1Xsec);
689 fOutputContainer->Add(fh1Trials);
690
691 PostOutputData();
692}
693
694
695//________________________________________________________________________
696void AliAnalysisTaskPID::UserExec(Option_t *)
697{
698 // Main loop
699 // Called for each event
700
701 // No processing of event, if input is fed in directly from another task
702 if (fInputFromOtherTask)
703 return;
704
705 fEvent = dynamic_cast<AliVEvent*>(InputEvent());
706 if (!fEvent) {
707 Printf("ERROR: fEvent not available");
708 return;
709 }
710
711 fMC = dynamic_cast<AliMCEvent*>(MCEvent());
712
713 if (!fPIDResponse || !fPIDcombined)
714 return;
715
716 if (!GetVertexIsOk(fEvent))
717 return;
718
719 fESD = dynamic_cast<AliESDEvent*>(fEvent);
720 const AliVVertex* primaryVertex = fESD ? fESD->GetPrimaryVertexTracks() : fEvent->GetPrimaryVertex();
721 if (!primaryVertex)
722 return;
723
724 if(primaryVertex->GetNContributors() <= 0)
725 return;
726
727 Double_t magField = fEvent->GetMagneticField();
728
729 //OLD with DeltaSpecies const Bool_t usePureGausForDelta = kTRUE;
730
731
732 Double_t centralityPercentile = -1;
733 if (fStoreCentralityPercentile)
734 centralityPercentile = fEvent->GetCentrality()->GetCentralityPercentile(fCentralityEstimator.Data());
735
736 if (fMC) {
737 for (Int_t iPart = 0; iPart < fMC->GetNumberOfTracks(); iPart++) {
738 AliMCParticle *mcPart = dynamic_cast<AliMCParticle*>(fMC->GetTrack(iPart));
739
740 if (!mcPart)
741 continue;
742
743 // Define clean MC sample with corresponding particle level track cuts:
744 // - MC-track must be in desired eta range
745 // - MC-track must be physical primary
746 // - Species must be one of those in question (everything else goes to the overflow bin of mcID)
747
748 // Geometrie should be the same as on the reconstructed level -> By definition analysis within this eta interval
749 if (TMath::Abs(mcPart->Eta()) < fEtaAbsCutLow || TMath::Abs(mcPart->Eta()) > fEtaAbsCutUp) continue;
750
751 Int_t mcID = PDGtoMCID(mcPart->PdgCode());
752
753 // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
754 Double_t chargeMC = mcPart->Charge() / 3.;
755
756 if (TMath::Abs(chargeMC) < 0.01)
757 continue; // Reject neutral particles (only relevant, if mcID is not used)
758
759 Double_t valuesGenYield[kGenYieldNumAxes] = { mcID, mcPart->Pt(), centralityPercentile, -1, -1, -1, -1 };
760 valuesGenYield[GetIndexOfChargeAxisGenYield()] = chargeMC;
761
762 if (!fMC->IsPhysicalPrimary(iPart))
763 continue;
764
765 fhMCgeneratedYieldsPrimaries->Fill(valuesGenYield);
766
767
768 Double_t valueEff[kEffNumAxes] = { mcID, mcPart->Pt(), mcPart->Eta(), chargeMC, centralityPercentile,
769 -1, -1, -1 };
770
771 fContainerEff->Fill(valueEff, kStepGenWithGenCuts);
772 }
773 }
774
775 // Track loop to fill a Train spectrum
776 for (Int_t iTracks = 0; iTracks < fEvent->GetNumberOfTracks(); iTracks++) {
777 AliVTrack* track = dynamic_cast<AliVTrack*>(fEvent->GetTrack(iTracks));
778 if (!track) {
779 Printf("ERROR: Could not retrieve track %d", iTracks);
780 continue;
781 }
782
783
784 // Apply detector level track cuts
785 //TODO NOW if (track->GetTPCsignalN() < 60)
786 //TODO NOW continue;//TODO was removed for a while
787
788
789 if(fTrackFilter && !fTrackFilter->IsSelected(track))
790 continue;
791
792 /*
793 if (fESD) {
794 AliESDtrack* trackESD = fESD->GetTrack(iTracks);
795 if (trackESD) {
796 if (trackESD->GetLengthInActiveZone(1, 1.8, 220, magField) < 120)
797 continue; //TODO NOW TEST
798 }
799 }*/
800
801 if(fUsePhiCut) {
802 if (!PhiPrimeCut(track, magField))
803 continue; // reject track
804 }
805
806 Double_t dEdxTPC = fPIDResponse->IsTunedOnData() ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
807
808 if (dEdxTPC <= 0)
809 continue;
810
811 Int_t pdg = 0; // = 0 indicates data for the moment
812 AliMCParticle* mcTrack = 0x0;
813 Int_t mcID = AliPID::kUnknown;
814 Int_t label = 0;
815
816 if (fMC) {
817 label = track->GetLabel();
818
819 //if (label < 0)
820 // continue;
821
822 mcTrack = dynamic_cast<AliMCParticle*>(fMC->GetTrack(TMath::Abs(label)));
823 if (!mcTrack) {
824 Printf("ERROR: Could not retrieve mcTrack with label %d for track %d", label, iTracks);
825 continue;
826 }
827
828 pdg = mcTrack->PdgCode();
829 mcID = PDGtoMCID(mcTrack->PdgCode());
830
831 // For efficiency: Reconstructed track has survived all cuts on the detector level (excluding acceptance)
832 // and has an associated MC track which is a physical primary and was generated inside the acceptance
833 if (fMC->IsPhysicalPrimary(TMath::Abs(label)) &&
834 (TMath::Abs(mcTrack->Eta()) >= fEtaAbsCutLow && TMath::Abs(mcTrack->Eta()) <= fEtaAbsCutUp)) {
835
836 // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
837 Double_t value[kEffNumAxes] = { mcID, mcTrack->Pt(), mcTrack->Eta(), mcTrack->Charge() / 3., centralityPercentile,
838 -1, -1, -1 };
839 fContainerEff->Fill(value, kStepRecWithGenCuts);
840
841 Double_t valueMeas[kEffNumAxes] = { mcID, track->Pt(), track->Eta(), track->Charge(), centralityPercentile,
842 -1, -1, -1 };
843 fContainerEff->Fill(valueMeas, kStepRecWithGenCutsMeasuredObs);
844 }
845 }
846
847 // Only process tracks inside the desired eta window
848 if (TMath::Abs(track->Eta()) < fEtaAbsCutLow || TMath::Abs(track->Eta()) > fEtaAbsCutUp) continue;
849
850 ProcessTrack(track, pdg, centralityPercentile, -1); // No jet information in this case -> Set jet pT to -1
851
852 if (mcTrack) {
853 Double_t valueRecAllCuts[kEffNumAxes] = { mcID, track->Pt(), track->Eta(), track->Charge(), centralityPercentile,
854 -1, -1, -1 };
855 fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObs);
856
857 Double_t weight = IsSecondaryWithStrangeMotherMC(fMC, TMath::Abs(label)) ? GetMCStrangenessFactorCMS(fMC, mcTrack) : 1.0;
858 fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObsStrangenessScaled, weight);
859
860 if (fMC->IsPhysicalPrimary(TMath::Abs(label)))
861 fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObsPrimaries);
862 }
863 } //track loop
864
865 IncrementEventsProcessed(centralityPercentile);
866
867 PostOutputData();
868}
869
870//________________________________________________________________________
871void AliAnalysisTaskPID::Terminate(const Option_t *)
872{
873 // Draw result to the screen
874 // Called once at the end of the query
875}
876
877
878//_____________________________________________________________________________
879void AliAnalysisTaskPID::CheckDoAnyStematicStudiesOnTheExpectedSignal()
880{
881 // Check whether at least one scale factor indicates the ussage of systematic studies
882 // and set internal flag accordingly.
883
884 fDoAnySystematicStudiesOnTheExpectedSignal = kFALSE;
885
886 if (TMath::Abs(fSystematicScalingSplines - 1.0) > fgkEpsilon) {
887 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
888 return;
889 }
890
891 if ((TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon) ||
892 (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon)) {
893 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
894 return;
895 }
896
897 if (TMath::Abs(fSystematicScalingEtaSigmaPara - 1.0) > fgkEpsilon) {
898 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
899 return;
900 }
901
902 if (TMath::Abs(fSystematicScalingMultCorrection - 1.0) > fgkEpsilon) {
903 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
904 return;
905 }
906}
907
908
909//_____________________________________________________________________________
910Int_t AliAnalysisTaskPID::PDGtoMCID(Int_t pdg)
911{
912 // Returns the corresponding AliPID index to the given pdg code.
913 // Returns AliPID::kUnkown if pdg belongs to a not considered species.
914
915 Int_t absPDGcode = TMath::Abs(pdg);
916 if (absPDGcode == 211) {//Pion
917 return AliPID::kPion;
918 }
919 else if (absPDGcode == 321) {//Kaon
920 return AliPID::kKaon;
921 }
922 else if (absPDGcode == 2212) {//Proton
923 return AliPID::kProton;
924 }
925 else if (absPDGcode == 11) {//Electron
926 return AliPID::kElectron;
927 }
928 else if (absPDGcode == 13) {//Muon
929 return AliPID::kMuon;
930 }
931
932 return AliPID::kUnknown;
933}
934
935
936//_____________________________________________________________________________
937void AliAnalysisTaskPID::GetJetTrackObservables(const Double_t trackPt, const Double_t jetPt, Double_t& z, Double_t& xi)
938{
939 // Uses trackPt and jetPt to obtain z and xi.
940
941 z = (jetPt > 0 && trackPt >= 0) ? (trackPt / jetPt) : -1;
942 xi = (z > 0) ? TMath::Log(1. / z) : -1;
943
944 if(trackPt > (1. - 1e-06) * jetPt && trackPt < (1. + 1e-06) * jetPt) { // case z=1 : move entry to last histo bin <1
945 z = 1. - 1e-06;
946 xi = 1e-06;
947 }
948}
949
950
951//_____________________________________________________________________________
952void AliAnalysisTaskPID::CleanupParticleFractionHistos()
953{
954 // Delete histos with particle fractions
955
956 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
957 delete fFractionHists[i];
958 fFractionHists[i] = 0x0;
959
960 delete fFractionSysErrorHists[i];
961 fFractionSysErrorHists[i] = 0x0;
962 }
963}
964
965
966//_____________________________________________________________________________
967Double_t AliAnalysisTaskPID::ConvolutedGaus(const Double_t* xx, const Double_t* par) const
968{
969 // Convolutes gauss with an exponential tail which describes dEdx-response better than pure gaussian
970
971 const Double_t mean = par[0];
972 const Double_t sigma = par[1];
973 const Double_t lambda = par[2];
974
975 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);
976}
977
978
979//_____________________________________________________________________________
980inline Double_t AliAnalysisTaskPID::FastGaus(const Double_t x, const Double_t mean, const Double_t sigma) const
981{
982 // Calculate an unnormalised gaussian function with mean and sigma.
983
984 if (sigma < fgkEpsilon)
985 return 1.e30;
986
987 const Double_t arg = (x - mean) / sigma;
988 return exp(-0.5 * arg * arg);
989}
990
991
992//_____________________________________________________________________________
993inline Double_t AliAnalysisTaskPID::FastNormalisedGaus(const Double_t x, const Double_t mean, const Double_t sigma) const
994{
995 // Calculate a normalised (divided by sqrt(2*Pi)*sigma) gaussian function with mean and sigma.
996
997 if (sigma < fgkEpsilon)
998 return 1.e30;
999
1000 const Double_t arg = (x - mean) / sigma;
1001 const Double_t res = exp(-0.5 * arg * arg);
1002 return res / (2.50662827463100024 * sigma); //sqrt(2*Pi)=2.50662827463100024
1003}
1004
1005
1006//_____________________________________________________________________________
1007Int_t AliAnalysisTaskPID::FindBinWithinRange(TAxis* axis, Double_t value) const
1008{
1009 // Find the corresponding bin of the axis. Values outside the range (also under and overflow) will be set to the first/last
1010 // available bin
1011
1012 Int_t bin = axis->FindFixBin(value);
1013
1014 if (bin <= 0)
1015 bin = 1;
1016 if (bin > axis->GetNbins())
1017 bin = axis->GetNbins();
1018
1019 return bin;
1020}
1021
1022
1023//_____________________________________________________________________________
1024Int_t AliAnalysisTaskPID::FindFirstBinAboveIn3dSubset(const TH3* hist, const Double_t threshold, const Double_t yBin, const Double_t zBin) const
1025{
1026 // Kind of projects a TH3 to 1 bin combination in y and z
1027 // and looks for the first x bin above a threshold for this projection.
1028 // If no such bin is found, -1 is returned.
1029
1030 if (!hist)
1031 return -1;
1032
1033 Int_t nBinsX = hist->GetNbinsX();
1034 for (Int_t xBin = 1; xBin <= nBinsX; xBin++) {
1035 if (hist->GetBinContent(xBin, yBin, zBin) > threshold)
1036 return xBin;
1037 }
1038
1039 return -1;
1040}
1041
1042
1043//_____________________________________________________________________________
1044Int_t AliAnalysisTaskPID::FindLastBinAboveIn3dSubset(const TH3* hist, const Double_t threshold, const Double_t yBin, const Double_t zBin) const
1045{
1046 // Kind of projects a TH3 to 1 bin combination in y and z
1047 // and looks for the last x bin above a threshold for this projection.
1048 // If no such bin is found, -1 is returned.
1049
1050 if (!hist)
1051 return -1;
1052
1053 Int_t nBinsX = hist->GetNbinsX();
1054 for (Int_t xBin = nBinsX; xBin >= 1; xBin--) {
1055 if (hist->GetBinContent(xBin, yBin, zBin) > threshold)
1056 return xBin;
1057 }
1058
1059 return -1;
1060}
1061
1062
1063//_____________________________________________________________________________
1064Bool_t AliAnalysisTaskPID::GetParticleFraction(const Double_t trackPt, const Double_t jetPt, const Double_t centralityPercentile,
1065 const AliPID::EParticleType species,
1066 Double_t& fraction, Double_t& fractionErrorStat, Double_t& fractionErrorSys) const
1067{
1068 // Computes the particle fraction for the corresponding species for the given trackPt, jetPt and centrality.
1069 // Use jetPt = -1 for inclusive spectra and centralityPercentile = -1 for pp.
1070 // On success (return value kTRUE), fraction contains the particle fraction, fractionErrorStat(Sys) the sigma of its
1071 // statistical (systematic) error
1072
1073 fraction = -999.;
1074 fractionErrorStat = 999.;
1075 fractionErrorSys = 999.;
1076
1077 if (species > AliPID::kProton || species < AliPID::kElectron) {
1078 AliError(Form("Only fractions for species index %d to %d availabe, but not for the requested one: %d", 0, AliPID::kProton, species));
1079 return kFALSE;
1080 }
1081
1082 if (!fFractionHists[species]) {
1083 AliError(Form("Histo with particle fractions for species %d not loaded!", species));
1084
1085 return kFALSE;
1086 }
1087
1088 Int_t jetPtBin = FindBinWithinRange(fFractionHists[species]->GetYaxis(), jetPt);
1089 Int_t centBin = FindBinWithinRange(fFractionHists[species]->GetZaxis(), centralityPercentile);
1090
1091 // The following interpolation takes the bin content of the first/last available bin,
1092 // if requested point lies beyond bin center of first/last bin.
1093 // The interpolation is only done for the x-axis (track pT), i.e. jetPtBin and centBin are fix,
1094 // because the analysis will anyhow be run in bins of jetPt and centrality and
1095 // it is not desired to correlate different jetPt bins via interpolation.
1096
1097 // The same procedure is used for the error of the fraction
1098 TAxis* xAxis = fFractionHists[species]->GetXaxis();
1099
1100 // No interpolation to values beyond the centers of the first/last bins (we don't know exactly where the spectra start or stop,
1101 // thus, search for the first and last bin above 0.0 to constrain the range
1102 Int_t firstBin = TMath::Max(1, FindFirstBinAboveIn3dSubset(fFractionHists[species], 0.0, jetPtBin, centBin));
1103 Int_t lastBin = TMath::Min(fFractionHists[species]->GetNbinsX(),
1104 FindLastBinAboveIn3dSubset(fFractionHists[species], 0.0, jetPtBin, centBin));
1105
1106 if (trackPt <= xAxis->GetBinCenter(firstBin)) {
1107 fraction = fFractionHists[species]->GetBinContent(firstBin, jetPtBin, centBin);
1108 fractionErrorStat = fFractionHists[species]->GetBinError(firstBin, jetPtBin, centBin);
1109 fractionErrorSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(firstBin, jetPtBin, centBin) : 0.;
1110 }
1111 else if (trackPt >= xAxis->GetBinCenter(lastBin)) {
1112 fraction = fFractionHists[species]->GetBinContent(lastBin, jetPtBin, centBin);
1113 fractionErrorStat = fFractionHists[species]->GetBinError(lastBin, jetPtBin, centBin);
1114 fractionErrorSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(lastBin, jetPtBin, centBin) : 0.;
1115 }
1116 else {
1117 Double_t x0 = 0., x1 = 0., y0 = 0., y1 = 0.;
1118 Double_t y0errStat = 0., y1errStat = 0., y0errSys = 0., y1errSys = 0.;
1119 Int_t trackPtBin = xAxis->FindBin(trackPt);
1120
1121 // Linear interpolation between nearest neighbours in trackPt
1122 if (trackPt <= xAxis->GetBinCenter(trackPtBin)) {
1123 y0 = fFractionHists[species]->GetBinContent(trackPtBin - 1, jetPtBin, centBin);
1124 y0errStat = fFractionHists[species]->GetBinError(trackPtBin - 1, jetPtBin, centBin);
1125 y0errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin - 1, jetPtBin, centBin)
1126 : 0.;
1127 x0 = xAxis->GetBinCenter(trackPtBin - 1);
1128 y1 = fFractionHists[species]->GetBinContent(trackPtBin, jetPtBin, centBin);
1129 y1errStat = fFractionHists[species]->GetBinError(trackPtBin, jetPtBin, centBin);
1130 y1errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin, jetPtBin, centBin)
1131 : 0.;
1132 x1 = xAxis->GetBinCenter(trackPtBin);
1133 }
1134 else {
1135 y0 = fFractionHists[species]->GetBinContent(trackPtBin, jetPtBin, centBin);
1136 y0errStat = fFractionHists[species]->GetBinError(trackPtBin, jetPtBin, centBin);
1137 y0errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin, jetPtBin, centBin)
1138 : 0.;
1139 x0 = xAxis->GetBinCenter(trackPtBin);
1140 y1 = fFractionHists[species]->GetBinContent(trackPtBin + 1, jetPtBin, centBin);
1141 y1errStat = fFractionHists[species]->GetBinError(trackPtBin + 1, jetPtBin, centBin);
1142 y1errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin + 1, jetPtBin, centBin)
1143 : 0.;
1144 x1 = xAxis->GetBinCenter(trackPtBin + 1);
1145 }
1146
1147 // Per construction: x0 < trackPt < x1
1148 fraction = y0 + (trackPt - x0) * ((y1 - y0) / (x1 - x0));
1149 fractionErrorStat = y0errStat + (trackPt - x0) * ((y1errStat - y0errStat) / (x1 - x0));
1150 fractionErrorSys = fFractionSysErrorHists[species] ? (y0errSys + (trackPt - x0) * ((y1errSys - y0errSys) / (x1 - x0))) : 0.;
1151 }
1152
1153 return kTRUE;
1154}
1155
1156
1157//_____________________________________________________________________________
1158Bool_t AliAnalysisTaskPID::GetParticleFractions(const Double_t trackPt, const Double_t jetPt, const Double_t centralityPercentile,
1159 Double_t* prob, const Int_t smearSpeciesByError,
1160 const Int_t takeIntoAccountSpeciesSysError, const Bool_t uniformSystematicError) const
1161{
1162 // Fills the particle fractions for the given trackPt, jetPt and centrality into "prob".
1163 // Use jetPt = -1 for inclusive spectra and centralityPercentile = -1 for pp.
1164 // If smearSpeciesByError is >= 0 && < AliPID::kSPECIES, the returned fractions will be a random number distributed
1165 // with a gauss with mean being the corresponding particle fraction and sigma it's error for the considered species
1166 // "smearSpeciesByError".
1167 // Note that in this case the fractions for all species will NOT sum up to 1!
1168 // Thus, all other species fractions will be re-scaled weighted with their corresponding statistical error.
1169 // A similar procedure is used for "takeIntoAccountSpeciesSysError": The systematic error of the corresponding species
1170 // is used to generate a random number with uniform distribution in [mean - sysError, mean + sysError] for the new mean
1171 // (in cace of uniformSystematicError = kTRUE, otherwise it will be a gaus(mean, sysError)),
1172 // then the other species will be re-scaled according to their systematic errors.
1173 // First, the systematic error uncertainty procedure will be performed (that is including re-scaling), then the statistical
1174 // uncertainty procedure.
1175 // On success, kTRUE is returned.
1176
1177 if (!prob || smearSpeciesByError >= AliPID::kSPECIES || takeIntoAccountSpeciesSysError >= AliPID::kSPECIES)
1178 return kFALSE;
1179
1180 Double_t probTemp[AliPID::kSPECIES];
1181 Double_t probErrorStat[AliPID::kSPECIES];
1182 Double_t probErrorSys[AliPID::kSPECIES];
1183
1184 Bool_t success = kTRUE;
1185 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kElectron,
1186 probTemp[AliPID::kElectron], probErrorStat[AliPID::kElectron],
1187 probErrorSys[AliPID::kElectron]);
1188 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kMuon,
1189 probTemp[AliPID::kMuon], probErrorStat[AliPID::kMuon], probErrorSys[AliPID::kMuon]);
1190 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kPion,
1191 probTemp[AliPID::kPion], probErrorStat[AliPID::kPion], probErrorSys[AliPID::kPion]);
1192 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kKaon,
1193 probTemp[AliPID::kKaon], probErrorStat[AliPID::kKaon], probErrorSys[AliPID::kKaon]);
1194 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kProton,
1195 probTemp[AliPID::kProton], probErrorStat[AliPID::kProton], probErrorSys[AliPID::kProton]);
1196
1197 if (!success)
1198 return kFALSE;
1199
1200 // If desired, take into account the systematic error of the corresponding species and re-generate probTemp accordingly
1201 if (takeIntoAccountSpeciesSysError >= 0) {
1202 // Generate random fraction of the considered species "smearSpeciesByError" according to mean and sigma
1203 Double_t generatedFraction = uniformSystematicError
1204 ? fRandom->Rndm() * 2. * probErrorSys[takeIntoAccountSpeciesSysError]
1205 - probErrorSys[takeIntoAccountSpeciesSysError]
1206 + probTemp[takeIntoAccountSpeciesSysError]
1207 : fRandom->Gaus(probTemp[takeIntoAccountSpeciesSysError],
1208 probErrorSys[takeIntoAccountSpeciesSysError]);
1209
1210 // Catch cases with invalid fraction (can happen for large errors), i.e. fraction < 0 or > 1
1211 if (generatedFraction < 0.)
1212 generatedFraction = 0.;
1213 else if (generatedFraction > 1.)
1214 generatedFraction = 1.;
1215
1216 // Calculate difference from original fraction (original fractions sum up to 1!)
1217 Double_t deltaFraction = generatedFraction - probTemp[takeIntoAccountSpeciesSysError];
1218
1219 // Fractions must (including errors) lie inside [0,1] -> Adapt weights accordingly by setting the errors
1220 if (deltaFraction > 0) {
1221 // Some part will be SUBTRACTED from the other fractions
1222 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1223 if (probTemp[i] - probErrorSys[i] < 0)
1224 probErrorSys[i] = probTemp[i];
1225 }
1226 }
1227 else {
1228 // Some part will be ADDED to the other fractions
1229 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1230 if (probTemp[i] + probErrorSys[i] > 1)
1231 probErrorSys[i] = 1. - probTemp[i];
1232 }
1233 }
1234
1235 // Compute summed weight of all fractions except for the considered one
1236 Double_t summedWeight = 0.;
1237 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1238 if (i != takeIntoAccountSpeciesSysError)
1239 summedWeight += probErrorSys[i];
1240 }
1241
1242 // Compute the weight for the other species
1243 /*
1244 if (summedWeight <= 1e-13) {
1245 // If this happens for some reason (it should not!), just assume flat weight
1246 printf("Error: summedWeight (sys error) ~ 0 for trackPt %f, jetPt %f, centralityPercentile %f. Setting flat weight!\n",
1247 trackPt, jetPt, centralityPercentile);
1248 }*/
1249
1250 Double_t weight[AliPID::kSPECIES];
1251 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1252 if (i != takeIntoAccountSpeciesSysError) {
1253 if (summedWeight > 1e-13)
1254 weight[i] = probErrorSys[i] / summedWeight;
1255 else
1256 weight[i] = probErrorSys[i] / (AliPID::kSPECIES - 1);
1257 }
1258 }
1259
1260 // For the final generated fractions, set the generated value for the considered species
1261 // and the generated value minus delta times statistical weight
1262 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1263 if (i != takeIntoAccountSpeciesSysError)
1264 probTemp[i] = probTemp[i] - weight[i] * deltaFraction;
1265 else
1266 probTemp[i] = generatedFraction;
1267 }
1268 }
1269
1270 // Using the values of probTemp (either the original ones or those after taking into account the systematic error),
1271 // calculate the final fractions - if the statistical error is to be taken into account, smear the corresponding
1272 // fraction. If not, just write probTemp to the final result array.
1273 if (smearSpeciesByError >= 0) {
1274 // Generate random fraction of the considered species "smearSpeciesByError" according to mean and sigma
1275 Double_t generatedFraction = fRandom->Gaus(probTemp[smearSpeciesByError], probErrorStat[smearSpeciesByError]);
1276
1277 // Catch cases with invalid fraction (can happen for large errors), i.e. fraction < 0 or > 1
1278 if (generatedFraction < 0.)
1279 generatedFraction = 0.;
1280 else if (generatedFraction > 1.)
1281 generatedFraction = 1.;
1282
1283 // Calculate difference from original fraction (original fractions sum up to 1!)
1284 Double_t deltaFraction = generatedFraction - probTemp[smearSpeciesByError];
1285
1286 // Fractions must (including errors) lie inside [0,1] -> Adapt weights accordingly by setting the errors
1287 if (deltaFraction > 0) {
1288 // Some part will be SUBTRACTED from the other fractions
1289 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1290 if (probTemp[i] - probErrorStat[i] < 0)
1291 probErrorStat[i] = probTemp[i];
1292 }
1293 }
1294 else {
1295 // Some part will be ADDED to the other fractions
1296 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1297 if (probTemp[i] + probErrorStat[i] > 1)
1298 probErrorStat[i] = 1. - probTemp[i];
1299 }
1300 }
1301
1302 // Compute summed weight of all fractions except for the considered one
1303 Double_t summedWeight = 0.;
1304 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1305 if (i != smearSpeciesByError)
1306 summedWeight += probErrorStat[i];
1307 }
1308
1309 // Compute the weight for the other species
1310 /*
1311 if (summedWeight <= 1e-13) {
1312 // If this happens for some reason (it should not!), just assume flat weight
1313 printf("Error: summedWeight (stat error) ~ 0 for trackPt %f, jetPt %f, centralityPercentile %f. Setting flat weight!\n",
1314 trackPt, jetPt, centralityPercentile);
1315 }*/
1316
1317 Double_t weight[AliPID::kSPECIES];
1318 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1319 if (i != smearSpeciesByError) {
1320 if (summedWeight > 1e-13)
1321 weight[i] = probErrorStat[i] / summedWeight;
1322 else
1323 weight[i] = probErrorStat[i] / (AliPID::kSPECIES - 1);
1324 }
1325 }
1326
1327 // For the final generated fractions, set the generated value for the considered species
1328 // and the generated value minus delta times statistical weight
1329 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1330 if (i != smearSpeciesByError)
1331 prob[i] = probTemp[i] - weight[i] * deltaFraction;
1332 else
1333 prob[i] = generatedFraction;
1334 }
1335 }
1336 else {
1337 // Just take the generated values
1338 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
1339 prob[i] = probTemp[i];
1340 }
1341
1342
1343 // Should already be normalised, but make sure that it really is:
1344 Double_t probSum = 0.;
1345 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1346 probSum += prob[i];
1347 }
1348
1349 if (probSum <= 0)
1350 return kFALSE;
1351
1352 if (TMath::Abs(probSum - 1.0) > 1e-4) {
1353 printf("Warning: Re-normalising sum of fractions: Sum is %e\n", probSum);
1354 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1355 prob[i] /= probSum;
1356 }
1357 }
1358
1359 return kTRUE;
1360}
1361
1362
1363//_____________________________________________________________________________
1364const TH3D* AliAnalysisTaskPID::GetParticleFractionHisto(const Int_t species, const Bool_t sysError) const
1365{
1366 if (species < AliPID::kElectron || species > AliPID::kProton)
1367 return 0x0;
1368
1369 return sysError ? fFractionSysErrorHists[species] : fFractionHists[species];
1370}
1371
1372
1373//_____________________________________________________________________________
1374Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(Int_t motherPDG, Double_t motherGenPt)
1375{
1376 // Strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
1377 // -> Based on function in PWGJE/AliAnalysisTaskFragmentationFunction, which uses
1378 // the following data from CMS pp @ 7 TeV inclusive (JHEP 05 (2011) 064)
1379
1380 Double_t fac = 1;
1381
1382 const Int_t absMotherPDG = TMath::Abs(motherPDG);
1383
1384 if (absMotherPDG == 310 || absMotherPDG == 321) { // K0s / K+ / K-
1385 if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.768049;
1386 else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.732933;
1387 else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.650298;
1388 else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.571332;
1389 else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.518734;
1390 else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.492543;
1391 else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.482704;
1392 else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.488056;
1393 else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.488861;
1394 else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.492862;
1395 else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.504332;
1396 else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.501858;
1397 else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.512970;
1398 else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.524131;
1399 else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.539130;
1400 else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.554101;
1401 else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.560348;
1402 else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.568869;
1403 else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.583310;
1404 else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.604818;
1405 else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.632630;
1406 else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.710070;
1407 else if(6.00 <= motherGenPt && motherGenPt < 8.00) fac = 0.736365;
1408 else if(8.00 <= motherGenPt && motherGenPt < 10.00) fac = 0.835865;
1409 }
1410
1411 if (absMotherPDG == 3122) { // Lambda
1412 if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.645162;
1413 else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.627431;
1414 else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.457136;
1415 else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.384369;
1416 else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.330597;
1417 else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.309571;
1418 else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.293620;
1419 else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.283709;
1420 else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.282047;
1421 else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.277261;
1422 else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.275772;
1423 else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.280726;
1424 else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.288540;
1425 else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.288315;
1426 else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.296619;
1427 else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.302993;
1428 else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.338121;
1429 else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.349800;
1430 else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.356802;
1431 else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.391202;
1432 else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.422573;
1433 else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.573815;
1434 else if(6.00 <= motherGenPt && motherGenPt < 8.00) fac = 0.786984;
1435 else if(8.00 <= motherGenPt && motherGenPt < 10.00) fac = 1.020021;
1436 }
1437
1438 if (absMotherPDG == 3312 || absMotherPDG == 3322) { // xi
1439 if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.666620;
1440 else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.575908;
1441 else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.433198;
1442 else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.340901;
1443 else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.290896;
1444 else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.236074;
1445 else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.218681;
1446 else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.207763;
1447 else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.222848;
1448 else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.208806;
1449 else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.197275;
1450 else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.183645;
1451 else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.188788;
1452 else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.188282;
1453 else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.207442;
1454 else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.240388;
1455 else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.241916;
1456 else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.208276;
1457 else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.234550;
1458 else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.251689;
1459 else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.310204;
1460 else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.343492;
1461 }
1462
1463 const Double_t weight = 1. / fac;
1464
1465 return weight;
1466}
1467
1468
1469//_____________________________________________________________________________
1470Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(AliMCEvent* mcEvent, AliMCParticle* daughter)
1471{
1472 // Strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
1473 // -> Based on function in PWGJE/AliAnalysisTaskFragmentationFunction
1474
1475 if (!mcEvent)
1476 return 1.;
1477
1478 AliMCParticle* currentMother = daughter;
1479 AliMCParticle* currentDaughter = daughter;
1480
1481
1482 // find first primary mother K0s, Lambda or Xi
1483 while(1) {
1484 Int_t daughterPDG = currentDaughter->PdgCode();
1485
1486 Int_t motherLabel = currentDaughter->GetMother();
1487 if(motherLabel >= mcEvent->GetNumberOfTracks()){ // protection
1488 currentMother = currentDaughter;
1489 break;
1490 }
1491
1492 currentMother = (AliMCParticle*)mcEvent->GetTrack(motherLabel);
1493
1494 if (!currentMother) {
1495 currentMother = currentDaughter;
1496 break;
1497 }
1498
1499 Int_t motherPDG = currentMother->PdgCode();
1500
1501 // phys. primary found ?
1502 if (mcEvent->IsPhysicalPrimary(motherLabel))
1503 break;
1504
1505 if (TMath::Abs(daughterPDG) == 321) {
1506 // K+/K- e.g. from phi (ref data not feeddown corrected)
1507 currentMother = currentDaughter;
1508 break;
1509 }
1510 if (TMath::Abs(motherPDG) == 310) {
1511 // K0s e.g. from phi (ref data not feeddown corrected)
1512 break;
1513 }
1514 if (TMath::Abs(motherPDG) == 3212 && TMath::Abs(daughterPDG) == 3122) {
1515 // Mother Sigma0, daughter Lambda (this case not included in feeddown corr.)
1516 currentMother = currentDaughter;
1517 break;
1518 }
1519
1520 currentDaughter = currentMother;
1521 }
1522
1523
1524 Int_t motherPDG = currentMother->PdgCode();
1525 Double_t motherGenPt = currentMother->Pt();
1526
1527 return GetMCStrangenessFactorCMS(motherPDG, motherGenPt);
1528}
1529
1530
1531// _________________________________________________________________________________
1532Bool_t AliAnalysisTaskPID::IsSecondaryWithStrangeMotherMC(AliMCEvent* mcEvent, Int_t partLabel)
1533{
1534 // Check whether particle is a secondary with strange mother, i.e. returns kTRUE if a strange mother is found
1535 // and the particle is NOT a physical primary. In all other cases kFALSE is returned
1536
1537 if (!mcEvent || partLabel < 0)
1538 return kFALSE;
1539
1540 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(partLabel);
1541
1542 if (!part)
1543 return kFALSE;
1544
1545 if (mcEvent->IsPhysicalPrimary(partLabel))
1546 return kFALSE;
1547
1548 Int_t iMother = part->GetMother();
1549 if (iMother < 0)
1550 return kFALSE;
1551
1552
1553 AliMCParticle* partM = (AliMCParticle*)mcEvent->GetTrack(iMother);
1554 if (!partM)
1555 return kFALSE;
1556
1557 Int_t codeM = TMath::Abs(partM->PdgCode());
1558 Int_t mfl = Int_t(codeM / TMath::Power(10, Int_t(TMath::Log10(codeM))));
1559 if (mfl == 3 && codeM != 3) // codeM = 3 is for s quark
1560 return kTRUE;
1561
1562 return kFALSE;
1563}
1564
1565
1566//_____________________________________________________________________________
1567Bool_t AliAnalysisTaskPID::SetParticleFractionHisto(const TH3D* hist, const Int_t species, const Bool_t sysError)
1568{
1569 // Store a clone of hist (containing the particle fractions of the corresponding species with statistical error (sysError = kFALSE)
1570 // or systematic error (sysError = kTRUE), respectively), internally
1571
1572 if (species < AliPID::kElectron || species > AliPID::kProton) {
1573 AliError(Form("Only fractions for species index %d to %d can be set, but not for the requested one: %d", 0,
1574 AliPID::kProton, species));
1575 return kFALSE;
1576 }
1577
1578 if (sysError) {
1579 delete fFractionSysErrorHists[species];
1580
1581 fFractionSysErrorHists[species] = new TH3D(*hist);
1582 }
1583 else {
1584 delete fFractionHists[species];
1585
1586 fFractionHists[species] = new TH3D(*hist);
1587 }
1588
1589 return kTRUE;
1590}
1591
1592
1593//_____________________________________________________________________________
1594Bool_t AliAnalysisTaskPID::SetParticleFractionHistosFromFile(const TString filePathName, const Bool_t sysError)
1595{
1596 // Loads particle fractions for all species from the desired file and returns kTRUE on success.
1597 // The maps are assumed to be of Type TH3D, to sit in the main directory and to have names
1598 // Form("hFraction_%e", AliPID::ParticleName(i)) for sysError = kFALSE and
1599 // Form("hFractionSysError_%e", AliPID::ParticleName(i)) for sysError = kTRUE.
1600
1601 TFile* f = TFile::Open(filePathName.Data());
1602 if (!f) {
1603 std::cout << "Failed to open file with particle fractions \"" << filePathName.Data() << "\"!" << std::endl;
1604 return kFALSE;
1605 }
1606
1607 TH3D* hist = 0x0;
1608 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1609 TString histName = Form("hFraction%s_%s", sysError ? "SysError" : "", AliPID::ParticleName(i));
1610 hist = dynamic_cast<TH3D*>(f->Get(histName.Data()));
1611 if (!hist) {
1612 std::cout << "Failed to load particle fractions for " << histName.Data() << "!";
1613 std::cout << std::endl << "Cleaning up particle fraction histos!" << std::endl;
1614 CleanupParticleFractionHistos();
1615 return kFALSE;
1616 }
1617
1618 if (!SetParticleFractionHisto(hist, i, sysError)) {
1619 std::cout << "Failed to load particle fractions for " << histName.Data() << "!";
1620 std::cout << std::endl << "Cleaning up particle fraction histos!" << std::endl;
1621 CleanupParticleFractionHistos();
1622 return kFALSE;
1623 }
1624 }
1625
1626 delete hist;
1627
1628 return kTRUE;
1629
1630}
1631
1632
1633//_____________________________________________________________________________
1634Int_t AliAnalysisTaskPID::GetRandomParticleTypeAccordingToParticleFractions(const Double_t trackPt, const Double_t jetPt,
1635 const Double_t centralityPercentile,
1636 const Bool_t smearByError,
1637 const Bool_t takeIntoAccountSysError) const
1638{
1639 // Uses the stored histograms with the particle fractions to generate a random particle type according to these fractions.
1640 // In case of problems (e.g. histo missing), AliPID::kUnknown is returned.
1641 // If smearByError is kTRUE, the used fractions will be random numbers distributed with a gauss with mean
1642 // being the corresponding particle fraction and sigma it's error.
1643 // Note that in this case only the fraction of a random species is varied in this way. The other fractions
1644 // will be re-normalised according their statistical errors.
1645 // The same holds for the systematic error of species "takeIntoAccountSpeciesSysError", but the random number will be
1646 // uniformly distributed within [mean - sys, mean + sys] and the re-normalisation will be weighted with the systematic errors.
1647 // Note that the fractions will be calculated first with only the systematic error taken into account (if desired), including
1648 // re-normalisation. Then, the resulting fractions will be used to calculate the final fractions - either with statistical error
1649 // or without. The species, for which the error will be used for smearing, is the same for sys and stat error.
1650
1651 Double_t prob[AliPID::kSPECIES];
1652 Int_t randomSpecies = (smearByError || takeIntoAccountSysError) ? (Int_t)(fRandom->Rndm() * AliPID::kSPECIES) : -1;
1653 Bool_t success = GetParticleFractions(trackPt, jetPt, centralityPercentile, prob, randomSpecies, randomSpecies);
1654
1655 if (!success)
1656 return AliPID::kUnknown;
1657
1658 Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
1659
1660 if (rnd <= prob[AliPID::kPion])
1661 return AliPID::kPion;
1662 else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon])
1663 return AliPID::kKaon;
1664 else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon] + prob[AliPID::kProton])
1665 return AliPID::kProton;
1666 else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon] + prob[AliPID::kProton] + prob[AliPID::kElectron])
1667 return AliPID::kElectron;
1668
1669 return AliPID::kMuon; //else it must be a muon (only species left)
1670}
1671
1672
1673//_____________________________________________________________________________
1674AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::GenerateDetectorResponse(const AliAnalysisTaskPID::ErrorCode errCode,
1675 const Double_t mean, const Double_t sigma,
1676 Double_t* responses, const Int_t nResponses,
1677 const Bool_t usePureGaus)
1678{
1679 // Generate detector response. If a previous generation was not successful or there is something wrong with this signal generation,
1680 // the function will return kFALSE
1681 if (!responses)
1682 return kError;
1683
1684 // Reset response array
1685 for (Int_t i = 0; i < nResponses; i++)
1686 responses[i] = -999;
1687
1688 if (errCode == kError)
1689 return kError;
1690
1691 ErrorCode ownErrCode = kNoErrors;
1692
1693 if (fUseConvolutedGaus && !usePureGaus) {
1694 // In case of convoluted gauss, calculate the probability density only once to save a lot of time!
1695
1696 TH1* hProbDensity = 0x0;
1697 ownErrCode = SetParamsForConvolutedGaus(mean, sigma);
1698 if (ownErrCode == kError)
1699 return kError;
1700
1701 hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
1702
1703 for (Int_t i = 0; i < nResponses; i++) {
1704 responses[i] = hProbDensity->GetRandom();
1705 //responses[i] fConvolutedGausDeltaPrime->GetRandom(); // MUCH slower than using the binned version via the histogram
1706 }
1707 }
1708 else {
1709 for (Int_t i = 0; i < nResponses; i++) {
1710 responses[i] = fRandom->Gaus(mean, sigma);
1711 }
1712 }
1713
1714 // If forwarded error code was a warning (error case has been handled before), return a warning
1715 if (errCode == kWarning)
1716 return kWarning;
1717
1718 return ownErrCode; // Forward success/warning
1719}
1720
1721
1722//_____________________________________________________________________________
1723void AliAnalysisTaskPID::PrintSettings(Bool_t printSystematicsSettings) const
1724{
1725 // Print current settings.
1726
1727 printf("\n\nSettings for task %s:\n", GetName());
1728 printf("Is pPb/Pbp: %d -> %s\n", GetIsPbpOrpPb(), GetIsPbpOrpPb() ? "Adapting vertex cuts" : "Using standard vertex cuts");
1729 printf("Track cuts: %s\n", fTrackFilter ? fTrackFilter->GetTitle() : "-");
1730 printf("Eta cut: %.2f <= |eta| <= %.2f\n", GetEtaAbsCutLow(), GetEtaAbsCutUp());
1731 printf("Phi' cut: %d\n", GetUsePhiCut());
1732
1733 printf("\n");
1734
1735 printf("Centrality estimator: %s\n", GetCentralityEstimator().Data());
1736
1737 printf("\n");
1738
1739 printf("Use MC-ID for signal generation: %d\n", GetUseMCidForGeneration());
1740 printf("Use ITS: %d\n", GetUseITS());
1741 printf("Use TOF: %d\n", GetUseTOF());
1742 printf("Use priors: %d\n", GetUsePriors());
1743 printf("Use TPC default priors: %d\n", GetUseTPCDefaultPriors());
1744 printf("Use convoluted Gauss: %d\n", GetUseConvolutedGaus());
1745 printf("Accuracy of non-Gaussian tail: %e\n", GetAccuracyNonGaussianTail());
1746 printf("Take into account muons: %d\n", GetTakeIntoAccountMuons());
1747 printf("\nParams for transition from gauss to asymmetric shape:\n");
1748 printf("[0]: %e\n", GetConvolutedGaussTransitionPar(0));
1749 printf("[1]: %e\n", GetConvolutedGaussTransitionPar(1));
1750 printf("[2]: %e\n", GetConvolutedGaussTransitionPar(2));
1751
1752 printf("\n");
1753
1754 printf("Input from other task: %d\n", GetInputFromOtherTask());
1755 printf("Store additional jet information: %d\n", GetStoreAdditionalJetInformation());
1756 printf("Store centrality percentile: %d", GetStoreCentralityPercentile());
1757
1758 if (printSystematicsSettings)
1759 PrintSystematicsSettings();
1760 else
1761 printf("\n\n\n");
1762}
1763
1764
1765//_____________________________________________________________________________
1766void AliAnalysisTaskPID::PrintSystematicsSettings() const
1767{
1768 // Print current settings for systematic studies.
1769
1770 printf("\n\nSettings for systematics for task %s:\n", GetName());
1771 printf("Splines:\t%f\n", GetSystematicScalingSplines());
1772 printf("EtaCorrMomThr:\t%f\n", GetSystematicScalingEtaCorrectionMomentumThr());
1773 printf("EtaCorrLowP:\t%f\n", GetSystematicScalingEtaCorrectionLowMomenta());
1774 printf("EtaCorrHighP:\t%f\n", GetSystematicScalingEtaCorrectionHighMomenta());
1775 printf("SigmaPara:\t%f\n", GetSystematicScalingEtaSigmaPara());
1776 printf("MultCorr:\t%f\n", GetSystematicScalingMultCorrection());
1777
1778 printf("\n\n");
1779}
1780
1781
1782//_____________________________________________________________________________
1783Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePDGcode, Double_t centralityPercentile,
1784 Double_t jetPt)
1785{
1786 // Process the track (generate expected response, fill histos, etc.).
1787 // particlePDGcode == 0 means data. Otherwise, the corresponding MC ID will be assumed.
1788
1789 //Printf("Debug: Task %s is starting to process track: dEdx %f, pTPC %f, eta %f, ncl %d\n", GetName(), track->GetTPCsignal(), track->GetTPCmomentum(),
1790 // track->Eta(), track->GetTPCsignalN());
1791
1792 const Bool_t isMC = (particlePDGcode == 0) ? kFALSE : kTRUE;
1793
1794 Int_t binMC = -1;
1795
1796 if (isMC) {
1797 if (TMath::Abs(particlePDGcode) == 211) {//Pion
1798 binMC = 3;
1799 }
1800 else if (TMath::Abs(particlePDGcode) == 321) {//Kaon
1801 binMC = 1;
1802 }
1803 else if (TMath::Abs(particlePDGcode) == 2212) {//Proton
1804 binMC = 4;
1805 }
1806 else if (TMath::Abs(particlePDGcode) == 11) {//Electron
1807 binMC = 0;
1808 }
1809 else if (TMath::Abs(particlePDGcode) == 13) {//Muon
1810 binMC = 2;
1811 }
1812 else // In MC-ID case, set to underflow bin such that the response from this track is only used for unidentified signal generation
1813 // or signal generation with PID response and the track is still there (as in data) - e.g. relevant w.r.t. deuterons.
1814 // This is important to be as much as possible consistent with data. And the tracks can still be removed by disabling the
1815 // underflow bin for the projections
1816 binMC = -1;
1817 }
1818
1819 // Momenta
1820 //Double_t p = track->GetP();
1821 //Double_t pTPC = track->GetTPCmomentum();
1822 Double_t pT = track->Pt();
1823
1824 Double_t z = -1, xi = -1;
1825 GetJetTrackObservables(pT, jetPt, z, xi);
1826
1827
1828 Double_t trackCharge = track->Charge();
1829
1830 // TPC signal
1831 Double_t dEdxTPC = fPIDResponse->IsTunedOnData() ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
1832
1833 if (dEdxTPC <= 0) {
1834 Printf("Skipping track with strange dEdx value: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), track->GetTPCmomentum(),
1835 track->Eta(), track->GetTPCsignalN());
1836 return kFALSE;
1837 }
1838
1839
1840
1841
1842 Double_t dEdxEl, dEdxKa, dEdxPi, dEdxMu, dEdxPr;
1843 Double_t sigmaEl, sigmaKa, sigmaPi, sigmaMu, sigmaPr;
1844
1845 if (fDoAnySystematicStudiesOnTheExpectedSignal) {
1846 // Get the uncorrected signal first and the corresponding correction factors.
1847 // Then modify the correction factors and properly recalculate the corrected dEdx
1848
1849 // Get pure spline values for dEdx_expected, without any correction
1850 dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
1851 dEdxKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
1852 dEdxPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
1853 dEdxMu = !fTakeIntoAccountMuons ? -1 :
1854 fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
1855 dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
1856
1857 // Scale splines, if desired
1858 if (TMath::Abs(fSystematicScalingSplines - 1.0) > fgkEpsilon) {
1859 dEdxEl *= fSystematicScalingSplines;
1860 dEdxKa *= fSystematicScalingSplines;
1861 dEdxPi *= fSystematicScalingSplines;
1862 dEdxMu *= fTakeIntoAccountMuons ? fSystematicScalingSplines : 1.;
1863 dEdxPr *= fSystematicScalingSplines;
1864 }
1865
1866 // Get the eta correction factors for the (modified) expected dEdx
1867 Double_t etaCorrEl = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxEl) : 1.;
1868 Double_t etaCorrKa = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxKa) : 1.;
1869 Double_t etaCorrPi = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPi) : 1.;
1870 Double_t etaCorrMu = fTakeIntoAccountMuons && !fPIDResponse->UseTPCEtaCorrection() ?
1871 fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxMu) : 1.;
1872 Double_t etaCorrPr = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPr) : 1.;
1873
1874 // Scale eta correction factors, if desired (and eta correction maps are to be used, otherwise it is not possible!)
1875 if (fPIDResponse->UseTPCEtaCorrection() &&
1876 (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon ||
1877 TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon)) {
1878 // Since we do not want to scale the splines with this, but only the eta variation, only scale the deviation of the correction factor!
1879 // 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!
1880
1881
1882 // Due to additional azimuthal effects, there is an additional eta dependence for low momenta which is not corrected successfully so far.
1883 // One can assign a different (higher) systematic scale factor for this low-p region and a threshold which separates low- and high-p.
1884 // 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
1885 Double_t usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionHighMomenta;
1886
1887 if (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - fSystematicScalingEtaCorrectionLowMomenta) > fgkEpsilon) {
1888 const Double_t pTPC = track->GetTPCmomentum();
1889 const Double_t fractionHighMomentumScaleFactor = 0.5 * (1. + TMath::Erf((pTPC - fSystematicScalingEtaCorrectionMomentumThr) / 0.1));
1890 usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionLowMomenta * (1 - fractionHighMomentumScaleFactor)
1891 + fSystematicScalingEtaCorrectionHighMomenta * fractionHighMomentumScaleFactor;
1892 }
1893
1894 etaCorrEl = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrEl - 1.0);
1895 etaCorrKa = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrKa - 1.0);
1896 etaCorrPi = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPi - 1.0);
1897 etaCorrMu = fTakeIntoAccountMuons ? (1.0 + usedSystematicScalingEtaCorrection * (etaCorrMu - 1.0)) : 1.0;
1898 etaCorrPr = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPr - 1.0);
1899 }
1900
1901 // Get the multiplicity correction factors for the (modified) expected dEdx
1902 const Int_t currEvtMultiplicity = fPIDResponse->GetTPCResponse().GetCurrentEventMultiplicity();
1903
1904 Double_t multiplicityCorrEl = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
1905 dEdxEl * etaCorrEl, currEvtMultiplicity) : 1.;
1906 Double_t multiplicityCorrKa = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
1907 dEdxKa * etaCorrKa, currEvtMultiplicity) : 1.;
1908 Double_t multiplicityCorrPi = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
1909 dEdxPi * etaCorrPi, currEvtMultiplicity) : 1.;
1910 Double_t multiplicityCorrMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCMultiplicityCorrection() ?
1911 fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track, dEdxMu * etaCorrMu, currEvtMultiplicity) : 1.;
1912 Double_t multiplicityCorrPr = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
1913 dEdxPr * etaCorrPr, currEvtMultiplicity) : 1.;
1914
1915 Double_t multiplicityCorrSigmaEl = fPIDResponse->UseTPCMultiplicityCorrection() ?
1916 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxEl * etaCorrEl, currEvtMultiplicity) : 1.;
1917 Double_t multiplicityCorrSigmaKa = fPIDResponse->UseTPCMultiplicityCorrection() ?
1918 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxKa * etaCorrKa, currEvtMultiplicity) : 1.;
1919 Double_t multiplicityCorrSigmaPi = fPIDResponse->UseTPCMultiplicityCorrection() ?
1920 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxPi * etaCorrPi, currEvtMultiplicity) : 1.;
1921 Double_t multiplicityCorrSigmaMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCMultiplicityCorrection() ?
1922 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxMu * etaCorrMu, currEvtMultiplicity) : 1.;
1923 Double_t multiplicityCorrSigmaPr = fPIDResponse->UseTPCMultiplicityCorrection() ?
1924 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxPr * etaCorrPr, currEvtMultiplicity) : 1.;
1925
1926 // Scale multiplicity correction factors, if desired (and multiplicity correction functions are to be used, otherwise it is not possible!)
1927 if (fPIDResponse->UseTPCMultiplicityCorrection() && TMath::Abs(fSystematicScalingMultCorrection - 1.0) > fgkEpsilon) {
1928 // Since we do not want to scale the splines with this, but only the multiplicity variation, only scale the deviation of the correction factor!
1929 // 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!
1930
1931 multiplicityCorrEl = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrEl - 1.0);
1932 multiplicityCorrKa = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrKa - 1.0);
1933 multiplicityCorrPi = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrPi - 1.0);
1934 multiplicityCorrMu = fTakeIntoAccountMuons ? (1.0 + fSystematicScalingMultCorrection * (multiplicityCorrMu - 1.0)) : 1.0;
1935 multiplicityCorrPr = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrPr - 1.0);
1936 }
1937
1938 // eta correction must be enabled in order to use the new sigma parametrisation maps. Since this is the absolute sigma
1939 // for a track calculated with the unscaled paramaters, we have to devide out dEdxExpectedEtaCorrected and then need
1940 // to scale with the multiplicitySigmaCorrFactor * fSystematicScalingEtaSigmaPara. In the end, one has to scale with the
1941 // (modified) dEdx to get the absolute sigma
1942 // This means there is no extra parameter for the multiplicitySigmaCorrFactor, but only for the sigma map itself.
1943 // This is valid, since it appears only as a product. One has to assume a larger systematic shift in case of additional
1944 // multiplicity dependence....
1945 Double_t sigmaRelEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1946 / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1947 * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaEl;
1948
1949 Double_t sigmaRelKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1950 / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1951 * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaKa;
1952
1953 Double_t sigmaRelPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1954 / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1955 * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaPi;
1956
1957 Double_t sigmaRelMu = fTakeIntoAccountMuons ?
1958 fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1959 / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1960 * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaMu
1961 : 999.;
1962
1963 Double_t sigmaRelPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1964 / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1965 * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaPr;
1966
1967 // Now scale the (possibly modified) spline values with the (possibly modified) correction factors
1968 dEdxEl *= etaCorrEl * multiplicityCorrEl;
1969 dEdxKa *= etaCorrKa * multiplicityCorrKa;
1970 dEdxPi *= etaCorrPi * multiplicityCorrPi;
1971 dEdxMu *= etaCorrMu * multiplicityCorrMu;
1972 dEdxPr *= etaCorrPr * multiplicityCorrPr;
1973
1974 // Finally, get the absolute sigma
1975 sigmaEl = sigmaRelEl * dEdxEl;
1976 sigmaKa = sigmaRelKa * dEdxKa;
1977 sigmaPi = sigmaRelPi * dEdxPi;
1978 sigmaMu = sigmaRelMu * dEdxMu;
1979 sigmaPr = sigmaRelPr * dEdxPr;
1980
1981 }
1982 else {
1983 // No systematic studies on expected signal - just take it as it comve from the TPCPIDResponse
1984 dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
1985 fPIDResponse->UseTPCEtaCorrection(),
1986 fPIDResponse->UseTPCMultiplicityCorrection());
1987 dEdxKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
1988 fPIDResponse->UseTPCEtaCorrection(),
1989 fPIDResponse->UseTPCMultiplicityCorrection());
1990 dEdxPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
1991 fPIDResponse->UseTPCEtaCorrection(),
1992 fPIDResponse->UseTPCMultiplicityCorrection());
1993 dEdxMu = !fTakeIntoAccountMuons ? -1 :
1994 fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
1995 fPIDResponse->UseTPCEtaCorrection(),
1996 fPIDResponse->UseTPCMultiplicityCorrection());
1997 dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
1998 fPIDResponse->UseTPCEtaCorrection(),
1999 fPIDResponse->UseTPCMultiplicityCorrection());
2000
2001 sigmaEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2002 fPIDResponse->UseTPCEtaCorrection(),
2003 fPIDResponse->UseTPCMultiplicityCorrection());
2004 sigmaKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2005 fPIDResponse->UseTPCEtaCorrection(),
2006 fPIDResponse->UseTPCMultiplicityCorrection());
2007 sigmaPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2008 fPIDResponse->UseTPCEtaCorrection(),
2009 fPIDResponse->UseTPCMultiplicityCorrection());
2010 sigmaMu = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2011 fPIDResponse->UseTPCEtaCorrection(),
2012 fPIDResponse->UseTPCMultiplicityCorrection());
2013 sigmaPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2014 fPIDResponse->UseTPCEtaCorrection(),
2015 fPIDResponse->UseTPCMultiplicityCorrection());
2016 }
2017 /*OLD with deltaSpecies
2018 Double_t deltaElectron = dEdxTPC - dEdxEl;
2019 Double_t deltaKaon = dEdxTPC - dEdxKa;
2020 Double_t deltaPion = dEdxTPC - dEdxPi;
2021 Double_t deltaProton = dEdxTPC - dEdxPr;
2022 */
2023
2024 Double_t deltaPrimeElectron = (dEdxEl > 0) ? dEdxTPC / dEdxEl : -1;
2025 if (dEdxEl <= 0) {
2026 Printf("Error: Expected TPC signal <= 0 for electron hypothesis");
2027 return kFALSE;
2028 }
2029
2030 Double_t deltaPrimeKaon = (dEdxKa > 0) ? dEdxTPC / dEdxKa : -1;
2031 if (dEdxKa <= 0) {
2032 Printf("Error: Expected TPC signal <= 0 for kaon hypothesis");
2033 return kFALSE;
2034 }
2035
2036 Double_t deltaPrimePion = (dEdxPi > 0) ? dEdxTPC / dEdxPi : -1;
2037 if (dEdxPi <= 0) {
2038 Printf("Error: Expected TPC signal <= 0 for pion hypothesis");
2039 return kFALSE;
2040 }
2041
2042 Double_t deltaPrimeProton = (dEdxPr > 0) ? dEdxTPC / dEdxPr : -1;
2043 if (dEdxPr <= 0) {
2044 Printf("Error: Expected TPC signal <= 0 for proton hypothesis");
2045 return kFALSE;
2046 }
2047
2048 /*TODO for TOF
2049 // TOF signal
2050 Double_t times[AliPID::kSPECIES];
2051 track->GetIntegratedTimes(times);
2052 AliTOFPIDResponse tofPIDResponse = fPIDResponse->GetTOFResponse();
2053 Float_t electronDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kElectron);
2054 Float_t pionDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kPion);
2055 Float_t kaonDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kKaon);
2056 Float_t protonDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kProton);
2057 */
2058
2059
2060 // Use probabilities to weigh the response generation for the different species.
2061 // Also determine the most probable particle type.
2062 Double_t prob[AliPID::kSPECIESC];
2063 for (Int_t i = 0; i < AliPID::kSPECIESC; i++)
2064 prob[i] = 0;
2065
2066 fPIDcombined->ComputeProbabilities(track, fPIDResponse, prob);
2067
2068 // Bug: One can set the number of species for PIDcombined, but PIDcombined will call PIDresponse, which writes without testing
2069 // the probs for kSPECIESC (including light nuclei) into the array.
2070 // In this case, when only kSPECIES are considered, the probabilities have to be rescaled!
2071 for (Int_t i = AliPID::kSPECIES; i < AliPID::kSPECIESC; i++)
2072 prob[i] = 0;
2073
2074 // If muons are not to be taken into account, just set their probability to zero and normalise the remaining probabilities
2075 if (!fTakeIntoAccountMuons)
2076 prob[AliPID::kMuon] = 0;
2077
2078 Double_t probSum = 0.;
2079 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2080 probSum += prob[i];
2081
2082 if (probSum > 0) {
2083 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2084 prob[i] /= probSum;
2085 }
2086
2087 if (!isMC) {
2088 // If there is no MC information, take the most probable species for the ID
2089 Float_t max = 0.;
2090 Int_t maxIndex = -1;
2091 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2092 if (prob[i] > max) {
2093 max = prob[i];
2094 maxIndex = i;
2095 }
2096 }
2097
2098 // Translate from AliPID numbering to numbering of this class
2099 if (max > 0) {
2100 if (maxIndex == AliPID::kElectron)
2101 binMC = 0;
2102 else if (maxIndex == AliPID::kKaon)
2103 binMC = 1;
2104 else if (maxIndex == AliPID::kMuon)
2105 binMC = 2;
2106 else if (maxIndex == AliPID::kPion)
2107 binMC = 3;
2108 else if (maxIndex == AliPID::kProton)
2109 binMC = 4;
2110 else
2111 binMC = -1;
2112 }
2113 else {
2114 // Only take track into account for expectation values, if valid pid response is available.. Otherwise: Set to underflow bin.
2115 binMC = -1;
2116 }
2117 }
2118
2119 /*
2120 //For testing: Swap p<->pT to analyse pure p-dependence => Needs to be removed later
2121 Double_t temp = pT;
2122 pT = pTPC;
2123 pTPC = temp;
2124 */
2125
2126
2127 Double_t entry[fStoreAdditionalJetInformation ? kDataNumAxes : kDataNumAxes - fgkNumJetAxes];
2128 entry[kDataMCID] = binMC;
2129 entry[kDataSelectSpecies] = 0;
2130 entry[kDataPt] = pT;
2131 entry[kDataDeltaPrimeSpecies] = deltaPrimeElectron;
2132 entry[kDataCentrality] = centralityPercentile;
2133
2134 if (fStoreAdditionalJetInformation) {
2135 entry[kDataJetPt] = jetPt;
2136 entry[kDataZ] = z;
2137 entry[kDataXi] = xi;
2138 }
2139
2140 entry[GetIndexOfChargeAxisData()] = trackCharge;
2141
2142 /*OLD with TOF, p_TPC_Inner and p_vertex and deltaSpecies
2143 // MC PID, SelectSpecies, P(TPC_inner), pT, p(Vertex), deltaSpecies, deltaPrimeSpecies, deltaTOFspecies
2144 Double_t entry[8] = { binMC, 0, pTPC, pT, p, deltaElectron, deltaPrimeElectron, electronDeltaTOF };
2145 */
2146
2147 fhPIDdataAll->Fill(entry);
2148
2149 entry[kDataSelectSpecies] = 1;
2150 //OLD with deltaSpecies entry[5] = deltaKaon;
2151 entry[kDataDeltaPrimeSpecies] = deltaPrimeKaon;
2152 //TODO for TOF entry[7] = kaonDeltaTOF;
2153 fhPIDdataAll->Fill(entry);
2154
2155 entry[kDataSelectSpecies] = 2;
2156 //OLD with deltaSpecies entry[5] = deltaPion;
2157 entry[kDataDeltaPrimeSpecies] = deltaPrimePion;
2158 //TODO for TOF entry[7] = pionDeltaTOF;
2159 fhPIDdataAll->Fill(entry);
2160
2161 entry[kDataSelectSpecies] = 3;
2162 //OLD with deltaSpecies entry[5] = deltaProton;
2163 entry[kDataDeltaPrimeSpecies] = deltaPrimeProton;
2164 //TODO for TOF entry[7] = protonDeltaTOF;
2165 fhPIDdataAll->Fill(entry);
2166
2167
2168 // Construct the expected shape for the transition p -> pT
2169
2170 Double_t genEntry[fStoreAdditionalJetInformation ? kGenNumAxes : kGenNumAxes - fgkNumJetAxes];
2171 genEntry[kGenMCID] = binMC;
2172 genEntry[kGenSelectSpecies] = 0;
2173 genEntry[kGenPt] = pT;
2174 genEntry[kGenDeltaPrimeSpecies] = -999;
2175 genEntry[kGenCentrality] = centralityPercentile;
2176
2177 if (fStoreAdditionalJetInformation) {
2178 genEntry[kGenJetPt] = jetPt;
2179 genEntry[kGenZ] = z;
2180 genEntry[kGenXi] = xi;
2181 }
2182
2183 genEntry[GetIndexOfChargeAxisGen()] = trackCharge;
2184
2185 //OLD with deltaSpecies Double_t genEntry[5] = { binMC, 0, pT, -999, -999 }; // MC PID, SelectSpecies, pT, deltaSpecies, deltaPrimeSpecies
2186
2187 // Generate numGenEntries "responses" with fluctuations around the expected values.
2188 // The higher the (transverse) momentum, the more "responses" will be generated in order not to run out of statistics too fast.
2189 Int_t numGenEntries = 10;
2190 if (pT >= 5)
2191 numGenEntries *= 5;
2192 else if (pT >= 2)
2193 numGenEntries *= 2;
2194
2195 // Jets have even worse statistics, therefore, scale numGenEntries further
2196 if (jetPt >= 40)
2197 numGenEntries *= 20;
2198 else if (jetPt >= 20)
2199 numGenEntries *= 10;
2200 else if (jetPt >= 10)
2201 numGenEntries *= 2;
2202
2203
2204 // Do not generate more entries than available in memory!
2205 if (numGenEntries > fgkMaxNumGenEntries)// fgkMaxNumGenEntries = 1000
2206 numGenEntries = fgkMaxNumGenEntries;
2207
2208 ErrorCode errCode = kNoErrors;
2209
2210 // Electrons
2211 errCode = GenerateDetectorResponse(errCode, 1., sigmaEl / dEdxEl, fGenRespElDeltaPrimeEl, numGenEntries);
2212 errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxKa, sigmaEl / dEdxKa, fGenRespElDeltaPrimeKa, numGenEntries);
2213 errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxPi, sigmaEl / dEdxPi, fGenRespElDeltaPrimePi, numGenEntries);
2214 errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxPr, sigmaEl / dEdxPr, fGenRespElDeltaPrimePr, numGenEntries);
2215
2216 // Kaons
2217 errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxEl, sigmaKa / dEdxEl, fGenRespKaDeltaPrimeEl, numGenEntries);
2218 errCode = GenerateDetectorResponse(errCode, 1., sigmaKa / dEdxKa, fGenRespKaDeltaPrimeKa, numGenEntries);
2219 errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxPi, sigmaKa / dEdxPi, fGenRespKaDeltaPrimePi, numGenEntries);
2220 errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxPr, sigmaKa / dEdxPr, fGenRespKaDeltaPrimePr, numGenEntries);
2221
2222 // Pions
2223 errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxEl, sigmaPi / dEdxEl, fGenRespPiDeltaPrimeEl, numGenEntries);
2224 errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxKa, sigmaPi / dEdxKa, fGenRespPiDeltaPrimeKa, numGenEntries);
2225 errCode = GenerateDetectorResponse(errCode, 1., sigmaPi / dEdxPi, fGenRespPiDeltaPrimePi, numGenEntries);
2226 errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxPr, sigmaPi / dEdxPr, fGenRespPiDeltaPrimePr, numGenEntries);
2227
2228 // Muons, if desired
2229 if (fTakeIntoAccountMuons) {
2230 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxEl, sigmaMu / dEdxEl, fGenRespMuDeltaPrimeEl, numGenEntries);
2231 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxKa, sigmaMu / dEdxKa, fGenRespMuDeltaPrimeKa, numGenEntries);
2232 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxPi, sigmaMu / dEdxPi, fGenRespMuDeltaPrimePi, numGenEntries);
2233 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxPr, sigmaMu / dEdxPr, fGenRespMuDeltaPrimePr, numGenEntries);
2234 }
2235
2236 // Protons
2237 errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxEl, sigmaPr / dEdxEl, fGenRespPrDeltaPrimeEl, numGenEntries);
2238 errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxKa, sigmaPr / dEdxKa, fGenRespPrDeltaPrimeKa, numGenEntries);
2239 errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxPi, sigmaPr / dEdxPi, fGenRespPrDeltaPrimePi, numGenEntries);
2240 errCode = GenerateDetectorResponse(errCode, 1., sigmaPr / dEdxPr, fGenRespPrDeltaPrimePr, numGenEntries);
2241
2242
2243 /*OLD with deltaSpecies
2244 // Delta
2245
2246 // Electrons
2247 errCode = GenerateDetectorResponse(errCode, 0., sigmaEl, fGenRespElDeltaEl, numGenEntries, usePureGausForDelta);
2248 errCode = GenerateDetectorResponse(errCode, dEdxEl - dEdxKa, sigmaEl, fGenRespElDeltaKa, numGenEntries, usePureGausForDelta);
2249 errCode = GenerateDetectorResponse(errCode, dEdxEl - dEdxPi, sigmaEl, fGenRespElDeltaPi, numGenEntries, usePureGausForDelta);
2250 errCode = GenerateDetectorResponse(errCode, dEdxEl - dEdxPr, sigmaEl, fGenRespElDeltaPr, numGenEntries, usePureGausForDelta);
2251
2252 // Kaons
2253 errCode = GenerateDetectorResponse(errCode, dEdxKa - dEdxEl, sigmaKa, fGenRespKaDeltaEl, numGenEntries, usePureGausForDelta);
2254 errCode = GenerateDetectorResponse(errCode, 0., sigmaKa, fGenRespKaDeltaKa, numGenEntries, usePureGausForDelta);
2255 errCode = GenerateDetectorResponse(errCode, dEdxKa - dEdxPi, sigmaKa, fGenRespKaDeltaPi, numGenEntries, usePureGausForDelta);
2256 errCode = GenerateDetectorResponse(errCode, dEdxKa - dEdxPr, sigmaKa, fGenRespKaDeltaPr, numGenEntries, usePureGausForDelta);
2257
2258 // Pions
2259 errCode = GenerateDetectorResponse(errCode, dEdxPi - dEdxEl, sigmaPi, fGenRespPiDeltaEl, numGenEntries, usePureGausForDelta);
2260 errCode = GenerateDetectorResponse(errCode, dEdxPi - dEdxKa, sigmaPi, fGenRespPiDeltaKa, numGenEntries, usePureGausForDelta);
2261 errCode = GenerateDetectorResponse(errCode, 0., sigmaPi, fGenRespPiDeltaPi, numGenEntries, usePureGausForDelta);
2262 errCode = GenerateDetectorResponse(errCode, dEdxPi - dEdxPr, sigmaPi, fGenRespPiDeltaPr, numGenEntries, usePureGausForDelta);
2263
2264 // Muons
2265 errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxEl, sigmaMu, fGenRespMuDeltaEl, numGenEntries, usePureGausForDelta);
2266 errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxKa, sigmaMu, fGenRespMuDeltaKa, numGenEntries, usePureGausForDelta);
2267 errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxPi, sigmaMu, fGenRespMuDeltaPi, numGenEntries, usePureGausForDelta);
2268 errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxPr, sigmaMu, fGenRespMuDeltaPr, numGenEntries, usePureGausForDelta);
2269
2270 // Protons
2271 errCode = GenerateDetectorResponse(errCode, dEdxPr - dEdxEl, sigmaPr, fGenRespPrDeltaEl, numGenEntries, usePureGausForDelta);
2272 errCode = GenerateDetectorResponse(errCode, dEdxPr - dEdxKa, sigmaPr, fGenRespPrDeltaKa, numGenEntries, usePureGausForDelta);
2273 errCode = GenerateDetectorResponse(errCode, dEdxPr - dEdxPi, sigmaPr, fGenRespPrDeltaPi, numGenEntries, usePureGausForDelta);
2274 errCode = GenerateDetectorResponse(errCode, 0., sigmaPr, fGenRespPrDeltaPr, numGenEntries, usePureGausForDelta);
2275 */
2276
2277 if (errCode != kNoErrors) {
2278 if (errCode == kWarning) {
2279 //Printf("Warning: Questionable detector response for track, most likely due to very low number of PID clusters! Debug output (dEdx_expected, sigma_expected):");
2280 }
2281 else
2282 Printf("Error: Failed to generate detector response for track - skipped! Debug output (dEdx_expected, sigma_expected):");
2283
2284 /*
2285 Printf("Pr: %e, %e", dEdxPr, sigmaPr);
2286 Printf("Pi: %e, %e", dEdxPi, sigmaPi);
2287 Printf("El: %e, %e", dEdxEl, sigmaEl);
2288 Printf("Mu: %e, %e", dEdxMu, sigmaMu);
2289 Printf("Ka: %e, %e", dEdxKa, sigmaKa);
2290 Printf("track: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), track->GetTPCmomentum(), track->Eta(),
2291 track->GetTPCsignalN());
2292 */
2293
2294 if (errCode != kWarning) {
2295 fhSkippedTracksForSignalGeneration->Fill(track->GetTPCmomentum(), track->GetTPCsignalN());
2296 return kFALSE;// Skip generated response in case of error
2297 }
2298 }
2299
2300 for (Int_t n = 0; n < numGenEntries; n++) {
2301 if (!isMC || !fUseMCidForGeneration) {
2302 // If no MC info is available or shall not be used, use weighting with priors to generate entries for the different species
2303 Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
2304
2305 // Consider generated response as originating from...
2306 if (rnd <= prob[AliPID::kElectron])
2307 genEntry[kGenMCID] = 0; // ... an electron
2308 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon])
2309 genEntry[kGenMCID] = 1; // ... a kaon
2310 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon])
2311 genEntry[kGenMCID] = 2; // ... a muon -> NOTE: prob[AliPID::kMuon] = 0 in case of fTakeIntoAccountMuons = kFALSE
2312 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon] + prob[AliPID::kPion])
2313 genEntry[kGenMCID] = 3; // ... a pion
2314 else
2315 genEntry[kGenMCID] = 4; // ... a proton
2316 }
2317
2318 // Electrons
2319 genEntry[kGenSelectSpecies] = 0;
2320 //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaEl[n];
2321 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeEl[n];
2322 fhGenEl->Fill(genEntry);
2323
2324 genEntry[kGenSelectSpecies] = 1;
2325 //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaKa[n];
2326 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeKa[n];
2327 fhGenEl->Fill(genEntry);
2328
2329 genEntry[kGenSelectSpecies] = 2;
2330 //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaPi[n];
2331 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePi[n];
2332 fhGenEl->Fill(genEntry);
2333
2334 genEntry[kGenSelectSpecies] = 3;
2335 //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaPr[n];
2336 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePr[n];
2337 fhGenEl->Fill(genEntry);
2338
2339 // Kaons
2340 genEntry[kGenSelectSpecies] = 0;
2341 //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaEl[n];
2342 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeEl[n];
2343 fhGenKa->Fill(genEntry);
2344
2345 genEntry[kGenSelectSpecies] = 1;
2346 //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaKa[n];
2347 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeKa[n];
2348 fhGenKa->Fill(genEntry);
2349
2350 genEntry[kGenSelectSpecies] = 2;
2351 //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaPi[n];
2352 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePi[n];
2353 fhGenKa->Fill(genEntry);
2354
2355 genEntry[kGenSelectSpecies] = 3;
2356 //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaPr[n];
2357 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePr[n];
2358 fhGenKa->Fill(genEntry);
2359
2360 // Pions
2361 genEntry[kGenSelectSpecies] = 0;
2362 //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaEl[n];
2363 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeEl[n];
2364 fhGenPi->Fill(genEntry);
2365
2366 genEntry[kGenSelectSpecies] = 1;
2367 //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaKa[n];
2368 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeKa[n];
2369 fhGenPi->Fill(genEntry);
2370
2371 genEntry[kGenSelectSpecies] = 2;
2372 //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaPi[n];
2373 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePi[n];
2374 fhGenPi->Fill(genEntry);
2375
2376 genEntry[kGenSelectSpecies] = 3;
2377 //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaPr[n];
2378 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePr[n];
2379 fhGenPi->Fill(genEntry);
2380
2381 if (fTakeIntoAccountMuons) {
2382 // Muons, if desired
2383 genEntry[kGenSelectSpecies] = 0;
2384 //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaEl[n];
2385 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeEl[n];
2386 fhGenMu->Fill(genEntry);
2387
2388 genEntry[kGenSelectSpecies] = 1;
2389 //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaKa[n];
2390 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeKa[n];
2391 fhGenMu->Fill(genEntry);
2392
2393 genEntry[kGenSelectSpecies] = 2;
2394 //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaPi[n];
2395 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePi[n];
2396 fhGenMu->Fill(genEntry);
2397
2398 genEntry[kGenSelectSpecies] = 3;
2399 //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaPr[n];
2400 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePr[n];
2401 fhGenMu->Fill(genEntry);
2402 }
2403
2404 // Protons
2405 genEntry[kGenSelectSpecies] = 0;
2406 //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaEl[n];
2407 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeEl[n];
2408 fhGenPr->Fill(genEntry);
2409
2410 genEntry[kGenSelectSpecies] = 1;
2411 //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaKa[n];
2412 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeKa[n];
2413 fhGenPr->Fill(genEntry);
2414
2415 genEntry[kGenSelectSpecies] = 2;
2416 //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaPi[n];
2417 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePi[n];
2418 fhGenPr->Fill(genEntry);
2419
2420 genEntry[kGenSelectSpecies] = 3;
2421 //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaPr[n];
2422 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePr[n];
2423 fhGenPr->Fill(genEntry);
2424 }
2425
2426 return kTRUE;
2427}
2428
2429
2430//_____________________________________________________________________________
2431Bool_t AliAnalysisTaskPID::SetConvolutedGaussLambdaParameter(Double_t lambda)
2432{
2433 // Set the lambda parameter of the convolution to the desired value. Automatically
2434 // calculates the parameters for the transition (restricted) gauss -> convoluted gauss.
2435 fConvolutedGaussTransitionPars[2] = lambda;
2436
2437 // Save old parameters and settings of function to restore them later:
2438 Double_t* oldFuncParams = new Double_t[fkConvolutedGausNPar];
2439 for (Int_t i = 0; i < fkConvolutedGausNPar; i++)
2440 oldFuncParams[i] = fConvolutedGausDeltaPrime->GetParameter(i);
2441 Int_t oldFuncNpx = fConvolutedGausDeltaPrime->GetNpx();
2442 Double_t oldFuncRangeLow = 0, oldFuncRangeUp = 100;
2443 fConvolutedGausDeltaPrime->GetRange(oldFuncRangeLow, oldFuncRangeUp);
2444
2445 // Choose some sufficiently large range
2446 const Double_t rangeStart = 0.5;
2447 const Double_t rangeEnd = 2.0;
2448
2449 // To get the parameters for the transition, just choose arbitrarily input parameters for mu and sigma
2450 // (it makes sense to choose typical values). The ratio sigma_gauss / sigma_convolution is independent
2451 // of mu and as well the difference mu_gauss - mu_convolution.
2452 Double_t muInput = 1.0;
2453 Double_t sigmaInput = fgkSigmaReferenceForTransitionPars;
2454
2455
2456 // Step 1: Generate distribution with input parameters
2457 const Int_t nPar = 3;
2458 Double_t inputPar[nPar] = { muInput, sigmaInput, lambda };
2459
2460 TH1D* hInput = new TH1D("hInput", "Input distribution", 2000, rangeStart, rangeEnd);
2461
2462 fConvolutedGausDeltaPrime->SetParameters(inputPar);
2463 fConvolutedGausDeltaPrime->SetRange(rangeStart, rangeEnd);
2464 fConvolutedGausDeltaPrime->SetNpx(2000);
2465
2466 /*OLD
2467 // The resolution and mean of the AliPIDResponse are determined in finite intervals
2468 // of ncl (also second order effects due to finite dEdx and tanTheta).
2469 // Take this into account for the transition parameters via assuming an approximately flat
2470 // distritubtion in ncl in this interval.
2471 // NOTE: The ncl interval should be the same as the one used for the sigma map creation!
2472 const Int_t nclStart = 151;
2473 const Int_t nclEnd = 160;
2474 const Int_t nclSteps = (nclEnd - nclStart) + 1;
2475 for (Int_t ncl = nclStart; ncl <= nclEnd; ncl++) {
2476 // Resolution scales with 1/sqrt(ncl)
2477 fConvolutedGausDeltaPrime->SetParameter(1, inputPar[1] * sqrt(nclEnd) / sqrt(ncl));
2478 TH1* hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
2479
2480 for (Int_t i = 0; i < 50000000 / nclSteps; i++)
2481 hInput->Fill(hProbDensity->GetRandom());
2482 }
2483 */
2484
2485 TH1* hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
2486
2487 for (Int_t i = 0; i < 50000000; i++)
2488 hInput->Fill(hProbDensity->GetRandom());
2489
2490 // Step 2: Fit generated distribution with restricted gaussian
2491 Int_t maxBin = hInput->GetMaximumBin();
2492 Double_t max = hInput->GetBinContent(maxBin);
2493
2494 UChar_t usedBins = 1;
2495 if (maxBin > 1) {
2496 max += hInput->GetBinContent(maxBin - 1);
2497 usedBins++;
2498 }
2499 if (maxBin < hInput->GetNbinsX()) {
2500 max += hInput->GetBinContent(maxBin + 1);
2501 usedBins++;
2502 }
2503 max /= usedBins;
2504
2505 // NOTE: The range (<-> fraction of maximum) should be the same
2506 // as for the spline and eta maps creation
2507 const Double_t lowThreshold = hInput->GetXaxis()->GetBinLowEdge(hInput->FindFirstBinAbove(0.1 * max));
2508 const Double_t highThreshold = hInput->GetXaxis()->GetBinUpEdge(hInput->FindLastBinAbove(0.1 * max));
2509
2510 TFitResultPtr fitResGaussFirstStep = hInput->Fit("gaus", "QNRS", "", lowThreshold, highThreshold);
2511
2512 TFitResultPtr fitResGauss;
2513
2514 if ((Int_t)fitResGaussFirstStep == 0) {
2515 TF1 fGauss("fGauss", "[0]*TMath::Gaus(x, [1], [2], 1)", rangeStart, rangeEnd);
2516 fGauss.SetParameter(0, fitResGaussFirstStep->GetParams()[0]);
2517 fGauss.SetParError(0, fitResGaussFirstStep->GetErrors()[0]);
2518 fGauss.SetParameter(2, fitResGaussFirstStep->GetParams()[2]);
2519 fGauss.SetParError(2, fitResGaussFirstStep->GetErrors()[2]);
2520
2521 fGauss.FixParameter(1, fitResGaussFirstStep->GetParams()[1]);
2522 fitResGauss = hInput->Fit(&fGauss, "QNS", "s", rangeStart, rangeEnd);
2523 }
2524 else {
2525 fitResGauss = hInput->Fit("gaus", "QNRS", "same", rangeStart, rangeEnd);
2526 }
2527 //OLD TFitResultPtr fitResGauss = hInput->Fit("gaus", "QNRS", "", hInput->GetXaxis()->GetBinLowEdge(hInput->FindFirstBinAbove(0.1 * max)),
2528 // hInput->GetXaxis()->GetBinUpEdge(hInput->FindLastBinAbove(0.1 * max)));
2529
2530
2531 // Step 3: Use parameters from gaussian fit to obtain parameters for the transition "restricted gauss" -> "convoluted gauss"
2532
2533 // 3.1 The ratio sigmaInput / sigma_gaussFit ONLY depends on lambda (which is fixed per period) -> Calculate this first
2534 // for an arbitrary (here: typical) sigma. The ratio is then ~the same for ALL sigma for given lambda!
2535 if ((Int_t)fitResGauss != 0) {
2536 AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Gauss Fit failed!\n");
2537 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2538 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2539 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2540
2541 delete hInput;
2542 delete oldFuncParams;
2543
2544 return kFALSE;
2545 }
2546
2547 if (fitResGauss->GetParams()[2] <= 0) {
2548 AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Sigma of gauss fit <= 0!\n");
2549 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2550 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2551 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2552
2553 delete hInput;
2554 delete oldFuncParams;
2555
2556 return kFALSE;
2557 }
2558
2559 // sigma correction factor
2560 fConvolutedGaussTransitionPars[1] = sigmaInput / fitResGauss->GetParams()[2];
2561
2562 // 3.2 Now that sigma und lambda are determined, one can calculate mu by shifting the maximum to the desired position,
2563 // i.e. the maximum of the input distribution should coincide with that of the re-calculated distribution,
2564 // which is achieved by getting the same mu for the same sigma.
2565 // NOTE: For fixed lambda, the shift is proportional to sigma and independent of mu!
2566 // So, one can calculate the shift for an arbitrary fixed (here: typical)
2567 // sigma and then simply use this shift for any other sigma by scaling it correspondingly!!!
2568
2569 // Mu shift correction:
2570 // Shift in mu (difference between mean of gauss and mean of convolution) is proportional to sigma!
2571 // Thus, choose a reference sigma (typical -> 0.05), such that for arbitrary sigma one can simple scale
2572 // this shift correction with sigma / referenceSigma.
2573 fConvolutedGaussTransitionPars[0] = (fitResGauss->GetParams()[1] - muInput);
2574
2575
2576 /*Changed on 03.07.2013
2577
2578 // Maximum of fConvolutedGausDeltaPrime should agree with maximum of input
2579 Double_t par[nPar] = { fitResGauss->GetParams()[1], // just as a guess of the maximum position
2580 sigmaInput,
2581 lambda };
2582
2583 fConvolutedGausDeltaPrime->SetParameters(par);
2584
2585 Double_t maxXInput = fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001, muInput - 3. * sigmaInput),
2586 muInput + 10. * sigmaInput,
2587 0.0001);
2588
2589 // Maximum shifts proportional to sigma and is linear in mu (~mean of gauss)
2590 // Maximum should be typically located within [gaussMean, gaussMean + 3 gaussSigma].
2591 // -> Larger search range for safety reasons (also: sigma and/or mean might be not 100% accurate).
2592 Double_t maxXconvoluted = fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001,
2593 fitResGauss->GetParams()[1] - 3. * fitResGauss->GetParams()[2]),
2594 fitResGauss->GetParams()[1] + 10. * fitResGauss->GetParams()[2],
2595 0.0001);
2596 if (maxXconvoluted <= 0) {
2597 AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Maximum of fConvolutedGausDeltaPrime <= 0!\n");
2598
2599 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2600 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2601 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2602
2603 delete hInput;
2604 delete oldFuncParams;
2605
2606 return kFALSE;
2607 }
2608
2609 // maxX perfectly shifts as par[0] (scaled by sigma) -> Can shift maxX to input value.
2610 // Mu shift correction:
2611 fConvolutedGaussTransitionPars[0] = maxXconvoluted - maxXInput;
2612 */
2613
2614
2615
2616 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2617 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2618 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2619
2620 delete hInput;
2621 delete oldFuncParams;
2622
2623 return kTRUE;
2624}
2625
2626
2627//_____________________________________________________________________________
2628AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::SetParamsForConvolutedGaus(const Double_t gausMean, const Double_t gausSigma)
2629{
2630 // Set parameters for convoluted gauss using parameters for a pure gaussian.
2631 // If SetConvolutedGaussLambdaParameter has not been called before to initialise the translation parameters,
2632 // some default parameters will be used and an error will show up.
2633
2634 if (fConvolutedGaussTransitionPars[1] < -998) {
2635 AliError("Transition parameters not initialised! Default parameters will be used. Please call SetConvolutedGaussLambdaParameter(...) before any calculations!");
2636 SetConvolutedGaussLambdaParameter(2.0);
2637 AliError(Form("Parameters set to:\n[0]: %f\n[1]: %f\n[2]: %f\n", fConvolutedGaussTransitionPars[0],
2638 fConvolutedGaussTransitionPars[1], fConvolutedGaussTransitionPars[2]));
2639 }
2640
2641 Double_t par[fkConvolutedGausNPar];
2642 par[2] = fConvolutedGaussTransitionPars[2];
2643 par[1] = fConvolutedGaussTransitionPars[1] * gausSigma;
2644 // maxX perfectly shifts as par[0] (scaled by sigma) -> Can shift maxX so that it sits at the right place.
2645 par[0] = gausMean - fConvolutedGaussTransitionPars[0] * par[1] / fgkSigmaReferenceForTransitionPars;
2646
2647 ErrorCode errCode = kNoErrors;
2648 fConvolutedGausDeltaPrime->SetParameters(par);
2649
2650 fConvolutedGausDeltaPrime->SetNpx(20); // Small value speeds up following algorithm (valid, since extrema far apart)
2651
2652 // Accuracy of 10^-5 is enough to get 0.1% precise peak for MIPS w.r.t. to dEdx = 2000 of protons
2653 // (should boost up the algorithm, because 10^-10 is the default value!)
2654 Double_t maxX= fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001, gausMean - 2. * gausSigma),
2655 gausMean + 6. * gausSigma, 1.0E-5);
2656
2657 const Double_t maximum = fConvolutedGausDeltaPrime->Eval(maxX);
2658 const Double_t maximumFraction = maximum * fAccuracyNonGaussianTail;
2659
2660 // Estimate lower boundary for subsequent search:
2661 Double_t lowBoundSearchBoundLow = TMath::Max(1e-4, maxX - 5. * gausSigma);
2662 Double_t lowBoundSearchBoundUp = maxX;
2663
2664 Bool_t lowerBoundaryFixedAtZero = kFALSE;
2665
2666 while (fConvolutedGausDeltaPrime->Eval(lowBoundSearchBoundLow) >= maximumFraction) {
2667 if (lowBoundSearchBoundLow <= 0) {
2668 // This should only happen to low dEdx particles with very few clusters and therefore large sigma, such that the gauss goes below zero deltaPrime
2669 if (maximum <= 0) { // Something is weired
2670 printf("Error generating signal: maximum is <= 0!\n");
2671 return kError;
2672 }
2673 else {
2674 const Double_t valueAtZero = fConvolutedGausDeltaPrime->Eval(0);
2675 if (valueAtZero / maximum > 0.05) {
2676 // Too large fraction below zero deltaPrime. Signal generation cannot be reliable in this case
2677 printf("Error generating signal: Too large fraction below zero deltaPrime: convGauss(0) / convGauss(max) = %e / %e = %e!\n",
2678 valueAtZero, maximum, valueAtZero / maximum);
2679 return kError;
2680 }
2681 }
2682
2683 /*
2684 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",
2685 fConvolutedGausDeltaPrime->Eval(maxX), fAccuracyNonGaussianTail, maximumFraction, maxX, par[0], par[1], par[2], gausMean, gausSigma);
2686 */
2687
2688 lowerBoundaryFixedAtZero = kTRUE;
2689
2690 if (errCode != kError)
2691 errCode = kWarning;
2692
2693 break;
2694 }
2695
2696 lowBoundSearchBoundUp -= gausSigma;
2697 lowBoundSearchBoundLow -= gausSigma;
2698
2699 if (lowBoundSearchBoundLow < 0) {
2700 lowBoundSearchBoundLow = 0;
2701 lowBoundSearchBoundUp += gausSigma;
2702 }
2703 }
2704
2705 // Determine lower boundary inside estimated range. For small values of the maximum: Need more precision, since finer binning!
2706 Double_t rangeStart = lowerBoundaryFixedAtZero ? 0 :
2707 fConvolutedGausDeltaPrime->GetX(maximumFraction, lowBoundSearchBoundLow, lowBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
2708
2709 // .. and the same for the upper boundary
2710 Double_t rangeEnd = 0;
2711 // If distribution starts beyond upper boundary, everything ends up in the overflow bin. So, just reduce range and Npx to minimum
2712 if (rangeStart > fkDeltaPrimeUpLimit) {
2713 rangeEnd = rangeStart + 0.00001;
2714 fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
2715 fConvolutedGausDeltaPrime->SetNpx(4);
2716 }
2717 else {
2718 // Estimate upper boundary for subsequent search:
2719 Double_t upBoundSearchBoundUp = maxX + 5 * gausSigma;
2720 Double_t upBoundSearchBoundLow = maxX;
2721 while (fConvolutedGausDeltaPrime->Eval(upBoundSearchBoundUp) >= maximumFraction) {
2722 upBoundSearchBoundUp += gausSigma;
2723 upBoundSearchBoundLow += gausSigma;
2724 }
2725
2726 // For small values of the maximum: Need more precision, since finer binning!
2727 rangeEnd = fConvolutedGausDeltaPrime->GetX(maximumFraction, upBoundSearchBoundLow, upBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
2728
2729 fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
2730 fConvolutedGausDeltaPrime->SetNpx(fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindBin(rangeEnd)
2731 - fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindBin(rangeStart) + 1);
2732 //fConvolutedGausDeltaPrime->SetNpx((rangeEnd - rangeStart) / fDeltaPrimeBinWidth + 1);
2733 }
2734
2735 return errCode;
2736}
2737
2738
2739//________________________________________________________________________
2740void AliAnalysisTaskPID::SetUpGenHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const
2741{
2742 // Sets bin limits for axes which are not standard binned and the axes titles.
2743
2744 hist->SetBinEdges(kGenPt, binsPt);
2745 hist->SetBinEdges(kGenDeltaPrimeSpecies, binsDeltaPrime);
2746 hist->SetBinEdges(kGenCentrality, binsCent);
2747
2748 if (fStoreAdditionalJetInformation)
2749 hist->SetBinEdges(kGenJetPt, binsJetPt);
2750
2751 // Set axes titles
2752 hist->GetAxis(kGenMCID)->SetTitle("MC PID");
2753 hist->GetAxis(kGenMCID)->SetBinLabel(1, "e");
2754 hist->GetAxis(kGenMCID)->SetBinLabel(2, "K");
2755 hist->GetAxis(kGenMCID)->SetBinLabel(3, "#mu");
2756 hist->GetAxis(kGenMCID)->SetBinLabel(4, "#pi");
2757 hist->GetAxis(kGenMCID)->SetBinLabel(5, "p");
2758
2759 hist->GetAxis(kGenSelectSpecies)->SetTitle("Select Species");
2760 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(1, "e");
2761 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(2, "K");
2762 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(3, "#pi");
2763 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(4, "p");
2764
2765 hist->GetAxis(kGenPt)->SetTitle("P_{T} (GeV/c)");
2766
2767 hist->GetAxis(kGenDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
2768
2769 hist->GetAxis(kGenCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
2770
2771 if (fStoreAdditionalJetInformation) {
2772 hist->GetAxis(kGenJetPt)->SetTitle("P_{T}^{jet} (GeV/c)");
2773
2774 hist->GetAxis(kGenZ)->SetTitle("z = P_{T}^{track} / P_{T}^{jet}");
2775
2776 hist->GetAxis(kGenXi)->SetTitle("#xi = ln(P_{T}^{jet} / P_{T}^{track})");
2777 }
2778
2779 hist->GetAxis(GetIndexOfChargeAxisGen())->SetTitle("Charge (e_{0})");
2780}
2781
2782
2783//________________________________________________________________________
2784void AliAnalysisTaskPID::SetUpGenYieldHist(THnSparse* hist, Double_t* binsPt, Double_t* binsCent, Double_t* binsJetPt) const
2785{
2786 // Sets bin limits for axes which are not standard binned and the axes titles.
2787
2788 hist->SetBinEdges(kGenYieldPt, binsPt);
2789 hist->SetBinEdges(kGenYieldCentrality, binsCent);
2790 if (fStoreAdditionalJetInformation)
2791 hist->SetBinEdges(kGenYieldJetPt, binsJetPt);
2792
2793 for (Int_t i = 0; i < 5; i++)
2794 hist->GetAxis(kGenYieldMCID)->SetBinLabel(i + 1, AliPID::ParticleLatexName(i));
2795
2796 // Set axes titles
2797 hist->GetAxis(kGenYieldMCID)->SetTitle("MC PID");
2798 hist->GetAxis(kGenYieldPt)->SetTitle("P_{T}^{gen} (GeV/c)");
2799 hist->GetAxis(kGenYieldCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
2800
2801 if (fStoreAdditionalJetInformation) {
2802 hist->GetAxis(kGenYieldJetPt)->SetTitle("P_{T}^{jet, gen} (GeV/c)");
2803
2804 hist->GetAxis(kGenYieldZ)->SetTitle("z = P_{T}^{track} / P_{T}^{jet}");
2805
2806 hist->GetAxis(kGenYieldXi)->SetTitle("#xi = ln(P_{T}^{jet} / P_{T}^{track})");
2807 }
2808
2809 hist->GetAxis(GetIndexOfChargeAxisGenYield())->SetTitle("Charge (e_{0})");
2810}
2811
2812
2813//________________________________________________________________________
2814void AliAnalysisTaskPID::SetUpHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const
2815{
2816 // Sets bin limits for axes which are not standard binned and the axes titles.
2817
2818 hist->SetBinEdges(kDataPt, binsPt);
2819 hist->SetBinEdges(kDataDeltaPrimeSpecies, binsDeltaPrime);
2820 hist->SetBinEdges(kDataCentrality, binsCent);
2821
2822 if (fStoreAdditionalJetInformation)
2823 hist->SetBinEdges(kDataJetPt, binsJetPt);
2824
2825 // Set axes titles
2826 hist->GetAxis(kDataMCID)->SetTitle("MC PID");
2827 hist->GetAxis(kDataMCID)->SetBinLabel(1, "e");
2828 hist->GetAxis(kDataMCID)->SetBinLabel(2, "K");
2829 hist->GetAxis(kDataMCID)->SetBinLabel(3, "#mu");
2830 hist->GetAxis(kDataMCID)->SetBinLabel(4, "#pi");
2831 hist->GetAxis(kDataMCID)->SetBinLabel(5, "p");
2832
2833 hist->GetAxis(kDataSelectSpecies)->SetTitle("Select Species");
2834 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(1, "e");
2835 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(2, "K");
2836 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(3, "#pi");
2837 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(4, "p");
2838
2839 hist->GetAxis(kDataPt)->SetTitle("P_{T} (GeV/c)");
2840
2841 hist->GetAxis(kDataDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
2842
2843 hist->GetAxis(kDataCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
2844
2845 if (fStoreAdditionalJetInformation) {
2846 hist->GetAxis(kDataJetPt)->SetTitle("P_{T}^{jet} (GeV/c)");
2847
2848 hist->GetAxis(kDataZ)->SetTitle("z = P_{T}^{track} / P_{T}^{jet}");
2849
2850 hist->GetAxis(kDataXi)->SetTitle("#xi = ln(P_{T}^{jet} / P_{T}^{track})");
2851 }
2852
2853 hist->GetAxis(GetIndexOfChargeAxisData())->SetTitle("Charge (e_{0})");
2854
2855 /*OLD with TOF, p_TPC_Inner and p_vertex
2856 // MC PID, SelectSpecies, P(TPC_inner), pT, p(Vertex), deltaSpecies, deltaPrimeSpecies, deltaTOFspecies
2857 hist->SetBinEdges(2, binsPt);
2858 hist->SetBinEdges(3, binsPt);
2859 hist->SetBinEdges(4, binsPt);
2860
2861 // Set axes titles
2862 hist->GetAxis(0)->SetTitle("MC PID");
2863 hist->GetAxis(0)->SetBinLabel(1, "e");
2864 hist->GetAxis(0)->SetBinLabel(2, "K");
2865 hist->GetAxis(0)->SetBinLabel(3, "#mu");
2866 hist->GetAxis(0)->SetBinLabel(4, "#pi");
2867 hist->GetAxis(0)->SetBinLabel(5, "p");
2868
2869 hist->GetAxis(1)->SetTitle("Select Species");
2870 hist->GetAxis(1)->SetBinLabel(1, "e");
2871 hist->GetAxis(1)->SetBinLabel(2, "K");
2872 hist->GetAxis(1)->SetBinLabel(3, "#pi");
2873 hist->GetAxis(1)->SetBinLabel(4, "p");
2874
2875 hist->GetAxis(2)->SetTitle("P_{TPC_inner} (GeV/c)");
2876 hist->GetAxis(3)->SetTitle("P_{T} (GeV/c)");
2877 hist->GetAxis(4)->SetTitle("P_{vertex} (GeV/c)");
2878
2879 hist->GetAxis(5)->SetTitle("TPC #Delta_{species} (arb. unit)");
2880
2881 hist->GetAxis(6)->SetTitle("TPC #Delta'_{species} (arb. unit)");
2882
2883 hist->GetAxis(7)->SetTitle("#Delta TOF_{species} (ps)");
2884 */
2885}