]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/UserTasks/AliAnalysisTaskPID.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Add task to handle the muon type in simulations. The class allows to identify, among...
[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//_____________________________________________________________________________
57fa828c 1024Int_t AliAnalysisTaskPID::FindFirstBinAboveIn3dSubset(const TH3* hist, const Double_t threshold, const Int_t yBin,
1025 const Int_t zBin) const
e131b05f 1026{
1027 // Kind of projects a TH3 to 1 bin combination in y and z
1028 // and looks for the first x bin above a threshold for this projection.
1029 // If no such bin is found, -1 is returned.
1030
1031 if (!hist)
1032 return -1;
1033
1034 Int_t nBinsX = hist->GetNbinsX();
1035 for (Int_t xBin = 1; xBin <= nBinsX; xBin++) {
1036 if (hist->GetBinContent(xBin, yBin, zBin) > threshold)
1037 return xBin;
1038 }
1039
1040 return -1;
1041}
1042
1043
1044//_____________________________________________________________________________
57fa828c 1045Int_t AliAnalysisTaskPID::FindLastBinAboveIn3dSubset(const TH3* hist, const Double_t threshold, const Int_t yBin,
1046 const Int_t zBin) const
e131b05f 1047{
1048 // Kind of projects a TH3 to 1 bin combination in y and z
1049 // and looks for the last x bin above a threshold for this projection.
1050 // If no such bin is found, -1 is returned.
1051
1052 if (!hist)
1053 return -1;
1054
1055 Int_t nBinsX = hist->GetNbinsX();
1056 for (Int_t xBin = nBinsX; xBin >= 1; xBin--) {
1057 if (hist->GetBinContent(xBin, yBin, zBin) > threshold)
1058 return xBin;
1059 }
1060
1061 return -1;
1062}
1063
1064
1065//_____________________________________________________________________________
1066Bool_t AliAnalysisTaskPID::GetParticleFraction(const Double_t trackPt, const Double_t jetPt, const Double_t centralityPercentile,
1067 const AliPID::EParticleType species,
1068 Double_t& fraction, Double_t& fractionErrorStat, Double_t& fractionErrorSys) const
1069{
1070 // Computes the particle fraction for the corresponding species for the given trackPt, jetPt and centrality.
1071 // Use jetPt = -1 for inclusive spectra and centralityPercentile = -1 for pp.
1072 // On success (return value kTRUE), fraction contains the particle fraction, fractionErrorStat(Sys) the sigma of its
1073 // statistical (systematic) error
1074
1075 fraction = -999.;
1076 fractionErrorStat = 999.;
1077 fractionErrorSys = 999.;
1078
1079 if (species > AliPID::kProton || species < AliPID::kElectron) {
1080 AliError(Form("Only fractions for species index %d to %d availabe, but not for the requested one: %d", 0, AliPID::kProton, species));
1081 return kFALSE;
1082 }
1083
1084 if (!fFractionHists[species]) {
1085 AliError(Form("Histo with particle fractions for species %d not loaded!", species));
1086
1087 return kFALSE;
1088 }
1089
1090 Int_t jetPtBin = FindBinWithinRange(fFractionHists[species]->GetYaxis(), jetPt);
1091 Int_t centBin = FindBinWithinRange(fFractionHists[species]->GetZaxis(), centralityPercentile);
1092
1093 // The following interpolation takes the bin content of the first/last available bin,
1094 // if requested point lies beyond bin center of first/last bin.
1095 // The interpolation is only done for the x-axis (track pT), i.e. jetPtBin and centBin are fix,
1096 // because the analysis will anyhow be run in bins of jetPt and centrality and
1097 // it is not desired to correlate different jetPt bins via interpolation.
1098
1099 // The same procedure is used for the error of the fraction
1100 TAxis* xAxis = fFractionHists[species]->GetXaxis();
1101
1102 // No interpolation to values beyond the centers of the first/last bins (we don't know exactly where the spectra start or stop,
1103 // thus, search for the first and last bin above 0.0 to constrain the range
1104 Int_t firstBin = TMath::Max(1, FindFirstBinAboveIn3dSubset(fFractionHists[species], 0.0, jetPtBin, centBin));
1105 Int_t lastBin = TMath::Min(fFractionHists[species]->GetNbinsX(),
1106 FindLastBinAboveIn3dSubset(fFractionHists[species], 0.0, jetPtBin, centBin));
1107
1108 if (trackPt <= xAxis->GetBinCenter(firstBin)) {
1109 fraction = fFractionHists[species]->GetBinContent(firstBin, jetPtBin, centBin);
1110 fractionErrorStat = fFractionHists[species]->GetBinError(firstBin, jetPtBin, centBin);
1111 fractionErrorSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(firstBin, jetPtBin, centBin) : 0.;
1112 }
1113 else if (trackPt >= xAxis->GetBinCenter(lastBin)) {
1114 fraction = fFractionHists[species]->GetBinContent(lastBin, jetPtBin, centBin);
1115 fractionErrorStat = fFractionHists[species]->GetBinError(lastBin, jetPtBin, centBin);
1116 fractionErrorSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(lastBin, jetPtBin, centBin) : 0.;
1117 }
1118 else {
1119 Double_t x0 = 0., x1 = 0., y0 = 0., y1 = 0.;
1120 Double_t y0errStat = 0., y1errStat = 0., y0errSys = 0., y1errSys = 0.;
1121 Int_t trackPtBin = xAxis->FindBin(trackPt);
1122
1123 // Linear interpolation between nearest neighbours in trackPt
1124 if (trackPt <= xAxis->GetBinCenter(trackPtBin)) {
1125 y0 = fFractionHists[species]->GetBinContent(trackPtBin - 1, jetPtBin, centBin);
1126 y0errStat = fFractionHists[species]->GetBinError(trackPtBin - 1, jetPtBin, centBin);
1127 y0errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin - 1, jetPtBin, centBin)
1128 : 0.;
1129 x0 = xAxis->GetBinCenter(trackPtBin - 1);
1130 y1 = fFractionHists[species]->GetBinContent(trackPtBin, jetPtBin, centBin);
1131 y1errStat = fFractionHists[species]->GetBinError(trackPtBin, jetPtBin, centBin);
1132 y1errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin, jetPtBin, centBin)
1133 : 0.;
1134 x1 = xAxis->GetBinCenter(trackPtBin);
1135 }
1136 else {
1137 y0 = fFractionHists[species]->GetBinContent(trackPtBin, jetPtBin, centBin);
1138 y0errStat = fFractionHists[species]->GetBinError(trackPtBin, jetPtBin, centBin);
1139 y0errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin, jetPtBin, centBin)
1140 : 0.;
1141 x0 = xAxis->GetBinCenter(trackPtBin);
1142 y1 = fFractionHists[species]->GetBinContent(trackPtBin + 1, jetPtBin, centBin);
1143 y1errStat = fFractionHists[species]->GetBinError(trackPtBin + 1, jetPtBin, centBin);
1144 y1errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin + 1, jetPtBin, centBin)
1145 : 0.;
1146 x1 = xAxis->GetBinCenter(trackPtBin + 1);
1147 }
1148
1149 // Per construction: x0 < trackPt < x1
1150 fraction = y0 + (trackPt - x0) * ((y1 - y0) / (x1 - x0));
1151 fractionErrorStat = y0errStat + (trackPt - x0) * ((y1errStat - y0errStat) / (x1 - x0));
1152 fractionErrorSys = fFractionSysErrorHists[species] ? (y0errSys + (trackPt - x0) * ((y1errSys - y0errSys) / (x1 - x0))) : 0.;
1153 }
1154
1155 return kTRUE;
1156}
1157
1158
1159//_____________________________________________________________________________
1160Bool_t AliAnalysisTaskPID::GetParticleFractions(const Double_t trackPt, const Double_t jetPt, const Double_t centralityPercentile,
1161 Double_t* prob, const Int_t smearSpeciesByError,
1162 const Int_t takeIntoAccountSpeciesSysError, const Bool_t uniformSystematicError) const
1163{
1164 // Fills the particle fractions for the given trackPt, jetPt and centrality into "prob".
1165 // Use jetPt = -1 for inclusive spectra and centralityPercentile = -1 for pp.
1166 // If smearSpeciesByError is >= 0 && < AliPID::kSPECIES, the returned fractions will be a random number distributed
1167 // with a gauss with mean being the corresponding particle fraction and sigma it's error for the considered species
1168 // "smearSpeciesByError".
1169 // Note that in this case the fractions for all species will NOT sum up to 1!
1170 // Thus, all other species fractions will be re-scaled weighted with their corresponding statistical error.
1171 // A similar procedure is used for "takeIntoAccountSpeciesSysError": The systematic error of the corresponding species
1172 // is used to generate a random number with uniform distribution in [mean - sysError, mean + sysError] for the new mean
1173 // (in cace of uniformSystematicError = kTRUE, otherwise it will be a gaus(mean, sysError)),
1174 // then the other species will be re-scaled according to their systematic errors.
1175 // First, the systematic error uncertainty procedure will be performed (that is including re-scaling), then the statistical
1176 // uncertainty procedure.
1177 // On success, kTRUE is returned.
1178
1179 if (!prob || smearSpeciesByError >= AliPID::kSPECIES || takeIntoAccountSpeciesSysError >= AliPID::kSPECIES)
1180 return kFALSE;
1181
1182 Double_t probTemp[AliPID::kSPECIES];
1183 Double_t probErrorStat[AliPID::kSPECIES];
1184 Double_t probErrorSys[AliPID::kSPECIES];
1185
1186 Bool_t success = kTRUE;
1187 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kElectron,
1188 probTemp[AliPID::kElectron], probErrorStat[AliPID::kElectron],
1189 probErrorSys[AliPID::kElectron]);
1190 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kMuon,
1191 probTemp[AliPID::kMuon], probErrorStat[AliPID::kMuon], probErrorSys[AliPID::kMuon]);
1192 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kPion,
1193 probTemp[AliPID::kPion], probErrorStat[AliPID::kPion], probErrorSys[AliPID::kPion]);
1194 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kKaon,
1195 probTemp[AliPID::kKaon], probErrorStat[AliPID::kKaon], probErrorSys[AliPID::kKaon]);
1196 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kProton,
1197 probTemp[AliPID::kProton], probErrorStat[AliPID::kProton], probErrorSys[AliPID::kProton]);
1198
1199 if (!success)
1200 return kFALSE;
1201
1202 // If desired, take into account the systematic error of the corresponding species and re-generate probTemp accordingly
1203 if (takeIntoAccountSpeciesSysError >= 0) {
1204 // Generate random fraction of the considered species "smearSpeciesByError" according to mean and sigma
1205 Double_t generatedFraction = uniformSystematicError
1206 ? fRandom->Rndm() * 2. * probErrorSys[takeIntoAccountSpeciesSysError]
1207 - probErrorSys[takeIntoAccountSpeciesSysError]
1208 + probTemp[takeIntoAccountSpeciesSysError]
1209 : fRandom->Gaus(probTemp[takeIntoAccountSpeciesSysError],
1210 probErrorSys[takeIntoAccountSpeciesSysError]);
1211
1212 // Catch cases with invalid fraction (can happen for large errors), i.e. fraction < 0 or > 1
1213 if (generatedFraction < 0.)
1214 generatedFraction = 0.;
1215 else if (generatedFraction > 1.)
1216 generatedFraction = 1.;
1217
1218 // Calculate difference from original fraction (original fractions sum up to 1!)
1219 Double_t deltaFraction = generatedFraction - probTemp[takeIntoAccountSpeciesSysError];
1220
1221 // Fractions must (including errors) lie inside [0,1] -> Adapt weights accordingly by setting the errors
1222 if (deltaFraction > 0) {
1223 // Some part will be SUBTRACTED from the other fractions
1224 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1225 if (probTemp[i] - probErrorSys[i] < 0)
1226 probErrorSys[i] = probTemp[i];
1227 }
1228 }
1229 else {
1230 // Some part will be ADDED to the other fractions
1231 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1232 if (probTemp[i] + probErrorSys[i] > 1)
1233 probErrorSys[i] = 1. - probTemp[i];
1234 }
1235 }
1236
1237 // Compute summed weight of all fractions except for the considered one
1238 Double_t summedWeight = 0.;
1239 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1240 if (i != takeIntoAccountSpeciesSysError)
1241 summedWeight += probErrorSys[i];
1242 }
1243
1244 // Compute the weight for the other species
1245 /*
1246 if (summedWeight <= 1e-13) {
1247 // If this happens for some reason (it should not!), just assume flat weight
1248 printf("Error: summedWeight (sys error) ~ 0 for trackPt %f, jetPt %f, centralityPercentile %f. Setting flat weight!\n",
1249 trackPt, jetPt, centralityPercentile);
1250 }*/
1251
1252 Double_t weight[AliPID::kSPECIES];
1253 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1254 if (i != takeIntoAccountSpeciesSysError) {
1255 if (summedWeight > 1e-13)
1256 weight[i] = probErrorSys[i] / summedWeight;
1257 else
1258 weight[i] = probErrorSys[i] / (AliPID::kSPECIES - 1);
1259 }
1260 }
1261
1262 // For the final generated fractions, set the generated value for the considered species
1263 // and the generated value minus delta times statistical weight
1264 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1265 if (i != takeIntoAccountSpeciesSysError)
1266 probTemp[i] = probTemp[i] - weight[i] * deltaFraction;
1267 else
1268 probTemp[i] = generatedFraction;
1269 }
1270 }
1271
1272 // Using the values of probTemp (either the original ones or those after taking into account the systematic error),
1273 // calculate the final fractions - if the statistical error is to be taken into account, smear the corresponding
1274 // fraction. If not, just write probTemp to the final result array.
1275 if (smearSpeciesByError >= 0) {
1276 // Generate random fraction of the considered species "smearSpeciesByError" according to mean and sigma
1277 Double_t generatedFraction = fRandom->Gaus(probTemp[smearSpeciesByError], probErrorStat[smearSpeciesByError]);
1278
1279 // Catch cases with invalid fraction (can happen for large errors), i.e. fraction < 0 or > 1
1280 if (generatedFraction < 0.)
1281 generatedFraction = 0.;
1282 else if (generatedFraction > 1.)
1283 generatedFraction = 1.;
1284
1285 // Calculate difference from original fraction (original fractions sum up to 1!)
1286 Double_t deltaFraction = generatedFraction - probTemp[smearSpeciesByError];
1287
1288 // Fractions must (including errors) lie inside [0,1] -> Adapt weights accordingly by setting the errors
1289 if (deltaFraction > 0) {
1290 // Some part will be SUBTRACTED from the other fractions
1291 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1292 if (probTemp[i] - probErrorStat[i] < 0)
1293 probErrorStat[i] = probTemp[i];
1294 }
1295 }
1296 else {
1297 // Some part will be ADDED to the other fractions
1298 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1299 if (probTemp[i] + probErrorStat[i] > 1)
1300 probErrorStat[i] = 1. - probTemp[i];
1301 }
1302 }
1303
1304 // Compute summed weight of all fractions except for the considered one
1305 Double_t summedWeight = 0.;
1306 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1307 if (i != smearSpeciesByError)
1308 summedWeight += probErrorStat[i];
1309 }
1310
1311 // Compute the weight for the other species
1312 /*
1313 if (summedWeight <= 1e-13) {
1314 // If this happens for some reason (it should not!), just assume flat weight
1315 printf("Error: summedWeight (stat error) ~ 0 for trackPt %f, jetPt %f, centralityPercentile %f. Setting flat weight!\n",
1316 trackPt, jetPt, centralityPercentile);
1317 }*/
1318
1319 Double_t weight[AliPID::kSPECIES];
1320 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1321 if (i != smearSpeciesByError) {
1322 if (summedWeight > 1e-13)
1323 weight[i] = probErrorStat[i] / summedWeight;
1324 else
1325 weight[i] = probErrorStat[i] / (AliPID::kSPECIES - 1);
1326 }
1327 }
1328
1329 // For the final generated fractions, set the generated value for the considered species
1330 // and the generated value minus delta times statistical weight
1331 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1332 if (i != smearSpeciesByError)
1333 prob[i] = probTemp[i] - weight[i] * deltaFraction;
1334 else
1335 prob[i] = generatedFraction;
1336 }
1337 }
1338 else {
1339 // Just take the generated values
1340 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
1341 prob[i] = probTemp[i];
1342 }
1343
1344
1345 // Should already be normalised, but make sure that it really is:
1346 Double_t probSum = 0.;
1347 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1348 probSum += prob[i];
1349 }
1350
1351 if (probSum <= 0)
1352 return kFALSE;
1353
1354 if (TMath::Abs(probSum - 1.0) > 1e-4) {
1355 printf("Warning: Re-normalising sum of fractions: Sum is %e\n", probSum);
1356 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1357 prob[i] /= probSum;
1358 }
1359 }
1360
1361 return kTRUE;
1362}
1363
1364
1365//_____________________________________________________________________________
1366const TH3D* AliAnalysisTaskPID::GetParticleFractionHisto(const Int_t species, const Bool_t sysError) const
1367{
1368 if (species < AliPID::kElectron || species > AliPID::kProton)
1369 return 0x0;
1370
1371 return sysError ? fFractionSysErrorHists[species] : fFractionHists[species];
1372}
1373
1374
1375//_____________________________________________________________________________
1376Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(Int_t motherPDG, Double_t motherGenPt)
1377{
1378 // Strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
1379 // -> Based on function in PWGJE/AliAnalysisTaskFragmentationFunction, which uses
1380 // the following data from CMS pp @ 7 TeV inclusive (JHEP 05 (2011) 064)
1381
1382 Double_t fac = 1;
1383
1384 const Int_t absMotherPDG = TMath::Abs(motherPDG);
1385
1386 if (absMotherPDG == 310 || absMotherPDG == 321) { // K0s / K+ / K-
1387 if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.768049;
1388 else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.732933;
1389 else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.650298;
1390 else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.571332;
1391 else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.518734;
1392 else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.492543;
1393 else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.482704;
1394 else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.488056;
1395 else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.488861;
1396 else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.492862;
1397 else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.504332;
1398 else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.501858;
1399 else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.512970;
1400 else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.524131;
1401 else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.539130;
1402 else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.554101;
1403 else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.560348;
1404 else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.568869;
1405 else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.583310;
1406 else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.604818;
1407 else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.632630;
1408 else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.710070;
1409 else if(6.00 <= motherGenPt && motherGenPt < 8.00) fac = 0.736365;
1410 else if(8.00 <= motherGenPt && motherGenPt < 10.00) fac = 0.835865;
1411 }
1412
1413 if (absMotherPDG == 3122) { // Lambda
1414 if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.645162;
1415 else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.627431;
1416 else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.457136;
1417 else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.384369;
1418 else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.330597;
1419 else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.309571;
1420 else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.293620;
1421 else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.283709;
1422 else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.282047;
1423 else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.277261;
1424 else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.275772;
1425 else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.280726;
1426 else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.288540;
1427 else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.288315;
1428 else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.296619;
1429 else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.302993;
1430 else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.338121;
1431 else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.349800;
1432 else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.356802;
1433 else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.391202;
1434 else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.422573;
1435 else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.573815;
1436 else if(6.00 <= motherGenPt && motherGenPt < 8.00) fac = 0.786984;
1437 else if(8.00 <= motherGenPt && motherGenPt < 10.00) fac = 1.020021;
1438 }
1439
1440 if (absMotherPDG == 3312 || absMotherPDG == 3322) { // xi
1441 if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.666620;
1442 else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.575908;
1443 else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.433198;
1444 else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.340901;
1445 else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.290896;
1446 else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.236074;
1447 else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.218681;
1448 else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.207763;
1449 else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.222848;
1450 else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.208806;
1451 else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.197275;
1452 else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.183645;
1453 else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.188788;
1454 else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.188282;
1455 else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.207442;
1456 else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.240388;
1457 else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.241916;
1458 else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.208276;
1459 else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.234550;
1460 else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.251689;
1461 else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.310204;
1462 else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.343492;
1463 }
1464
1465 const Double_t weight = 1. / fac;
1466
1467 return weight;
1468}
1469
1470
1471//_____________________________________________________________________________
1472Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(AliMCEvent* mcEvent, AliMCParticle* daughter)
1473{
1474 // Strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
1475 // -> Based on function in PWGJE/AliAnalysisTaskFragmentationFunction
1476
1477 if (!mcEvent)
1478 return 1.;
1479
1480 AliMCParticle* currentMother = daughter;
1481 AliMCParticle* currentDaughter = daughter;
1482
1483
1484 // find first primary mother K0s, Lambda or Xi
1485 while(1) {
1486 Int_t daughterPDG = currentDaughter->PdgCode();
1487
1488 Int_t motherLabel = currentDaughter->GetMother();
1489 if(motherLabel >= mcEvent->GetNumberOfTracks()){ // protection
1490 currentMother = currentDaughter;
1491 break;
1492 }
1493
1494 currentMother = (AliMCParticle*)mcEvent->GetTrack(motherLabel);
1495
1496 if (!currentMother) {
1497 currentMother = currentDaughter;
1498 break;
1499 }
1500
1501 Int_t motherPDG = currentMother->PdgCode();
1502
1503 // phys. primary found ?
1504 if (mcEvent->IsPhysicalPrimary(motherLabel))
1505 break;
1506
1507 if (TMath::Abs(daughterPDG) == 321) {
1508 // K+/K- e.g. from phi (ref data not feeddown corrected)
1509 currentMother = currentDaughter;
1510 break;
1511 }
1512 if (TMath::Abs(motherPDG) == 310) {
1513 // K0s e.g. from phi (ref data not feeddown corrected)
1514 break;
1515 }
1516 if (TMath::Abs(motherPDG) == 3212 && TMath::Abs(daughterPDG) == 3122) {
1517 // Mother Sigma0, daughter Lambda (this case not included in feeddown corr.)
1518 currentMother = currentDaughter;
1519 break;
1520 }
1521
1522 currentDaughter = currentMother;
1523 }
1524
1525
1526 Int_t motherPDG = currentMother->PdgCode();
1527 Double_t motherGenPt = currentMother->Pt();
1528
1529 return GetMCStrangenessFactorCMS(motherPDG, motherGenPt);
1530}
1531
1532
1533// _________________________________________________________________________________
1534Bool_t AliAnalysisTaskPID::IsSecondaryWithStrangeMotherMC(AliMCEvent* mcEvent, Int_t partLabel)
1535{
1536 // Check whether particle is a secondary with strange mother, i.e. returns kTRUE if a strange mother is found
1537 // and the particle is NOT a physical primary. In all other cases kFALSE is returned
1538
1539 if (!mcEvent || partLabel < 0)
1540 return kFALSE;
1541
1542 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(partLabel);
1543
1544 if (!part)
1545 return kFALSE;
1546
1547 if (mcEvent->IsPhysicalPrimary(partLabel))
1548 return kFALSE;
1549
1550 Int_t iMother = part->GetMother();
1551 if (iMother < 0)
1552 return kFALSE;
1553
1554
1555 AliMCParticle* partM = (AliMCParticle*)mcEvent->GetTrack(iMother);
1556 if (!partM)
1557 return kFALSE;
1558
1559 Int_t codeM = TMath::Abs(partM->PdgCode());
1560 Int_t mfl = Int_t(codeM / TMath::Power(10, Int_t(TMath::Log10(codeM))));
1561 if (mfl == 3 && codeM != 3) // codeM = 3 is for s quark
1562 return kTRUE;
1563
1564 return kFALSE;
1565}
1566
1567
1568//_____________________________________________________________________________
1569Bool_t AliAnalysisTaskPID::SetParticleFractionHisto(const TH3D* hist, const Int_t species, const Bool_t sysError)
1570{
1571 // Store a clone of hist (containing the particle fractions of the corresponding species with statistical error (sysError = kFALSE)
1572 // or systematic error (sysError = kTRUE), respectively), internally
1573
1574 if (species < AliPID::kElectron || species > AliPID::kProton) {
1575 AliError(Form("Only fractions for species index %d to %d can be set, but not for the requested one: %d", 0,
1576 AliPID::kProton, species));
1577 return kFALSE;
1578 }
1579
1580 if (sysError) {
1581 delete fFractionSysErrorHists[species];
1582
1583 fFractionSysErrorHists[species] = new TH3D(*hist);
1584 }
1585 else {
1586 delete fFractionHists[species];
1587
1588 fFractionHists[species] = new TH3D(*hist);
1589 }
1590
1591 return kTRUE;
1592}
1593
1594
1595//_____________________________________________________________________________
1596Bool_t AliAnalysisTaskPID::SetParticleFractionHistosFromFile(const TString filePathName, const Bool_t sysError)
1597{
1598 // Loads particle fractions for all species from the desired file and returns kTRUE on success.
1599 // The maps are assumed to be of Type TH3D, to sit in the main directory and to have names
1600 // Form("hFraction_%e", AliPID::ParticleName(i)) for sysError = kFALSE and
1601 // Form("hFractionSysError_%e", AliPID::ParticleName(i)) for sysError = kTRUE.
1602
1603 TFile* f = TFile::Open(filePathName.Data());
1604 if (!f) {
1605 std::cout << "Failed to open file with particle fractions \"" << filePathName.Data() << "\"!" << std::endl;
1606 return kFALSE;
1607 }
1608
1609 TH3D* hist = 0x0;
1610 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1611 TString histName = Form("hFraction%s_%s", sysError ? "SysError" : "", AliPID::ParticleName(i));
1612 hist = dynamic_cast<TH3D*>(f->Get(histName.Data()));
1613 if (!hist) {
1614 std::cout << "Failed to load particle fractions for " << histName.Data() << "!";
1615 std::cout << std::endl << "Cleaning up particle fraction histos!" << std::endl;
1616 CleanupParticleFractionHistos();
1617 return kFALSE;
1618 }
1619
1620 if (!SetParticleFractionHisto(hist, i, sysError)) {
1621 std::cout << "Failed to load particle fractions for " << histName.Data() << "!";
1622 std::cout << std::endl << "Cleaning up particle fraction histos!" << std::endl;
1623 CleanupParticleFractionHistos();
1624 return kFALSE;
1625 }
1626 }
1627
1628 delete hist;
1629
1630 return kTRUE;
1631
1632}
1633
1634
1635//_____________________________________________________________________________
1636Int_t AliAnalysisTaskPID::GetRandomParticleTypeAccordingToParticleFractions(const Double_t trackPt, const Double_t jetPt,
1637 const Double_t centralityPercentile,
1638 const Bool_t smearByError,
1639 const Bool_t takeIntoAccountSysError) const
1640{
1641 // Uses the stored histograms with the particle fractions to generate a random particle type according to these fractions.
1642 // In case of problems (e.g. histo missing), AliPID::kUnknown is returned.
1643 // If smearByError is kTRUE, the used fractions will be random numbers distributed with a gauss with mean
1644 // being the corresponding particle fraction and sigma it's error.
1645 // Note that in this case only the fraction of a random species is varied in this way. The other fractions
1646 // will be re-normalised according their statistical errors.
1647 // The same holds for the systematic error of species "takeIntoAccountSpeciesSysError", but the random number will be
1648 // uniformly distributed within [mean - sys, mean + sys] and the re-normalisation will be weighted with the systematic errors.
1649 // Note that the fractions will be calculated first with only the systematic error taken into account (if desired), including
1650 // re-normalisation. Then, the resulting fractions will be used to calculate the final fractions - either with statistical error
1651 // or without. The species, for which the error will be used for smearing, is the same for sys and stat error.
1652
1653 Double_t prob[AliPID::kSPECIES];
1654 Int_t randomSpecies = (smearByError || takeIntoAccountSysError) ? (Int_t)(fRandom->Rndm() * AliPID::kSPECIES) : -1;
1655 Bool_t success = GetParticleFractions(trackPt, jetPt, centralityPercentile, prob, randomSpecies, randomSpecies);
1656
1657 if (!success)
1658 return AliPID::kUnknown;
1659
1660 Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
1661
1662 if (rnd <= prob[AliPID::kPion])
1663 return AliPID::kPion;
1664 else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon])
1665 return AliPID::kKaon;
1666 else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon] + prob[AliPID::kProton])
1667 return AliPID::kProton;
1668 else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon] + prob[AliPID::kProton] + prob[AliPID::kElectron])
1669 return AliPID::kElectron;
1670
1671 return AliPID::kMuon; //else it must be a muon (only species left)
1672}
1673
1674
1675//_____________________________________________________________________________
1676AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::GenerateDetectorResponse(const AliAnalysisTaskPID::ErrorCode errCode,
1677 const Double_t mean, const Double_t sigma,
1678 Double_t* responses, const Int_t nResponses,
1679 const Bool_t usePureGaus)
1680{
1681 // Generate detector response. If a previous generation was not successful or there is something wrong with this signal generation,
1682 // the function will return kFALSE
1683 if (!responses)
1684 return kError;
1685
1686 // Reset response array
1687 for (Int_t i = 0; i < nResponses; i++)
1688 responses[i] = -999;
1689
1690 if (errCode == kError)
1691 return kError;
1692
1693 ErrorCode ownErrCode = kNoErrors;
1694
1695 if (fUseConvolutedGaus && !usePureGaus) {
1696 // In case of convoluted gauss, calculate the probability density only once to save a lot of time!
1697
1698 TH1* hProbDensity = 0x0;
1699 ownErrCode = SetParamsForConvolutedGaus(mean, sigma);
1700 if (ownErrCode == kError)
1701 return kError;
1702
1703 hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
1704
1705 for (Int_t i = 0; i < nResponses; i++) {
1706 responses[i] = hProbDensity->GetRandom();
1707 //responses[i] fConvolutedGausDeltaPrime->GetRandom(); // MUCH slower than using the binned version via the histogram
1708 }
1709 }
1710 else {
1711 for (Int_t i = 0; i < nResponses; i++) {
1712 responses[i] = fRandom->Gaus(mean, sigma);
1713 }
1714 }
1715
1716 // If forwarded error code was a warning (error case has been handled before), return a warning
1717 if (errCode == kWarning)
1718 return kWarning;
1719
1720 return ownErrCode; // Forward success/warning
1721}
1722
1723
1724//_____________________________________________________________________________
1725void AliAnalysisTaskPID::PrintSettings(Bool_t printSystematicsSettings) const
1726{
1727 // Print current settings.
1728
1729 printf("\n\nSettings for task %s:\n", GetName());
1730 printf("Is pPb/Pbp: %d -> %s\n", GetIsPbpOrpPb(), GetIsPbpOrpPb() ? "Adapting vertex cuts" : "Using standard vertex cuts");
1731 printf("Track cuts: %s\n", fTrackFilter ? fTrackFilter->GetTitle() : "-");
1732 printf("Eta cut: %.2f <= |eta| <= %.2f\n", GetEtaAbsCutLow(), GetEtaAbsCutUp());
1733 printf("Phi' cut: %d\n", GetUsePhiCut());
1734
1735 printf("\n");
1736
1737 printf("Centrality estimator: %s\n", GetCentralityEstimator().Data());
1738
1739 printf("\n");
1740
1741 printf("Use MC-ID for signal generation: %d\n", GetUseMCidForGeneration());
1742 printf("Use ITS: %d\n", GetUseITS());
1743 printf("Use TOF: %d\n", GetUseTOF());
1744 printf("Use priors: %d\n", GetUsePriors());
1745 printf("Use TPC default priors: %d\n", GetUseTPCDefaultPriors());
1746 printf("Use convoluted Gauss: %d\n", GetUseConvolutedGaus());
1747 printf("Accuracy of non-Gaussian tail: %e\n", GetAccuracyNonGaussianTail());
1748 printf("Take into account muons: %d\n", GetTakeIntoAccountMuons());
1749 printf("\nParams for transition from gauss to asymmetric shape:\n");
1750 printf("[0]: %e\n", GetConvolutedGaussTransitionPar(0));
1751 printf("[1]: %e\n", GetConvolutedGaussTransitionPar(1));
1752 printf("[2]: %e\n", GetConvolutedGaussTransitionPar(2));
1753
1754 printf("\n");
1755
1756 printf("Input from other task: %d\n", GetInputFromOtherTask());
1757 printf("Store additional jet information: %d\n", GetStoreAdditionalJetInformation());
1758 printf("Store centrality percentile: %d", GetStoreCentralityPercentile());
1759
1760 if (printSystematicsSettings)
1761 PrintSystematicsSettings();
1762 else
1763 printf("\n\n\n");
1764}
1765
1766
1767//_____________________________________________________________________________
1768void AliAnalysisTaskPID::PrintSystematicsSettings() const
1769{
1770 // Print current settings for systematic studies.
1771
1772 printf("\n\nSettings for systematics for task %s:\n", GetName());
1773 printf("Splines:\t%f\n", GetSystematicScalingSplines());
1774 printf("EtaCorrMomThr:\t%f\n", GetSystematicScalingEtaCorrectionMomentumThr());
1775 printf("EtaCorrLowP:\t%f\n", GetSystematicScalingEtaCorrectionLowMomenta());
1776 printf("EtaCorrHighP:\t%f\n", GetSystematicScalingEtaCorrectionHighMomenta());
1777 printf("SigmaPara:\t%f\n", GetSystematicScalingEtaSigmaPara());
1778 printf("MultCorr:\t%f\n", GetSystematicScalingMultCorrection());
1779
1780 printf("\n\n");
1781}
1782
1783
1784//_____________________________________________________________________________
1785Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePDGcode, Double_t centralityPercentile,
1786 Double_t jetPt)
1787{
1788 // Process the track (generate expected response, fill histos, etc.).
1789 // particlePDGcode == 0 means data. Otherwise, the corresponding MC ID will be assumed.
1790
1791 //Printf("Debug: Task %s is starting to process track: dEdx %f, pTPC %f, eta %f, ncl %d\n", GetName(), track->GetTPCsignal(), track->GetTPCmomentum(),
1792 // track->Eta(), track->GetTPCsignalN());
1793
1794 const Bool_t isMC = (particlePDGcode == 0) ? kFALSE : kTRUE;
1795
1796 Int_t binMC = -1;
1797
1798 if (isMC) {
1799 if (TMath::Abs(particlePDGcode) == 211) {//Pion
1800 binMC = 3;
1801 }
1802 else if (TMath::Abs(particlePDGcode) == 321) {//Kaon
1803 binMC = 1;
1804 }
1805 else if (TMath::Abs(particlePDGcode) == 2212) {//Proton
1806 binMC = 4;
1807 }
1808 else if (TMath::Abs(particlePDGcode) == 11) {//Electron
1809 binMC = 0;
1810 }
1811 else if (TMath::Abs(particlePDGcode) == 13) {//Muon
1812 binMC = 2;
1813 }
1814 else // In MC-ID case, set to underflow bin such that the response from this track is only used for unidentified signal generation
1815 // or signal generation with PID response and the track is still there (as in data) - e.g. relevant w.r.t. deuterons.
1816 // This is important to be as much as possible consistent with data. And the tracks can still be removed by disabling the
1817 // underflow bin for the projections
1818 binMC = -1;
1819 }
1820
1821 // Momenta
1822 //Double_t p = track->GetP();
1823 //Double_t pTPC = track->GetTPCmomentum();
1824 Double_t pT = track->Pt();
1825
1826 Double_t z = -1, xi = -1;
1827 GetJetTrackObservables(pT, jetPt, z, xi);
1828
1829
1830 Double_t trackCharge = track->Charge();
1831
1832 // TPC signal
1833 Double_t dEdxTPC = fPIDResponse->IsTunedOnData() ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
1834
1835 if (dEdxTPC <= 0) {
1836 Printf("Skipping track with strange dEdx value: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), track->GetTPCmomentum(),
1837 track->Eta(), track->GetTPCsignalN());
1838 return kFALSE;
1839 }
1840
1841
1842
1843
1844 Double_t dEdxEl, dEdxKa, dEdxPi, dEdxMu, dEdxPr;
1845 Double_t sigmaEl, sigmaKa, sigmaPi, sigmaMu, sigmaPr;
1846
1847 if (fDoAnySystematicStudiesOnTheExpectedSignal) {
1848 // Get the uncorrected signal first and the corresponding correction factors.
1849 // Then modify the correction factors and properly recalculate the corrected dEdx
1850
1851 // Get pure spline values for dEdx_expected, without any correction
1852 dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
1853 dEdxKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
1854 dEdxPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
1855 dEdxMu = !fTakeIntoAccountMuons ? -1 :
1856 fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
1857 dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
1858
1859 // Scale splines, if desired
1860 if (TMath::Abs(fSystematicScalingSplines - 1.0) > fgkEpsilon) {
1861 dEdxEl *= fSystematicScalingSplines;
1862 dEdxKa *= fSystematicScalingSplines;
1863 dEdxPi *= fSystematicScalingSplines;
1864 dEdxMu *= fTakeIntoAccountMuons ? fSystematicScalingSplines : 1.;
1865 dEdxPr *= fSystematicScalingSplines;
1866 }
1867
1868 // Get the eta correction factors for the (modified) expected dEdx
1869 Double_t etaCorrEl = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxEl) : 1.;
1870 Double_t etaCorrKa = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxKa) : 1.;
1871 Double_t etaCorrPi = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPi) : 1.;
1872 Double_t etaCorrMu = fTakeIntoAccountMuons && !fPIDResponse->UseTPCEtaCorrection() ?
1873 fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxMu) : 1.;
1874 Double_t etaCorrPr = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPr) : 1.;
1875
1876 // Scale eta correction factors, if desired (and eta correction maps are to be used, otherwise it is not possible!)
1877 if (fPIDResponse->UseTPCEtaCorrection() &&
1878 (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon ||
1879 TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon)) {
1880 // Since we do not want to scale the splines with this, but only the eta variation, only scale the deviation of the correction factor!
1881 // 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!
1882
1883
1884 // Due to additional azimuthal effects, there is an additional eta dependence for low momenta which is not corrected successfully so far.
1885 // One can assign a different (higher) systematic scale factor for this low-p region and a threshold which separates low- and high-p.
1886 // 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
1887 Double_t usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionHighMomenta;
1888
1889 if (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - fSystematicScalingEtaCorrectionLowMomenta) > fgkEpsilon) {
1890 const Double_t pTPC = track->GetTPCmomentum();
1891 const Double_t fractionHighMomentumScaleFactor = 0.5 * (1. + TMath::Erf((pTPC - fSystematicScalingEtaCorrectionMomentumThr) / 0.1));
1892 usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionLowMomenta * (1 - fractionHighMomentumScaleFactor)
1893 + fSystematicScalingEtaCorrectionHighMomenta * fractionHighMomentumScaleFactor;
1894 }
1895
1896 etaCorrEl = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrEl - 1.0);
1897 etaCorrKa = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrKa - 1.0);
1898 etaCorrPi = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPi - 1.0);
1899 etaCorrMu = fTakeIntoAccountMuons ? (1.0 + usedSystematicScalingEtaCorrection * (etaCorrMu - 1.0)) : 1.0;
1900 etaCorrPr = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPr - 1.0);
1901 }
1902
1903 // Get the multiplicity correction factors for the (modified) expected dEdx
1904 const Int_t currEvtMultiplicity = fPIDResponse->GetTPCResponse().GetCurrentEventMultiplicity();
1905
1906 Double_t multiplicityCorrEl = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
1907 dEdxEl * etaCorrEl, currEvtMultiplicity) : 1.;
1908 Double_t multiplicityCorrKa = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
1909 dEdxKa * etaCorrKa, currEvtMultiplicity) : 1.;
1910 Double_t multiplicityCorrPi = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
1911 dEdxPi * etaCorrPi, currEvtMultiplicity) : 1.;
1912 Double_t multiplicityCorrMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCMultiplicityCorrection() ?
1913 fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track, dEdxMu * etaCorrMu, currEvtMultiplicity) : 1.;
1914 Double_t multiplicityCorrPr = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
1915 dEdxPr * etaCorrPr, currEvtMultiplicity) : 1.;
1916
1917 Double_t multiplicityCorrSigmaEl = fPIDResponse->UseTPCMultiplicityCorrection() ?
1918 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxEl * etaCorrEl, currEvtMultiplicity) : 1.;
1919 Double_t multiplicityCorrSigmaKa = fPIDResponse->UseTPCMultiplicityCorrection() ?
1920 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxKa * etaCorrKa, currEvtMultiplicity) : 1.;
1921 Double_t multiplicityCorrSigmaPi = fPIDResponse->UseTPCMultiplicityCorrection() ?
1922 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxPi * etaCorrPi, currEvtMultiplicity) : 1.;
1923 Double_t multiplicityCorrSigmaMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCMultiplicityCorrection() ?
1924 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxMu * etaCorrMu, currEvtMultiplicity) : 1.;
1925 Double_t multiplicityCorrSigmaPr = fPIDResponse->UseTPCMultiplicityCorrection() ?
1926 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxPr * etaCorrPr, currEvtMultiplicity) : 1.;
1927
1928 // Scale multiplicity correction factors, if desired (and multiplicity correction functions are to be used, otherwise it is not possible!)
1929 if (fPIDResponse->UseTPCMultiplicityCorrection() && TMath::Abs(fSystematicScalingMultCorrection - 1.0) > fgkEpsilon) {
1930 // Since we do not want to scale the splines with this, but only the multiplicity variation, only scale the deviation of the correction factor!
1931 // 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!
1932
1933 multiplicityCorrEl = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrEl - 1.0);
1934 multiplicityCorrKa = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrKa - 1.0);
1935 multiplicityCorrPi = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrPi - 1.0);
1936 multiplicityCorrMu = fTakeIntoAccountMuons ? (1.0 + fSystematicScalingMultCorrection * (multiplicityCorrMu - 1.0)) : 1.0;
1937 multiplicityCorrPr = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrPr - 1.0);
1938 }
1939
1940 // eta correction must be enabled in order to use the new sigma parametrisation maps. Since this is the absolute sigma
1941 // for a track calculated with the unscaled paramaters, we have to devide out dEdxExpectedEtaCorrected and then need
1942 // to scale with the multiplicitySigmaCorrFactor * fSystematicScalingEtaSigmaPara. In the end, one has to scale with the
1943 // (modified) dEdx to get the absolute sigma
1944 // This means there is no extra parameter for the multiplicitySigmaCorrFactor, but only for the sigma map itself.
1945 // This is valid, since it appears only as a product. One has to assume a larger systematic shift in case of additional
1946 // multiplicity dependence....
1947 Double_t sigmaRelEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1948 / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1949 * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaEl;
1950
1951 Double_t sigmaRelKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1952 / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1953 * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaKa;
1954
1955 Double_t sigmaRelPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1956 / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1957 * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaPi;
1958
1959 Double_t sigmaRelMu = fTakeIntoAccountMuons ?
1960 fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1961 / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1962 * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaMu
1963 : 999.;
1964
1965 Double_t sigmaRelPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1966 / fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
1967 * fSystematicScalingEtaSigmaPara * multiplicityCorrSigmaPr;
1968
1969 // Now scale the (possibly modified) spline values with the (possibly modified) correction factors
1970 dEdxEl *= etaCorrEl * multiplicityCorrEl;
1971 dEdxKa *= etaCorrKa * multiplicityCorrKa;
1972 dEdxPi *= etaCorrPi * multiplicityCorrPi;
1973 dEdxMu *= etaCorrMu * multiplicityCorrMu;
1974 dEdxPr *= etaCorrPr * multiplicityCorrPr;
1975
1976 // Finally, get the absolute sigma
1977 sigmaEl = sigmaRelEl * dEdxEl;
1978 sigmaKa = sigmaRelKa * dEdxKa;
1979 sigmaPi = sigmaRelPi * dEdxPi;
1980 sigmaMu = sigmaRelMu * dEdxMu;
1981 sigmaPr = sigmaRelPr * dEdxPr;
1982
1983 }
1984 else {
1985 // No systematic studies on expected signal - just take it as it comve from the TPCPIDResponse
1986 dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
1987 fPIDResponse->UseTPCEtaCorrection(),
1988 fPIDResponse->UseTPCMultiplicityCorrection());
1989 dEdxKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
1990 fPIDResponse->UseTPCEtaCorrection(),
1991 fPIDResponse->UseTPCMultiplicityCorrection());
1992 dEdxPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
1993 fPIDResponse->UseTPCEtaCorrection(),
1994 fPIDResponse->UseTPCMultiplicityCorrection());
1995 dEdxMu = !fTakeIntoAccountMuons ? -1 :
1996 fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
1997 fPIDResponse->UseTPCEtaCorrection(),
1998 fPIDResponse->UseTPCMultiplicityCorrection());
1999 dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2000 fPIDResponse->UseTPCEtaCorrection(),
2001 fPIDResponse->UseTPCMultiplicityCorrection());
2002
2003 sigmaEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2004 fPIDResponse->UseTPCEtaCorrection(),
2005 fPIDResponse->UseTPCMultiplicityCorrection());
2006 sigmaKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2007 fPIDResponse->UseTPCEtaCorrection(),
2008 fPIDResponse->UseTPCMultiplicityCorrection());
2009 sigmaPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2010 fPIDResponse->UseTPCEtaCorrection(),
2011 fPIDResponse->UseTPCMultiplicityCorrection());
2012 sigmaMu = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2013 fPIDResponse->UseTPCEtaCorrection(),
2014 fPIDResponse->UseTPCMultiplicityCorrection());
2015 sigmaPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2016 fPIDResponse->UseTPCEtaCorrection(),
2017 fPIDResponse->UseTPCMultiplicityCorrection());
2018 }
2019 /*OLD with deltaSpecies
2020 Double_t deltaElectron = dEdxTPC - dEdxEl;
2021 Double_t deltaKaon = dEdxTPC - dEdxKa;
2022 Double_t deltaPion = dEdxTPC - dEdxPi;
2023 Double_t deltaProton = dEdxTPC - dEdxPr;
2024 */
2025
2026 Double_t deltaPrimeElectron = (dEdxEl > 0) ? dEdxTPC / dEdxEl : -1;
2027 if (dEdxEl <= 0) {
2028 Printf("Error: Expected TPC signal <= 0 for electron hypothesis");
2029 return kFALSE;
2030 }
2031
2032 Double_t deltaPrimeKaon = (dEdxKa > 0) ? dEdxTPC / dEdxKa : -1;
2033 if (dEdxKa <= 0) {
2034 Printf("Error: Expected TPC signal <= 0 for kaon hypothesis");
2035 return kFALSE;
2036 }
2037
2038 Double_t deltaPrimePion = (dEdxPi > 0) ? dEdxTPC / dEdxPi : -1;
2039 if (dEdxPi <= 0) {
2040 Printf("Error: Expected TPC signal <= 0 for pion hypothesis");
2041 return kFALSE;
2042 }
2043
2044 Double_t deltaPrimeProton = (dEdxPr > 0) ? dEdxTPC / dEdxPr : -1;
2045 if (dEdxPr <= 0) {
2046 Printf("Error: Expected TPC signal <= 0 for proton hypothesis");
2047 return kFALSE;
2048 }
2049
2050 /*TODO for TOF
2051 // TOF signal
2052 Double_t times[AliPID::kSPECIES];
2053 track->GetIntegratedTimes(times);
2054 AliTOFPIDResponse tofPIDResponse = fPIDResponse->GetTOFResponse();
2055 Float_t electronDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kElectron);
2056 Float_t pionDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kPion);
2057 Float_t kaonDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kKaon);
2058 Float_t protonDeltaTOF = GetDeltaTOF(track, &tofPIDResponse, times, AliPID::kProton);
2059 */
2060
2061
2062 // Use probabilities to weigh the response generation for the different species.
2063 // Also determine the most probable particle type.
2064 Double_t prob[AliPID::kSPECIESC];
2065 for (Int_t i = 0; i < AliPID::kSPECIESC; i++)
2066 prob[i] = 0;
2067
2068 fPIDcombined->ComputeProbabilities(track, fPIDResponse, prob);
2069
2070 // Bug: One can set the number of species for PIDcombined, but PIDcombined will call PIDresponse, which writes without testing
2071 // the probs for kSPECIESC (including light nuclei) into the array.
2072 // In this case, when only kSPECIES are considered, the probabilities have to be rescaled!
2073 for (Int_t i = AliPID::kSPECIES; i < AliPID::kSPECIESC; i++)
2074 prob[i] = 0;
2075
2076 // If muons are not to be taken into account, just set their probability to zero and normalise the remaining probabilities
2077 if (!fTakeIntoAccountMuons)
2078 prob[AliPID::kMuon] = 0;
2079
2080 Double_t probSum = 0.;
2081 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2082 probSum += prob[i];
2083
2084 if (probSum > 0) {
2085 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2086 prob[i] /= probSum;
2087 }
2088
2089 if (!isMC) {
2090 // If there is no MC information, take the most probable species for the ID
2091 Float_t max = 0.;
2092 Int_t maxIndex = -1;
2093 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2094 if (prob[i] > max) {
2095 max = prob[i];
2096 maxIndex = i;
2097 }
2098 }
2099
2100 // Translate from AliPID numbering to numbering of this class
2101 if (max > 0) {
2102 if (maxIndex == AliPID::kElectron)
2103 binMC = 0;
2104 else if (maxIndex == AliPID::kKaon)
2105 binMC = 1;
2106 else if (maxIndex == AliPID::kMuon)
2107 binMC = 2;
2108 else if (maxIndex == AliPID::kPion)
2109 binMC = 3;
2110 else if (maxIndex == AliPID::kProton)
2111 binMC = 4;
2112 else
2113 binMC = -1;
2114 }
2115 else {
2116 // Only take track into account for expectation values, if valid pid response is available.. Otherwise: Set to underflow bin.
2117 binMC = -1;
2118 }
2119 }
2120
2121 /*
2122 //For testing: Swap p<->pT to analyse pure p-dependence => Needs to be removed later
2123 Double_t temp = pT;
2124 pT = pTPC;
2125 pTPC = temp;
2126 */
2127
2128
2129 Double_t entry[fStoreAdditionalJetInformation ? kDataNumAxes : kDataNumAxes - fgkNumJetAxes];
2130 entry[kDataMCID] = binMC;
2131 entry[kDataSelectSpecies] = 0;
2132 entry[kDataPt] = pT;
2133 entry[kDataDeltaPrimeSpecies] = deltaPrimeElectron;
2134 entry[kDataCentrality] = centralityPercentile;
2135
2136 if (fStoreAdditionalJetInformation) {
2137 entry[kDataJetPt] = jetPt;
2138 entry[kDataZ] = z;
2139 entry[kDataXi] = xi;
2140 }
2141
2142 entry[GetIndexOfChargeAxisData()] = trackCharge;
2143
2144 /*OLD with TOF, p_TPC_Inner and p_vertex and deltaSpecies
2145 // MC PID, SelectSpecies, P(TPC_inner), pT, p(Vertex), deltaSpecies, deltaPrimeSpecies, deltaTOFspecies
2146 Double_t entry[8] = { binMC, 0, pTPC, pT, p, deltaElectron, deltaPrimeElectron, electronDeltaTOF };
2147 */
2148
2149 fhPIDdataAll->Fill(entry);
2150
2151 entry[kDataSelectSpecies] = 1;
2152 //OLD with deltaSpecies entry[5] = deltaKaon;
2153 entry[kDataDeltaPrimeSpecies] = deltaPrimeKaon;
2154 //TODO for TOF entry[7] = kaonDeltaTOF;
2155 fhPIDdataAll->Fill(entry);
2156
2157 entry[kDataSelectSpecies] = 2;
2158 //OLD with deltaSpecies entry[5] = deltaPion;
2159 entry[kDataDeltaPrimeSpecies] = deltaPrimePion;
2160 //TODO for TOF entry[7] = pionDeltaTOF;
2161 fhPIDdataAll->Fill(entry);
2162
2163 entry[kDataSelectSpecies] = 3;
2164 //OLD with deltaSpecies entry[5] = deltaProton;
2165 entry[kDataDeltaPrimeSpecies] = deltaPrimeProton;
2166 //TODO for TOF entry[7] = protonDeltaTOF;
2167 fhPIDdataAll->Fill(entry);
2168
2169
2170 // Construct the expected shape for the transition p -> pT
2171
2172 Double_t genEntry[fStoreAdditionalJetInformation ? kGenNumAxes : kGenNumAxes - fgkNumJetAxes];
2173 genEntry[kGenMCID] = binMC;
2174 genEntry[kGenSelectSpecies] = 0;
2175 genEntry[kGenPt] = pT;
2176 genEntry[kGenDeltaPrimeSpecies] = -999;
2177 genEntry[kGenCentrality] = centralityPercentile;
2178
2179 if (fStoreAdditionalJetInformation) {
2180 genEntry[kGenJetPt] = jetPt;
2181 genEntry[kGenZ] = z;
2182 genEntry[kGenXi] = xi;
2183 }
2184
2185 genEntry[GetIndexOfChargeAxisGen()] = trackCharge;
2186
2187 //OLD with deltaSpecies Double_t genEntry[5] = { binMC, 0, pT, -999, -999 }; // MC PID, SelectSpecies, pT, deltaSpecies, deltaPrimeSpecies
2188
2189 // Generate numGenEntries "responses" with fluctuations around the expected values.
2190 // The higher the (transverse) momentum, the more "responses" will be generated in order not to run out of statistics too fast.
2191 Int_t numGenEntries = 10;
2192 if (pT >= 5)
2193 numGenEntries *= 5;
2194 else if (pT >= 2)
2195 numGenEntries *= 2;
2196
2197 // Jets have even worse statistics, therefore, scale numGenEntries further
2198 if (jetPt >= 40)
2199 numGenEntries *= 20;
2200 else if (jetPt >= 20)
2201 numGenEntries *= 10;
2202 else if (jetPt >= 10)
2203 numGenEntries *= 2;
2204
2205
2206 // Do not generate more entries than available in memory!
2207 if (numGenEntries > fgkMaxNumGenEntries)// fgkMaxNumGenEntries = 1000
2208 numGenEntries = fgkMaxNumGenEntries;
2209
2210 ErrorCode errCode = kNoErrors;
2211
2212 // Electrons
2213 errCode = GenerateDetectorResponse(errCode, 1., sigmaEl / dEdxEl, fGenRespElDeltaPrimeEl, numGenEntries);
2214 errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxKa, sigmaEl / dEdxKa, fGenRespElDeltaPrimeKa, numGenEntries);
2215 errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxPi, sigmaEl / dEdxPi, fGenRespElDeltaPrimePi, numGenEntries);
2216 errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxPr, sigmaEl / dEdxPr, fGenRespElDeltaPrimePr, numGenEntries);
2217
2218 // Kaons
2219 errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxEl, sigmaKa / dEdxEl, fGenRespKaDeltaPrimeEl, numGenEntries);
2220 errCode = GenerateDetectorResponse(errCode, 1., sigmaKa / dEdxKa, fGenRespKaDeltaPrimeKa, numGenEntries);
2221 errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxPi, sigmaKa / dEdxPi, fGenRespKaDeltaPrimePi, numGenEntries);
2222 errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxPr, sigmaKa / dEdxPr, fGenRespKaDeltaPrimePr, numGenEntries);
2223
2224 // Pions
2225 errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxEl, sigmaPi / dEdxEl, fGenRespPiDeltaPrimeEl, numGenEntries);
2226 errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxKa, sigmaPi / dEdxKa, fGenRespPiDeltaPrimeKa, numGenEntries);
2227 errCode = GenerateDetectorResponse(errCode, 1., sigmaPi / dEdxPi, fGenRespPiDeltaPrimePi, numGenEntries);
2228 errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxPr, sigmaPi / dEdxPr, fGenRespPiDeltaPrimePr, numGenEntries);
2229
2230 // Muons, if desired
2231 if (fTakeIntoAccountMuons) {
2232 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxEl, sigmaMu / dEdxEl, fGenRespMuDeltaPrimeEl, numGenEntries);
2233 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxKa, sigmaMu / dEdxKa, fGenRespMuDeltaPrimeKa, numGenEntries);
2234 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxPi, sigmaMu / dEdxPi, fGenRespMuDeltaPrimePi, numGenEntries);
2235 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxPr, sigmaMu / dEdxPr, fGenRespMuDeltaPrimePr, numGenEntries);
2236 }
2237
2238 // Protons
2239 errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxEl, sigmaPr / dEdxEl, fGenRespPrDeltaPrimeEl, numGenEntries);
2240 errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxKa, sigmaPr / dEdxKa, fGenRespPrDeltaPrimeKa, numGenEntries);
2241 errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxPi, sigmaPr / dEdxPi, fGenRespPrDeltaPrimePi, numGenEntries);
2242 errCode = GenerateDetectorResponse(errCode, 1., sigmaPr / dEdxPr, fGenRespPrDeltaPrimePr, numGenEntries);
2243
2244
2245 /*OLD with deltaSpecies
2246 // Delta
2247
2248 // Electrons
2249 errCode = GenerateDetectorResponse(errCode, 0., sigmaEl, fGenRespElDeltaEl, numGenEntries, usePureGausForDelta);
2250 errCode = GenerateDetectorResponse(errCode, dEdxEl - dEdxKa, sigmaEl, fGenRespElDeltaKa, numGenEntries, usePureGausForDelta);
2251 errCode = GenerateDetectorResponse(errCode, dEdxEl - dEdxPi, sigmaEl, fGenRespElDeltaPi, numGenEntries, usePureGausForDelta);
2252 errCode = GenerateDetectorResponse(errCode, dEdxEl - dEdxPr, sigmaEl, fGenRespElDeltaPr, numGenEntries, usePureGausForDelta);
2253
2254 // Kaons
2255 errCode = GenerateDetectorResponse(errCode, dEdxKa - dEdxEl, sigmaKa, fGenRespKaDeltaEl, numGenEntries, usePureGausForDelta);
2256 errCode = GenerateDetectorResponse(errCode, 0., sigmaKa, fGenRespKaDeltaKa, numGenEntries, usePureGausForDelta);
2257 errCode = GenerateDetectorResponse(errCode, dEdxKa - dEdxPi, sigmaKa, fGenRespKaDeltaPi, numGenEntries, usePureGausForDelta);
2258 errCode = GenerateDetectorResponse(errCode, dEdxKa - dEdxPr, sigmaKa, fGenRespKaDeltaPr, numGenEntries, usePureGausForDelta);
2259
2260 // Pions
2261 errCode = GenerateDetectorResponse(errCode, dEdxPi - dEdxEl, sigmaPi, fGenRespPiDeltaEl, numGenEntries, usePureGausForDelta);
2262 errCode = GenerateDetectorResponse(errCode, dEdxPi - dEdxKa, sigmaPi, fGenRespPiDeltaKa, numGenEntries, usePureGausForDelta);
2263 errCode = GenerateDetectorResponse(errCode, 0., sigmaPi, fGenRespPiDeltaPi, numGenEntries, usePureGausForDelta);
2264 errCode = GenerateDetectorResponse(errCode, dEdxPi - dEdxPr, sigmaPi, fGenRespPiDeltaPr, numGenEntries, usePureGausForDelta);
2265
2266 // Muons
2267 errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxEl, sigmaMu, fGenRespMuDeltaEl, numGenEntries, usePureGausForDelta);
2268 errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxKa, sigmaMu, fGenRespMuDeltaKa, numGenEntries, usePureGausForDelta);
2269 errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxPi, sigmaMu, fGenRespMuDeltaPi, numGenEntries, usePureGausForDelta);
2270 errCode = GenerateDetectorResponse(errCode, dEdxMu - dEdxPr, sigmaMu, fGenRespMuDeltaPr, numGenEntries, usePureGausForDelta);
2271
2272 // Protons
2273 errCode = GenerateDetectorResponse(errCode, dEdxPr - dEdxEl, sigmaPr, fGenRespPrDeltaEl, numGenEntries, usePureGausForDelta);
2274 errCode = GenerateDetectorResponse(errCode, dEdxPr - dEdxKa, sigmaPr, fGenRespPrDeltaKa, numGenEntries, usePureGausForDelta);
2275 errCode = GenerateDetectorResponse(errCode, dEdxPr - dEdxPi, sigmaPr, fGenRespPrDeltaPi, numGenEntries, usePureGausForDelta);
2276 errCode = GenerateDetectorResponse(errCode, 0., sigmaPr, fGenRespPrDeltaPr, numGenEntries, usePureGausForDelta);
2277 */
2278
2279 if (errCode != kNoErrors) {
2280 if (errCode == kWarning) {
2281 //Printf("Warning: Questionable detector response for track, most likely due to very low number of PID clusters! Debug output (dEdx_expected, sigma_expected):");
2282 }
2283 else
2284 Printf("Error: Failed to generate detector response for track - skipped! Debug output (dEdx_expected, sigma_expected):");
2285
2286 /*
2287 Printf("Pr: %e, %e", dEdxPr, sigmaPr);
2288 Printf("Pi: %e, %e", dEdxPi, sigmaPi);
2289 Printf("El: %e, %e", dEdxEl, sigmaEl);
2290 Printf("Mu: %e, %e", dEdxMu, sigmaMu);
2291 Printf("Ka: %e, %e", dEdxKa, sigmaKa);
2292 Printf("track: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), track->GetTPCmomentum(), track->Eta(),
2293 track->GetTPCsignalN());
2294 */
2295
2296 if (errCode != kWarning) {
2297 fhSkippedTracksForSignalGeneration->Fill(track->GetTPCmomentum(), track->GetTPCsignalN());
2298 return kFALSE;// Skip generated response in case of error
2299 }
2300 }
2301
2302 for (Int_t n = 0; n < numGenEntries; n++) {
2303 if (!isMC || !fUseMCidForGeneration) {
2304 // If no MC info is available or shall not be used, use weighting with priors to generate entries for the different species
2305 Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
2306
2307 // Consider generated response as originating from...
2308 if (rnd <= prob[AliPID::kElectron])
2309 genEntry[kGenMCID] = 0; // ... an electron
2310 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon])
2311 genEntry[kGenMCID] = 1; // ... a kaon
2312 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon])
2313 genEntry[kGenMCID] = 2; // ... a muon -> NOTE: prob[AliPID::kMuon] = 0 in case of fTakeIntoAccountMuons = kFALSE
2314 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon] + prob[AliPID::kPion])
2315 genEntry[kGenMCID] = 3; // ... a pion
2316 else
2317 genEntry[kGenMCID] = 4; // ... a proton
2318 }
2319
2320 // Electrons
2321 genEntry[kGenSelectSpecies] = 0;
2322 //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaEl[n];
2323 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeEl[n];
2324 fhGenEl->Fill(genEntry);
2325
2326 genEntry[kGenSelectSpecies] = 1;
2327 //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaKa[n];
2328 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeKa[n];
2329 fhGenEl->Fill(genEntry);
2330
2331 genEntry[kGenSelectSpecies] = 2;
2332 //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaPi[n];
2333 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePi[n];
2334 fhGenEl->Fill(genEntry);
2335
2336 genEntry[kGenSelectSpecies] = 3;
2337 //OLD with deltaSpecies genEntry[3] = fGenRespElDeltaPr[n];
2338 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePr[n];
2339 fhGenEl->Fill(genEntry);
2340
2341 // Kaons
2342 genEntry[kGenSelectSpecies] = 0;
2343 //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaEl[n];
2344 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeEl[n];
2345 fhGenKa->Fill(genEntry);
2346
2347 genEntry[kGenSelectSpecies] = 1;
2348 //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaKa[n];
2349 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeKa[n];
2350 fhGenKa->Fill(genEntry);
2351
2352 genEntry[kGenSelectSpecies] = 2;
2353 //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaPi[n];
2354 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePi[n];
2355 fhGenKa->Fill(genEntry);
2356
2357 genEntry[kGenSelectSpecies] = 3;
2358 //OLD with deltaSpecies genEntry[3] = fGenRespKaDeltaPr[n];
2359 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePr[n];
2360 fhGenKa->Fill(genEntry);
2361
2362 // Pions
2363 genEntry[kGenSelectSpecies] = 0;
2364 //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaEl[n];
2365 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeEl[n];
2366 fhGenPi->Fill(genEntry);
2367
2368 genEntry[kGenSelectSpecies] = 1;
2369 //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaKa[n];
2370 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeKa[n];
2371 fhGenPi->Fill(genEntry);
2372
2373 genEntry[kGenSelectSpecies] = 2;
2374 //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaPi[n];
2375 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePi[n];
2376 fhGenPi->Fill(genEntry);
2377
2378 genEntry[kGenSelectSpecies] = 3;
2379 //OLD with deltaSpecies genEntry[3] = fGenRespPiDeltaPr[n];
2380 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePr[n];
2381 fhGenPi->Fill(genEntry);
2382
2383 if (fTakeIntoAccountMuons) {
2384 // Muons, if desired
2385 genEntry[kGenSelectSpecies] = 0;
2386 //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaEl[n];
2387 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeEl[n];
2388 fhGenMu->Fill(genEntry);
2389
2390 genEntry[kGenSelectSpecies] = 1;
2391 //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaKa[n];
2392 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeKa[n];
2393 fhGenMu->Fill(genEntry);
2394
2395 genEntry[kGenSelectSpecies] = 2;
2396 //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaPi[n];
2397 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePi[n];
2398 fhGenMu->Fill(genEntry);
2399
2400 genEntry[kGenSelectSpecies] = 3;
2401 //OLD with deltaSpecies genEntry[3] = fGenRespMuDeltaPr[n];
2402 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePr[n];
2403 fhGenMu->Fill(genEntry);
2404 }
2405
2406 // Protons
2407 genEntry[kGenSelectSpecies] = 0;
2408 //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaEl[n];
2409 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeEl[n];
2410 fhGenPr->Fill(genEntry);
2411
2412 genEntry[kGenSelectSpecies] = 1;
2413 //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaKa[n];
2414 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeKa[n];
2415 fhGenPr->Fill(genEntry);
2416
2417 genEntry[kGenSelectSpecies] = 2;
2418 //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaPi[n];
2419 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePi[n];
2420 fhGenPr->Fill(genEntry);
2421
2422 genEntry[kGenSelectSpecies] = 3;
2423 //OLD with deltaSpecies genEntry[3] = fGenRespPrDeltaPr[n];
2424 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePr[n];
2425 fhGenPr->Fill(genEntry);
2426 }
2427
2428 return kTRUE;
2429}
2430
2431
2432//_____________________________________________________________________________
2433Bool_t AliAnalysisTaskPID::SetConvolutedGaussLambdaParameter(Double_t lambda)
2434{
2435 // Set the lambda parameter of the convolution to the desired value. Automatically
2436 // calculates the parameters for the transition (restricted) gauss -> convoluted gauss.
2437 fConvolutedGaussTransitionPars[2] = lambda;
2438
2439 // Save old parameters and settings of function to restore them later:
2440 Double_t* oldFuncParams = new Double_t[fkConvolutedGausNPar];
2441 for (Int_t i = 0; i < fkConvolutedGausNPar; i++)
2442 oldFuncParams[i] = fConvolutedGausDeltaPrime->GetParameter(i);
2443 Int_t oldFuncNpx = fConvolutedGausDeltaPrime->GetNpx();
2444 Double_t oldFuncRangeLow = 0, oldFuncRangeUp = 100;
2445 fConvolutedGausDeltaPrime->GetRange(oldFuncRangeLow, oldFuncRangeUp);
2446
2447 // Choose some sufficiently large range
2448 const Double_t rangeStart = 0.5;
2449 const Double_t rangeEnd = 2.0;
2450
2451 // To get the parameters for the transition, just choose arbitrarily input parameters for mu and sigma
2452 // (it makes sense to choose typical values). The ratio sigma_gauss / sigma_convolution is independent
2453 // of mu and as well the difference mu_gauss - mu_convolution.
2454 Double_t muInput = 1.0;
2455 Double_t sigmaInput = fgkSigmaReferenceForTransitionPars;
2456
2457
2458 // Step 1: Generate distribution with input parameters
2459 const Int_t nPar = 3;
2460 Double_t inputPar[nPar] = { muInput, sigmaInput, lambda };
2461
2462 TH1D* hInput = new TH1D("hInput", "Input distribution", 2000, rangeStart, rangeEnd);
2463
2464 fConvolutedGausDeltaPrime->SetParameters(inputPar);
2465 fConvolutedGausDeltaPrime->SetRange(rangeStart, rangeEnd);
2466 fConvolutedGausDeltaPrime->SetNpx(2000);
2467
2468 /*OLD
2469 // The resolution and mean of the AliPIDResponse are determined in finite intervals
2470 // of ncl (also second order effects due to finite dEdx and tanTheta).
2471 // Take this into account for the transition parameters via assuming an approximately flat
2472 // distritubtion in ncl in this interval.
2473 // NOTE: The ncl interval should be the same as the one used for the sigma map creation!
2474 const Int_t nclStart = 151;
2475 const Int_t nclEnd = 160;
2476 const Int_t nclSteps = (nclEnd - nclStart) + 1;
2477 for (Int_t ncl = nclStart; ncl <= nclEnd; ncl++) {
2478 // Resolution scales with 1/sqrt(ncl)
2479 fConvolutedGausDeltaPrime->SetParameter(1, inputPar[1] * sqrt(nclEnd) / sqrt(ncl));
2480 TH1* hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
2481
2482 for (Int_t i = 0; i < 50000000 / nclSteps; i++)
2483 hInput->Fill(hProbDensity->GetRandom());
2484 }
2485 */
2486
2487 TH1* hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
2488
2489 for (Int_t i = 0; i < 50000000; i++)
2490 hInput->Fill(hProbDensity->GetRandom());
2491
2492 // Step 2: Fit generated distribution with restricted gaussian
2493 Int_t maxBin = hInput->GetMaximumBin();
2494 Double_t max = hInput->GetBinContent(maxBin);
2495
2496 UChar_t usedBins = 1;
2497 if (maxBin > 1) {
2498 max += hInput->GetBinContent(maxBin - 1);
2499 usedBins++;
2500 }
2501 if (maxBin < hInput->GetNbinsX()) {
2502 max += hInput->GetBinContent(maxBin + 1);
2503 usedBins++;
2504 }
2505 max /= usedBins;
2506
2507 // NOTE: The range (<-> fraction of maximum) should be the same
2508 // as for the spline and eta maps creation
2509 const Double_t lowThreshold = hInput->GetXaxis()->GetBinLowEdge(hInput->FindFirstBinAbove(0.1 * max));
2510 const Double_t highThreshold = hInput->GetXaxis()->GetBinUpEdge(hInput->FindLastBinAbove(0.1 * max));
2511
2512 TFitResultPtr fitResGaussFirstStep = hInput->Fit("gaus", "QNRS", "", lowThreshold, highThreshold);
2513
2514 TFitResultPtr fitResGauss;
2515
2516 if ((Int_t)fitResGaussFirstStep == 0) {
2517 TF1 fGauss("fGauss", "[0]*TMath::Gaus(x, [1], [2], 1)", rangeStart, rangeEnd);
2518 fGauss.SetParameter(0, fitResGaussFirstStep->GetParams()[0]);
2519 fGauss.SetParError(0, fitResGaussFirstStep->GetErrors()[0]);
2520 fGauss.SetParameter(2, fitResGaussFirstStep->GetParams()[2]);
2521 fGauss.SetParError(2, fitResGaussFirstStep->GetErrors()[2]);
2522
2523 fGauss.FixParameter(1, fitResGaussFirstStep->GetParams()[1]);
2524 fitResGauss = hInput->Fit(&fGauss, "QNS", "s", rangeStart, rangeEnd);
2525 }
2526 else {
2527 fitResGauss = hInput->Fit("gaus", "QNRS", "same", rangeStart, rangeEnd);
2528 }
2529 //OLD TFitResultPtr fitResGauss = hInput->Fit("gaus", "QNRS", "", hInput->GetXaxis()->GetBinLowEdge(hInput->FindFirstBinAbove(0.1 * max)),
2530 // hInput->GetXaxis()->GetBinUpEdge(hInput->FindLastBinAbove(0.1 * max)));
2531
2532
2533 // Step 3: Use parameters from gaussian fit to obtain parameters for the transition "restricted gauss" -> "convoluted gauss"
2534
2535 // 3.1 The ratio sigmaInput / sigma_gaussFit ONLY depends on lambda (which is fixed per period) -> Calculate this first
2536 // for an arbitrary (here: typical) sigma. The ratio is then ~the same for ALL sigma for given lambda!
2537 if ((Int_t)fitResGauss != 0) {
2538 AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Gauss Fit failed!\n");
2539 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2540 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2541 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2542
2543 delete hInput;
2544 delete oldFuncParams;
2545
2546 return kFALSE;
2547 }
2548
2549 if (fitResGauss->GetParams()[2] <= 0) {
2550 AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Sigma of gauss fit <= 0!\n");
2551 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2552 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2553 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2554
2555 delete hInput;
2556 delete oldFuncParams;
2557
2558 return kFALSE;
2559 }
2560
2561 // sigma correction factor
2562 fConvolutedGaussTransitionPars[1] = sigmaInput / fitResGauss->GetParams()[2];
2563
2564 // 3.2 Now that sigma und lambda are determined, one can calculate mu by shifting the maximum to the desired position,
2565 // i.e. the maximum of the input distribution should coincide with that of the re-calculated distribution,
2566 // which is achieved by getting the same mu for the same sigma.
2567 // NOTE: For fixed lambda, the shift is proportional to sigma and independent of mu!
2568 // So, one can calculate the shift for an arbitrary fixed (here: typical)
2569 // sigma and then simply use this shift for any other sigma by scaling it correspondingly!!!
2570
2571 // Mu shift correction:
2572 // Shift in mu (difference between mean of gauss and mean of convolution) is proportional to sigma!
2573 // Thus, choose a reference sigma (typical -> 0.05), such that for arbitrary sigma one can simple scale
2574 // this shift correction with sigma / referenceSigma.
2575 fConvolutedGaussTransitionPars[0] = (fitResGauss->GetParams()[1] - muInput);
2576
2577
2578 /*Changed on 03.07.2013
2579
2580 // Maximum of fConvolutedGausDeltaPrime should agree with maximum of input
2581 Double_t par[nPar] = { fitResGauss->GetParams()[1], // just as a guess of the maximum position
2582 sigmaInput,
2583 lambda };
2584
2585 fConvolutedGausDeltaPrime->SetParameters(par);
2586
2587 Double_t maxXInput = fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001, muInput - 3. * sigmaInput),
2588 muInput + 10. * sigmaInput,
2589 0.0001);
2590
2591 // Maximum shifts proportional to sigma and is linear in mu (~mean of gauss)
2592 // Maximum should be typically located within [gaussMean, gaussMean + 3 gaussSigma].
2593 // -> Larger search range for safety reasons (also: sigma and/or mean might be not 100% accurate).
2594 Double_t maxXconvoluted = fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001,
2595 fitResGauss->GetParams()[1] - 3. * fitResGauss->GetParams()[2]),
2596 fitResGauss->GetParams()[1] + 10. * fitResGauss->GetParams()[2],
2597 0.0001);
2598 if (maxXconvoluted <= 0) {
2599 AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Maximum of fConvolutedGausDeltaPrime <= 0!\n");
2600
2601 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2602 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2603 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2604
2605 delete hInput;
2606 delete oldFuncParams;
2607
2608 return kFALSE;
2609 }
2610
2611 // maxX perfectly shifts as par[0] (scaled by sigma) -> Can shift maxX to input value.
2612 // Mu shift correction:
2613 fConvolutedGaussTransitionPars[0] = maxXconvoluted - maxXInput;
2614 */
2615
2616
2617
2618 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2619 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2620 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2621
2622 delete hInput;
2623 delete oldFuncParams;
2624
2625 return kTRUE;
2626}
2627
2628
2629//_____________________________________________________________________________
2630AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::SetParamsForConvolutedGaus(const Double_t gausMean, const Double_t gausSigma)
2631{
2632 // Set parameters for convoluted gauss using parameters for a pure gaussian.
2633 // If SetConvolutedGaussLambdaParameter has not been called before to initialise the translation parameters,
2634 // some default parameters will be used and an error will show up.
2635
2636 if (fConvolutedGaussTransitionPars[1] < -998) {
2637 AliError("Transition parameters not initialised! Default parameters will be used. Please call SetConvolutedGaussLambdaParameter(...) before any calculations!");
2638 SetConvolutedGaussLambdaParameter(2.0);
2639 AliError(Form("Parameters set to:\n[0]: %f\n[1]: %f\n[2]: %f\n", fConvolutedGaussTransitionPars[0],
2640 fConvolutedGaussTransitionPars[1], fConvolutedGaussTransitionPars[2]));
2641 }
2642
2643 Double_t par[fkConvolutedGausNPar];
2644 par[2] = fConvolutedGaussTransitionPars[2];
2645 par[1] = fConvolutedGaussTransitionPars[1] * gausSigma;
2646 // maxX perfectly shifts as par[0] (scaled by sigma) -> Can shift maxX so that it sits at the right place.
2647 par[0] = gausMean - fConvolutedGaussTransitionPars[0] * par[1] / fgkSigmaReferenceForTransitionPars;
2648
2649 ErrorCode errCode = kNoErrors;
2650 fConvolutedGausDeltaPrime->SetParameters(par);
2651
2652 fConvolutedGausDeltaPrime->SetNpx(20); // Small value speeds up following algorithm (valid, since extrema far apart)
2653
2654 // Accuracy of 10^-5 is enough to get 0.1% precise peak for MIPS w.r.t. to dEdx = 2000 of protons
2655 // (should boost up the algorithm, because 10^-10 is the default value!)
2656 Double_t maxX= fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001, gausMean - 2. * gausSigma),
2657 gausMean + 6. * gausSigma, 1.0E-5);
2658
2659 const Double_t maximum = fConvolutedGausDeltaPrime->Eval(maxX);
2660 const Double_t maximumFraction = maximum * fAccuracyNonGaussianTail;
2661
2662 // Estimate lower boundary for subsequent search:
2663 Double_t lowBoundSearchBoundLow = TMath::Max(1e-4, maxX - 5. * gausSigma);
2664 Double_t lowBoundSearchBoundUp = maxX;
2665
2666 Bool_t lowerBoundaryFixedAtZero = kFALSE;
2667
2668 while (fConvolutedGausDeltaPrime->Eval(lowBoundSearchBoundLow) >= maximumFraction) {
2669 if (lowBoundSearchBoundLow <= 0) {
2670 // This should only happen to low dEdx particles with very few clusters and therefore large sigma, such that the gauss goes below zero deltaPrime
2671 if (maximum <= 0) { // Something is weired
2672 printf("Error generating signal: maximum is <= 0!\n");
2673 return kError;
2674 }
2675 else {
2676 const Double_t valueAtZero = fConvolutedGausDeltaPrime->Eval(0);
2677 if (valueAtZero / maximum > 0.05) {
2678 // Too large fraction below zero deltaPrime. Signal generation cannot be reliable in this case
2679 printf("Error generating signal: Too large fraction below zero deltaPrime: convGauss(0) / convGauss(max) = %e / %e = %e!\n",
2680 valueAtZero, maximum, valueAtZero / maximum);
2681 return kError;
2682 }
2683 }
2684
2685 /*
2686 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",
2687 fConvolutedGausDeltaPrime->Eval(maxX), fAccuracyNonGaussianTail, maximumFraction, maxX, par[0], par[1], par[2], gausMean, gausSigma);
2688 */
2689
2690 lowerBoundaryFixedAtZero = kTRUE;
2691
2692 if (errCode != kError)
2693 errCode = kWarning;
2694
2695 break;
2696 }
2697
2698 lowBoundSearchBoundUp -= gausSigma;
2699 lowBoundSearchBoundLow -= gausSigma;
2700
2701 if (lowBoundSearchBoundLow < 0) {
2702 lowBoundSearchBoundLow = 0;
2703 lowBoundSearchBoundUp += gausSigma;
2704 }
2705 }
2706
2707 // Determine lower boundary inside estimated range. For small values of the maximum: Need more precision, since finer binning!
2708 Double_t rangeStart = lowerBoundaryFixedAtZero ? 0 :
2709 fConvolutedGausDeltaPrime->GetX(maximumFraction, lowBoundSearchBoundLow, lowBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
2710
2711 // .. and the same for the upper boundary
2712 Double_t rangeEnd = 0;
2713 // If distribution starts beyond upper boundary, everything ends up in the overflow bin. So, just reduce range and Npx to minimum
2714 if (rangeStart > fkDeltaPrimeUpLimit) {
2715 rangeEnd = rangeStart + 0.00001;
2716 fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
2717 fConvolutedGausDeltaPrime->SetNpx(4);
2718 }
2719 else {
2720 // Estimate upper boundary for subsequent search:
2721 Double_t upBoundSearchBoundUp = maxX + 5 * gausSigma;
2722 Double_t upBoundSearchBoundLow = maxX;
2723 while (fConvolutedGausDeltaPrime->Eval(upBoundSearchBoundUp) >= maximumFraction) {
2724 upBoundSearchBoundUp += gausSigma;
2725 upBoundSearchBoundLow += gausSigma;
2726 }
2727
2728 // For small values of the maximum: Need more precision, since finer binning!
2729 rangeEnd = fConvolutedGausDeltaPrime->GetX(maximumFraction, upBoundSearchBoundLow, upBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
2730
2731 fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
2732 fConvolutedGausDeltaPrime->SetNpx(fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindBin(rangeEnd)
2733 - fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindBin(rangeStart) + 1);
2734 //fConvolutedGausDeltaPrime->SetNpx((rangeEnd - rangeStart) / fDeltaPrimeBinWidth + 1);
2735 }
2736
2737 return errCode;
2738}
2739
2740
2741//________________________________________________________________________
2742void AliAnalysisTaskPID::SetUpGenHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const
2743{
2744 // Sets bin limits for axes which are not standard binned and the axes titles.
2745
2746 hist->SetBinEdges(kGenPt, binsPt);
2747 hist->SetBinEdges(kGenDeltaPrimeSpecies, binsDeltaPrime);
2748 hist->SetBinEdges(kGenCentrality, binsCent);
2749
2750 if (fStoreAdditionalJetInformation)
2751 hist->SetBinEdges(kGenJetPt, binsJetPt);
2752
2753 // Set axes titles
2754 hist->GetAxis(kGenMCID)->SetTitle("MC PID");
2755 hist->GetAxis(kGenMCID)->SetBinLabel(1, "e");
2756 hist->GetAxis(kGenMCID)->SetBinLabel(2, "K");
2757 hist->GetAxis(kGenMCID)->SetBinLabel(3, "#mu");
2758 hist->GetAxis(kGenMCID)->SetBinLabel(4, "#pi");
2759 hist->GetAxis(kGenMCID)->SetBinLabel(5, "p");
2760
2761 hist->GetAxis(kGenSelectSpecies)->SetTitle("Select Species");
2762 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(1, "e");
2763 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(2, "K");
2764 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(3, "#pi");
2765 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(4, "p");
2766
2767 hist->GetAxis(kGenPt)->SetTitle("P_{T} (GeV/c)");
2768
2769 hist->GetAxis(kGenDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
2770
2771 hist->GetAxis(kGenCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
2772
2773 if (fStoreAdditionalJetInformation) {
2774 hist->GetAxis(kGenJetPt)->SetTitle("P_{T}^{jet} (GeV/c)");
2775
2776 hist->GetAxis(kGenZ)->SetTitle("z = P_{T}^{track} / P_{T}^{jet}");
2777
2778 hist->GetAxis(kGenXi)->SetTitle("#xi = ln(P_{T}^{jet} / P_{T}^{track})");
2779 }
2780
2781 hist->GetAxis(GetIndexOfChargeAxisGen())->SetTitle("Charge (e_{0})");
2782}
2783
2784
2785//________________________________________________________________________
2786void AliAnalysisTaskPID::SetUpGenYieldHist(THnSparse* hist, Double_t* binsPt, Double_t* binsCent, Double_t* binsJetPt) const
2787{
2788 // Sets bin limits for axes which are not standard binned and the axes titles.
2789
2790 hist->SetBinEdges(kGenYieldPt, binsPt);
2791 hist->SetBinEdges(kGenYieldCentrality, binsCent);
2792 if (fStoreAdditionalJetInformation)
2793 hist->SetBinEdges(kGenYieldJetPt, binsJetPt);
2794
2795 for (Int_t i = 0; i < 5; i++)
2796 hist->GetAxis(kGenYieldMCID)->SetBinLabel(i + 1, AliPID::ParticleLatexName(i));
2797
2798 // Set axes titles
2799 hist->GetAxis(kGenYieldMCID)->SetTitle("MC PID");
2800 hist->GetAxis(kGenYieldPt)->SetTitle("P_{T}^{gen} (GeV/c)");
2801 hist->GetAxis(kGenYieldCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
2802
2803 if (fStoreAdditionalJetInformation) {
2804 hist->GetAxis(kGenYieldJetPt)->SetTitle("P_{T}^{jet, gen} (GeV/c)");
2805
2806 hist->GetAxis(kGenYieldZ)->SetTitle("z = P_{T}^{track} / P_{T}^{jet}");
2807
2808 hist->GetAxis(kGenYieldXi)->SetTitle("#xi = ln(P_{T}^{jet} / P_{T}^{track})");
2809 }
2810
2811 hist->GetAxis(GetIndexOfChargeAxisGenYield())->SetTitle("Charge (e_{0})");
2812}
2813
2814
2815//________________________________________________________________________
2816void AliAnalysisTaskPID::SetUpHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const
2817{
2818 // Sets bin limits for axes which are not standard binned and the axes titles.
2819
2820 hist->SetBinEdges(kDataPt, binsPt);
2821 hist->SetBinEdges(kDataDeltaPrimeSpecies, binsDeltaPrime);
2822 hist->SetBinEdges(kDataCentrality, binsCent);
2823
2824 if (fStoreAdditionalJetInformation)
2825 hist->SetBinEdges(kDataJetPt, binsJetPt);
2826
2827 // Set axes titles
2828 hist->GetAxis(kDataMCID)->SetTitle("MC PID");
2829 hist->GetAxis(kDataMCID)->SetBinLabel(1, "e");
2830 hist->GetAxis(kDataMCID)->SetBinLabel(2, "K");
2831 hist->GetAxis(kDataMCID)->SetBinLabel(3, "#mu");
2832 hist->GetAxis(kDataMCID)->SetBinLabel(4, "#pi");
2833 hist->GetAxis(kDataMCID)->SetBinLabel(5, "p");
2834
2835 hist->GetAxis(kDataSelectSpecies)->SetTitle("Select Species");
2836 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(1, "e");
2837 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(2, "K");
2838 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(3, "#pi");
2839 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(4, "p");
2840
2841 hist->GetAxis(kDataPt)->SetTitle("P_{T} (GeV/c)");
2842
2843 hist->GetAxis(kDataDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
2844
2845 hist->GetAxis(kDataCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
2846
2847 if (fStoreAdditionalJetInformation) {
2848 hist->GetAxis(kDataJetPt)->SetTitle("P_{T}^{jet} (GeV/c)");
2849
2850 hist->GetAxis(kDataZ)->SetTitle("z = P_{T}^{track} / P_{T}^{jet}");
2851
2852 hist->GetAxis(kDataXi)->SetTitle("#xi = ln(P_{T}^{jet} / P_{T}^{track})");
2853 }
2854
2855 hist->GetAxis(GetIndexOfChargeAxisData())->SetTitle("Charge (e_{0})");
2856
2857 /*OLD with TOF, p_TPC_Inner and p_vertex
2858 // MC PID, SelectSpecies, P(TPC_inner), pT, p(Vertex), deltaSpecies, deltaPrimeSpecies, deltaTOFspecies
2859 hist->SetBinEdges(2, binsPt);
2860 hist->SetBinEdges(3, binsPt);
2861 hist->SetBinEdges(4, binsPt);
2862
2863 // Set axes titles
2864 hist->GetAxis(0)->SetTitle("MC PID");
2865 hist->GetAxis(0)->SetBinLabel(1, "e");
2866 hist->GetAxis(0)->SetBinLabel(2, "K");
2867 hist->GetAxis(0)->SetBinLabel(3, "#mu");
2868 hist->GetAxis(0)->SetBinLabel(4, "#pi");
2869 hist->GetAxis(0)->SetBinLabel(5, "p");
2870
2871 hist->GetAxis(1)->SetTitle("Select Species");
2872 hist->GetAxis(1)->SetBinLabel(1, "e");
2873 hist->GetAxis(1)->SetBinLabel(2, "K");
2874 hist->GetAxis(1)->SetBinLabel(3, "#pi");
2875 hist->GetAxis(1)->SetBinLabel(4, "p");
2876
2877 hist->GetAxis(2)->SetTitle("P_{TPC_inner} (GeV/c)");
2878 hist->GetAxis(3)->SetTitle("P_{T} (GeV/c)");
2879 hist->GetAxis(4)->SetTitle("P_{vertex} (GeV/c)");
2880
2881 hist->GetAxis(5)->SetTitle("TPC #Delta_{species} (arb. unit)");
2882
2883 hist->GetAxis(6)->SetTitle("TPC #Delta'_{species} (arb. unit)");
2884
2885 hist->GetAxis(7)->SetTitle("#Delta TOF_{species} (ps)");
2886 */
2887}