7 #include "TFitResultPtr.h"
8 #include "TFitResult.h"
10 #include "AliMCParticle.h"
12 #include "AliAnalysisTask.h"
13 #include "AliAnalysisManager.h"
15 #include "AliESDEvent.h"
16 #include "AliMCEvent.h"
17 #include "AliESDtrackCuts.h"
18 #include "AliAnalysisFilter.h"
19 #include "AliInputEventHandler.h"
21 #include "AliVVertex.h"
23 #include "AliPIDCombined.h"
24 #include "AliPIDResponse.h"
25 #include "AliTPCPIDResponse.h"
27 #include "AliAnalysisTaskPID.h"
30 This task collects PID output from different detectors.
31 Only tracks fulfilling some standard quality cuts are taken into account.
32 At the moment, only data from TPC and TOF is collected. But in future,
33 data from e.g. HMPID is also foreseen.
35 Contact: bhess@cern.ch
38 ClassImp(AliAnalysisTaskPID)
40 const Int_t AliAnalysisTaskPID::fgkNumJetAxes = 3; // Number of additional axes for jets
41 const Double_t AliAnalysisTaskPID::fgkEpsilon = 1e-8; // Double_t threshold above zero
42 const Int_t AliAnalysisTaskPID::fgkMaxNumGenEntries = 500; // Maximum number of generated detector responses per track and delta(Prime) and associated species
44 const Double_t AliAnalysisTaskPID::fgkOneOverSqrt2 = 0.707106781186547462; // = 1. / TMath::Sqrt2();
46 const Double_t AliAnalysisTaskPID::fgkSigmaReferenceForTransitionPars = 0.05; // Reference sigma chosen to calculate transition
48 //________________________________________________________________________
49 AliAnalysisTaskPID::AliAnalysisTaskPID()
50 : AliAnalysisTaskPIDV0base()
52 , fPIDcombined(new AliPIDCombined())
53 , fInputFromOtherTask(kFALSE)
55 , fDoEfficiency(kTRUE)
56 , fDoPtResolution(kFALSE)
57 , fDoDeDxCheck(kFALSE)
58 , fDoBinZeroStudy(kFALSE)
59 , fStoreCentralityPercentile(kFALSE)
60 , fStoreAdditionalJetInformation(kFALSE)
61 , fTakeIntoAccountMuons(kFALSE)
65 , fTPCDefaultPriors(kFALSE)
66 , fUseMCidForGeneration(kTRUE)
67 , fUseConvolutedGaus(kFALSE)
68 , fkConvolutedGausNPar(3)
69 , fAccuracyNonGaussianTail(1e-8)
70 , fkDeltaPrimeLowLimit(0.02)
71 , fkDeltaPrimeUpLimit(40.0)
72 , fConvolutedGausDeltaPrime(0x0)
76 , fPileUpRejectionType(AliAnalysisTaskPIDV0base::kPileUpRejectionOff)
77 , fDoAnySystematicStudiesOnTheExpectedSignal(kFALSE)
78 , fSystematicScalingSplinesThreshold(50.)
79 , fSystematicScalingSplinesBelowThreshold(1.0)
80 , fSystematicScalingSplinesAboveThreshold(1.0)
81 , fSystematicScalingEtaCorrectionMomentumThr(0.35)
82 , fSystematicScalingEtaCorrectionLowMomenta(1.0)
83 , fSystematicScalingEtaCorrectionHighMomenta(1.0)
84 , fSystematicScalingEtaSigmaParaThreshold(250.)
85 , fSystematicScalingEtaSigmaParaBelowThreshold(1.0)
86 , fSystematicScalingEtaSigmaParaAboveThreshold(1.0)
87 , fSystematicScalingMultCorrection(1.0)
88 , fCentralityEstimator("V0M")
95 , fGenRespElDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
96 , fGenRespElDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
97 , fGenRespElDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
98 , fGenRespElDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
99 , fGenRespKaDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
100 , fGenRespKaDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
101 , fGenRespKaDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
102 , fGenRespKaDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
103 , fGenRespPiDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
104 , fGenRespPiDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
105 , fGenRespPiDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
106 , fGenRespPiDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
107 , fGenRespMuDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
108 , fGenRespMuDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
109 , fGenRespMuDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
110 , fGenRespMuDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
111 , fGenRespPrDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
112 , fGenRespPrDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
113 , fGenRespPrDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
114 , fGenRespPrDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
116 , fGenRespElDeltaEl(new Double_t[fgkMaxNumGenEntries])
117 , fGenRespElDeltaKa(new Double_t[fgkMaxNumGenEntries])
118 , fGenRespElDeltaPi(new Double_t[fgkMaxNumGenEntries])
119 , fGenRespElDeltaPr(new Double_t[fgkMaxNumGenEntries])
120 , fGenRespKaDeltaEl(new Double_t[fgkMaxNumGenEntries])
121 , fGenRespKaDeltaKa(new Double_t[fgkMaxNumGenEntries])
122 , fGenRespKaDeltaPi(new Double_t[fgkMaxNumGenEntries])
123 , fGenRespKaDeltaPr(new Double_t[fgkMaxNumGenEntries])
124 , fGenRespPiDeltaEl(new Double_t[fgkMaxNumGenEntries])
125 , fGenRespPiDeltaKa(new Double_t[fgkMaxNumGenEntries])
126 , fGenRespPiDeltaPi(new Double_t[fgkMaxNumGenEntries])
127 , fGenRespPiDeltaPr(new Double_t[fgkMaxNumGenEntries])
128 , fGenRespMuDeltaEl(new Double_t[fgkMaxNumGenEntries])
129 , fGenRespMuDeltaKa(new Double_t[fgkMaxNumGenEntries])
130 , fGenRespMuDeltaPi(new Double_t[fgkMaxNumGenEntries])
131 , fGenRespMuDeltaPr(new Double_t[fgkMaxNumGenEntries])
132 , fGenRespPrDeltaEl(new Double_t[fgkMaxNumGenEntries])
133 , fGenRespPrDeltaKa(new Double_t[fgkMaxNumGenEntries])
134 , fGenRespPrDeltaPi(new Double_t[fgkMaxNumGenEntries])
135 , fGenRespPrDeltaPr(new Double_t[fgkMaxNumGenEntries])
137 , fDeltaPrimeAxis(0x0)
138 , fhMaxEtaVariation(0x0)
139 , fhEventsProcessed(0x0)
140 , fhEventsTriggerSel(0x0)
141 , fhEventsTriggerSelVtxCut(0x0)
142 , fhEventsProcessedNoPileUpRejection(0x0)
143 , fChargedGenPrimariesTriggerSel(0x0)
144 , fChargedGenPrimariesTriggerSelVtxCut(0x0)
145 , fChargedGenPrimariesTriggerSelVtxCutZ(0x0)
146 , fChargedGenPrimariesTriggerSelVtxCutZPileUpRej(0x0)
147 , fhMCgeneratedYieldsPrimaries(0x0)
155 , fOutputContainer(0x0)
158 // default Constructor
160 AliLog::SetClassDebugLevel("AliAnalysisTaskPID", AliLog::kInfo);
162 fConvolutedGausDeltaPrime = new TF1("convolutedGausDeltaPrime", this, &AliAnalysisTaskPID::ConvolutedGaus,
163 fkDeltaPrimeLowLimit, fkDeltaPrimeUpLimit,
164 fkConvolutedGausNPar, "AliAnalysisTaskPID", "ConvolutedGaus");
166 // Set some arbitrary parameteres, such that the function call will not crash
167 // (although it should not be called with these parameters...)
168 fConvolutedGausDeltaPrime->SetParameter(0, 0);
169 fConvolutedGausDeltaPrime->SetParameter(1, 1);
170 fConvolutedGausDeltaPrime->SetParameter(2, 2);
173 // Initialisation of translation parameters is time consuming.
174 // Therefore, default values will only be initialised if they are really needed.
175 // Otherwise, it is left to the user to set the parameter properly.
176 fConvolutedGaussTransitionPars[0] = -999;
177 fConvolutedGaussTransitionPars[1] = -999;
178 fConvolutedGaussTransitionPars[2] = -999;
181 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
182 fFractionHists[i] = 0x0;
183 fFractionSysErrorHists[i] = 0x0;
185 fPtResolution[i] = 0x0;
189 //________________________________________________________________________
190 AliAnalysisTaskPID::AliAnalysisTaskPID(const char *name)
191 : AliAnalysisTaskPIDV0base(name)
193 , fPIDcombined(new AliPIDCombined())
194 , fInputFromOtherTask(kFALSE)
196 , fDoEfficiency(kTRUE)
197 , fDoPtResolution(kFALSE)
198 , fDoDeDxCheck(kFALSE)
199 , fDoBinZeroStudy(kFALSE)
200 , fStoreCentralityPercentile(kFALSE)
201 , fStoreAdditionalJetInformation(kFALSE)
202 , fTakeIntoAccountMuons(kFALSE)
206 , fTPCDefaultPriors(kFALSE)
207 , fUseMCidForGeneration(kTRUE)
208 , fUseConvolutedGaus(kFALSE)
209 , fkConvolutedGausNPar(3)
210 , fAccuracyNonGaussianTail(1e-8)
211 , fkDeltaPrimeLowLimit(0.02)
212 , fkDeltaPrimeUpLimit(40.0)
213 , fConvolutedGausDeltaPrime(0x0)
217 , fPileUpRejectionType(AliAnalysisTaskPIDV0base::kPileUpRejectionOff)
218 , fDoAnySystematicStudiesOnTheExpectedSignal(kFALSE)
219 , fSystematicScalingSplinesThreshold(50.)
220 , fSystematicScalingSplinesBelowThreshold(1.0)
221 , fSystematicScalingSplinesAboveThreshold(1.0)
222 , fSystematicScalingEtaCorrectionMomentumThr(0.35)
223 , fSystematicScalingEtaCorrectionLowMomenta(1.0)
224 , fSystematicScalingEtaCorrectionHighMomenta(1.0)
225 , fSystematicScalingEtaSigmaParaThreshold(250.)
226 , fSystematicScalingEtaSigmaParaBelowThreshold(1.0)
227 , fSystematicScalingEtaSigmaParaAboveThreshold(1.0)
228 , fSystematicScalingMultCorrection(1.0)
229 , fCentralityEstimator("V0M")
236 , fGenRespElDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
237 , fGenRespElDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
238 , fGenRespElDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
239 , fGenRespElDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
240 , fGenRespKaDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
241 , fGenRespKaDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
242 , fGenRespKaDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
243 , fGenRespKaDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
244 , fGenRespPiDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
245 , fGenRespPiDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
246 , fGenRespPiDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
247 , fGenRespPiDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
248 , fGenRespMuDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
249 , fGenRespMuDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
250 , fGenRespMuDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
251 , fGenRespMuDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
252 , fGenRespPrDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
253 , fGenRespPrDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
254 , fGenRespPrDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
255 , fGenRespPrDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
257 , fGenRespElDeltaEl(new Double_t[fgkMaxNumGenEntries])
258 , fGenRespElDeltaKa(new Double_t[fgkMaxNumGenEntries])
259 , fGenRespElDeltaPi(new Double_t[fgkMaxNumGenEntries])
260 , fGenRespElDeltaPr(new Double_t[fgkMaxNumGenEntries])
261 , fGenRespKaDeltaEl(new Double_t[fgkMaxNumGenEntries])
262 , fGenRespKaDeltaKa(new Double_t[fgkMaxNumGenEntries])
263 , fGenRespKaDeltaPi(new Double_t[fgkMaxNumGenEntries])
264 , fGenRespKaDeltaPr(new Double_t[fgkMaxNumGenEntries])
265 , fGenRespPiDeltaEl(new Double_t[fgkMaxNumGenEntries])
266 , fGenRespPiDeltaKa(new Double_t[fgkMaxNumGenEntries])
267 , fGenRespPiDeltaPi(new Double_t[fgkMaxNumGenEntries])
268 , fGenRespPiDeltaPr(new Double_t[fgkMaxNumGenEntries])
269 , fGenRespMuDeltaEl(new Double_t[fgkMaxNumGenEntries])
270 , fGenRespMuDeltaKa(new Double_t[fgkMaxNumGenEntries])
271 , fGenRespMuDeltaPi(new Double_t[fgkMaxNumGenEntries])
272 , fGenRespMuDeltaPr(new Double_t[fgkMaxNumGenEntries])
273 , fGenRespPrDeltaEl(new Double_t[fgkMaxNumGenEntries])
274 , fGenRespPrDeltaKa(new Double_t[fgkMaxNumGenEntries])
275 , fGenRespPrDeltaPi(new Double_t[fgkMaxNumGenEntries])
276 , fGenRespPrDeltaPr(new Double_t[fgkMaxNumGenEntries])
278 , fDeltaPrimeAxis(0x0)
279 , fhMaxEtaVariation(0x0)
280 , fhEventsProcessed(0x0)
281 , fhEventsTriggerSel(0x0)
282 , fhEventsTriggerSelVtxCut(0x0)
283 , fhEventsProcessedNoPileUpRejection(0x0)
284 , fChargedGenPrimariesTriggerSel(0x0)
285 , fChargedGenPrimariesTriggerSelVtxCut(0x0)
286 , fChargedGenPrimariesTriggerSelVtxCutZ(0x0)
287 , fChargedGenPrimariesTriggerSelVtxCutZPileUpRej(0x0)
288 , fhMCgeneratedYieldsPrimaries(0x0)
296 , fOutputContainer(0x0)
301 AliLog::SetClassDebugLevel("AliAnalysisTaskPID", AliLog::kInfo);
303 fConvolutedGausDeltaPrime = new TF1("convolutedGausDeltaPrime", this, &AliAnalysisTaskPID::ConvolutedGaus,
304 fkDeltaPrimeLowLimit, fkDeltaPrimeUpLimit,
305 fkConvolutedGausNPar, "AliAnalysisTaskPID", "ConvolutedGaus");
307 // Set some arbitrary parameteres, such that the function call will not crash
308 // (although it should not be called with these parameters...)
309 fConvolutedGausDeltaPrime->SetParameter(0, 0);
310 fConvolutedGausDeltaPrime->SetParameter(1, 1);
311 fConvolutedGausDeltaPrime->SetParameter(2, 2);
314 // Initialisation of translation parameters is time consuming.
315 // Therefore, default values will only be initialised if they are really needed.
316 // Otherwise, it is left to the user to set the parameter properly.
317 fConvolutedGaussTransitionPars[0] = -999;
318 fConvolutedGaussTransitionPars[1] = -999;
319 fConvolutedGaussTransitionPars[2] = -999;
322 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
323 fFractionHists[i] = 0x0;
324 fFractionSysErrorHists[i] = 0x0;
326 fPtResolution[i] = 0x0;
329 // Define input and output slots here
330 // Input slot #0 works with a TChain
331 DefineInput(0, TChain::Class());
333 DefineOutput(1, TObjArray::Class());
335 DefineOutput(2, AliCFContainer::Class());
337 DefineOutput(3, TObjArray::Class());
341 //________________________________________________________________________
342 AliAnalysisTaskPID::~AliAnalysisTaskPID()
346 CleanupParticleFractionHistos();
348 delete fOutputContainer;
349 fOutputContainer = 0x0;
354 delete fConvolutedGausDeltaPrime;
355 fConvolutedGausDeltaPrime = 0x0;
357 delete fDeltaPrimeAxis;
358 fDeltaPrimeAxis = 0x0;
360 delete [] fGenRespElDeltaPrimeEl;
361 delete [] fGenRespElDeltaPrimeKa;
362 delete [] fGenRespElDeltaPrimePi;
363 delete [] fGenRespElDeltaPrimePr;
365 fGenRespElDeltaPrimeEl = 0x0;
366 fGenRespElDeltaPrimeKa = 0x0;
367 fGenRespElDeltaPrimePi = 0x0;
368 fGenRespElDeltaPrimePr = 0x0;
370 delete [] fGenRespKaDeltaPrimeEl;
371 delete [] fGenRespKaDeltaPrimeKa;
372 delete [] fGenRespKaDeltaPrimePi;
373 delete [] fGenRespKaDeltaPrimePr;
375 fGenRespKaDeltaPrimeEl = 0x0;
376 fGenRespKaDeltaPrimeKa = 0x0;
377 fGenRespKaDeltaPrimePi = 0x0;
378 fGenRespKaDeltaPrimePr = 0x0;
380 delete [] fGenRespPiDeltaPrimeEl;
381 delete [] fGenRespPiDeltaPrimeKa;
382 delete [] fGenRespPiDeltaPrimePi;
383 delete [] fGenRespPiDeltaPrimePr;
385 fGenRespPiDeltaPrimeEl = 0x0;
386 fGenRespPiDeltaPrimeKa = 0x0;
387 fGenRespPiDeltaPrimePi = 0x0;
388 fGenRespPiDeltaPrimePr = 0x0;
390 delete [] fGenRespMuDeltaPrimeEl;
391 delete [] fGenRespMuDeltaPrimeKa;
392 delete [] fGenRespMuDeltaPrimePi;
393 delete [] fGenRespMuDeltaPrimePr;
395 fGenRespMuDeltaPrimeEl = 0x0;
396 fGenRespMuDeltaPrimeKa = 0x0;
397 fGenRespMuDeltaPrimePi = 0x0;
398 fGenRespMuDeltaPrimePr = 0x0;
400 delete [] fGenRespPrDeltaPrimeEl;
401 delete [] fGenRespPrDeltaPrimeKa;
402 delete [] fGenRespPrDeltaPrimePi;
403 delete [] fGenRespPrDeltaPrimePr;
405 fGenRespPrDeltaPrimeEl = 0x0;
406 fGenRespPrDeltaPrimeKa = 0x0;
407 fGenRespPrDeltaPrimePi = 0x0;
408 fGenRespPrDeltaPrimePr = 0x0;
410 delete fhMaxEtaVariation;
411 fhMaxEtaVariation = 0x0;
413 /*OLD with deltaSpecies
414 delete [] fGenRespElDeltaEl;
415 delete [] fGenRespElDeltaKa;
416 delete [] fGenRespElDeltaPi;
417 delete [] fGenRespElDeltaPr;
419 fGenRespElDeltaEl = 0x0;
420 fGenRespElDeltaKa = 0x0;
421 fGenRespElDeltaPi = 0x0;
422 fGenRespElDeltaPr = 0x0;
424 delete [] fGenRespKaDeltaEl;
425 delete [] fGenRespKaDeltaKa;
426 delete [] fGenRespKaDeltaPi;
427 delete [] fGenRespKaDeltaPr;
429 fGenRespKaDeltaEl = 0x0;
430 fGenRespKaDeltaKa = 0x0;
431 fGenRespKaDeltaPi = 0x0;
432 fGenRespKaDeltaPr = 0x0;
434 delete [] fGenRespPiDeltaEl;
435 delete [] fGenRespPiDeltaKa;
436 delete [] fGenRespPiDeltaPi;
437 delete [] fGenRespPiDeltaPr;
439 fGenRespPiDeltaEl = 0x0;
440 fGenRespPiDeltaKa = 0x0;
441 fGenRespPiDeltaPi = 0x0;
442 fGenRespPiDeltaPr = 0x0;
444 delete [] fGenRespMuDeltaEl;
445 delete [] fGenRespMuDeltaKa;
446 delete [] fGenRespMuDeltaPi;
447 delete [] fGenRespMuDeltaPr;
449 fGenRespMuDeltaEl = 0x0;
450 fGenRespMuDeltaKa = 0x0;
451 fGenRespMuDeltaPi = 0x0;
452 fGenRespMuDeltaPr = 0x0;
454 delete [] fGenRespPrDeltaEl;
455 delete [] fGenRespPrDeltaKa;
456 delete [] fGenRespPrDeltaPi;
457 delete [] fGenRespPrDeltaPr;
459 fGenRespPrDeltaEl = 0x0;
460 fGenRespPrDeltaKa = 0x0;
461 fGenRespPrDeltaPi = 0x0;
462 fGenRespPrDeltaPr = 0x0;
467 //________________________________________________________________________
468 void AliAnalysisTaskPID::SetUpPIDcombined()
470 // Initialise the PIDcombined object
472 if (!fDoPID && !fDoDeDxCheck)
476 printf("File: %s, Line: %d: SetUpPIDcombined\n", (char*)__FILE__, __LINE__);
479 AliFatal("No PIDcombined object!\n");
483 fPIDcombined->SetSelectedSpecies(AliPID::kSPECIESC);
484 fPIDcombined->SetEnablePriors(fUsePriors);
486 if (fTPCDefaultPriors)
487 fPIDcombined->SetDefaultTPCPriors();
489 //TODO use individual priors...
491 // Change detector mask (TPC,TOF,ITS)
492 Int_t detectorMask = AliPIDResponse::kDetTPC;
494 // Reject mismatch mask - mismatch only relevant for TOF at the moment - other detectors do not use it
495 Int_t rejectMismatchMask = AliPIDResponse::kDetTPC;
499 detectorMask = detectorMask | AliPIDResponse::kDetITS;
500 rejectMismatchMask = rejectMismatchMask | AliPIDResponse::kDetITS;
503 detectorMask = detectorMask | AliPIDResponse::kDetTOF;
504 rejectMismatchMask = rejectMismatchMask | AliPIDResponse::kDetTOF;
507 fPIDcombined->SetDetectorMask(detectorMask);
508 fPIDcombined->SetRejectMismatchMask(rejectMismatchMask);
511 printf("File: %s, Line: %d: SetUpPIDcombined done\n", (char*)__FILE__, __LINE__);
515 //________________________________________________________________________
516 Bool_t AliAnalysisTaskPID::CalculateMaxEtaVariationMapFromPIDResponse()
518 // Calculate the maximum deviation from unity of the eta correction factors for each row in 1/dEdx(splines)
519 // from the eta correction map of the TPCPIDResponse. The result is stored in fhMaxEtaVariation.
522 AliError("No PID response!");
526 delete fhMaxEtaVariation;
528 const TH2D* hEta = fPIDResponse->GetTPCResponse().GetEtaCorrMap();
530 AliError("No eta correction map!");
534 // Take binning from hEta in Y for fhMaxEtaVariation
535 fhMaxEtaVariation = hEta->ProjectionY("hMaxEtaVariation");
536 fhMaxEtaVariation->SetDirectory(0);
537 fhMaxEtaVariation->Reset();
539 // For each bin in 1/dEdx, loop of all tanTheta bins and find the maximum deviation from unity.
540 // Store the result in fhMaxEtaVariation
542 for (Int_t binY = 1; binY <= fhMaxEtaVariation->GetNbinsX(); binY++) {
543 Double_t maxAbs = -1;
544 for (Int_t binX = 1; binX <= hEta->GetNbinsX(); binX++) {
545 Double_t curr = TMath::Abs(hEta->GetBinContent(binX, binY) - 1.);
550 if (maxAbs < 1e-12) {
551 AliError(Form("Maximum deviation from unity is zero for 1/dEdx = %f (bin %d)", hEta->GetYaxis()->GetBinCenter(binY), binY));
552 delete fhMaxEtaVariation;
556 fhMaxEtaVariation->SetBinContent(binY, maxAbs);
559 printf("AliAnalysisTaskPID: Calculated max eta variation from map \"%s\".\n", hEta->GetTitle());
565 //________________________________________________________________________
566 void AliAnalysisTaskPID::UserCreateOutputObjects()
572 printf("File: %s, Line: %d: UserCreateOutputObjects\n", (char*)__FILE__, __LINE__);
574 // Setup basic things, like PIDResponse
575 AliAnalysisTaskPIDV0base::UserCreateOutputObjects();
578 AliFatal("PIDResponse object was not created");
585 printf("File: %s, Line: %d: UserCreateOutputObjects -> OpenFile(1) successful\n", (char*)__FILE__, __LINE__);
587 fOutputContainer = new TObjArray(1);
588 fOutputContainer->SetName(GetName());
589 fOutputContainer->SetOwner(kTRUE);
591 const Int_t nPtBins = 68;
592 Double_t binsPt[nPtBins+1] = {0. , 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45,
593 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95,
594 1.0, 1.1 , 1.2, 1.3 , 1.4, 1.5 , 1.6, 1.7 , 1.8, 1.9 ,
595 2.0, 2.2 , 2.4, 2.6 , 2.8, 3.0 , 3.2, 3.4 , 3.6, 3.8 ,
596 4.0, 4.5 , 5.0, 5.5 , 6.0, 6.5 , 7.0, 8.0 , 9.0, 10.0,
597 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0,
598 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 };
600 const Int_t nCentBins = 12;
601 //-1 for pp (unless explicitely requested); 90-100 has huge electro-magnetic impurities
602 Double_t binsCent[nCentBins+1] = {-1, 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 };
604 // This centrality estimator deals with integers! This implies that the ranges are always [lowlim, uplim - 1]
605 Double_t binsCentITSTPCTracklets[nCentBins+1] = { 0, 7, 13, 20, 29, 40, 50, 60, 72, 83, 95, 105, 115 };
607 // Special centrality binning for pp
608 Double_t binsCentpp[nCentBins+1] = { 0, 0.01, 0.1, 1, 5, 10, 15, 20, 30, 40, 50, 70, 100};
610 if (fCentralityEstimator.CompareTo("ITSTPCtracklets", TString::kIgnoreCase) == 0 && fStoreCentralityPercentile) {
611 // Special binning for this centrality estimator; but keep number of bins!
612 for (Int_t i = 0; i < nCentBins+1; i++)
613 binsCent[i] = binsCentITSTPCTracklets[i];
615 else if (fCentralityEstimator.Contains("ppMult", TString::kIgnoreCase) && fStoreCentralityPercentile) {
616 // Special binning for this pp centrality estimator; but keep number of bins!
617 for (Int_t i = 0; i < nCentBins+1; i++)
618 binsCent[i] = binsCentpp[i];
621 const Int_t nJetPtBins = 11;
622 Double_t binsJetPt[nJetPtBins+1] = {0, 2, 5, 10, 15, 20, 30, 40, 60, 80, 120, 200};
624 const Int_t nChargeBins = 2;
625 const Double_t binsCharge[nChargeBins+1] = { -1.0 - 1e-4, 0.0, 1.0 + 1e-4 };
627 const Int_t nBinsJets = kDataNumAxes;
628 const Int_t nBinsNoJets = nBinsJets - fgkNumJetAxes;
630 const Int_t nBins = fStoreAdditionalJetInformation ? nBinsJets : nBinsNoJets;
632 // deltaPrimeSpecies binning
633 const Int_t deltaPrimeNBins = 600;
634 Double_t deltaPrimeBins[deltaPrimeNBins + 1];
636 const Double_t fromLow = fkDeltaPrimeLowLimit;
637 const Double_t toHigh = fkDeltaPrimeUpLimit;
638 const Double_t factor = TMath::Power(toHigh/fromLow, 1./deltaPrimeNBins);
640 // Log binning for whole deltaPrime range
641 deltaPrimeBins[0] = fromLow;
642 for (Int_t i = 0 + 1; i <= deltaPrimeNBins; i++) {
643 deltaPrimeBins[i] = factor * deltaPrimeBins[i - 1];
646 fDeltaPrimeAxis = new TAxis(deltaPrimeNBins, deltaPrimeBins);
648 const Int_t nMCPIDbins = 5;
649 const Double_t mcPIDmin = 0.;
650 const Double_t mcPIDmax = 5.;
652 const Int_t nSelSpeciesBins = 4;
653 const Double_t selSpeciesMin = 0.;
654 const Double_t selSpeciesMax = 4.;
656 const Int_t nZBins = 20;
657 const Double_t zMin = 0.;
658 const Double_t zMax = 1.;
660 const Int_t nXiBins = 70;
661 const Double_t xiMin = 0.;
662 const Double_t xiMax = 7.;
664 const Int_t nTOFpidInfoBins = kNumTOFpidInfoBins;
665 const Double_t tofPIDinfoMin = kNoTOFinfo;
666 const Double_t tofPIDinfoMax = kNoTOFinfo + kNumTOFpidInfoBins;
668 // MC PID, SelectSpecies, pT, deltaPrimeSpecies, centrality percentile, jet pT, z = track_pT/jet_pT, xi = log(1/z)
669 Int_t binsNoJets[nBinsNoJets] = { nMCPIDbins,
677 Int_t binsJets[nBinsJets] = { nMCPIDbins,
688 Int_t *bins = fStoreAdditionalJetInformation ? &binsJets[0] : &binsNoJets[0];
690 Double_t xminNoJets[nBinsNoJets] = { mcPIDmin,
698 Double_t xminJets[nBinsJets] = { mcPIDmin,
709 Double_t *xmin = fStoreAdditionalJetInformation? &xminJets[0] : &xminNoJets[0];
711 Double_t xmaxNoJets[nBinsNoJets] = { mcPIDmax,
714 deltaPrimeBins[deltaPrimeNBins],
716 binsCharge[nChargeBins],
719 Double_t xmaxJets[nBinsJets] = { mcPIDmax,
722 deltaPrimeBins[deltaPrimeNBins],
724 binsJetPt[nJetPtBins],
727 binsCharge[nChargeBins],
730 Double_t *xmax = fStoreAdditionalJetInformation? &xmaxJets[0] : &xmaxNoJets[0];
732 fConvolutedGausDeltaPrime->SetNpx(deltaPrimeNBins);
735 fhPIDdataAll = new THnSparseD("hPIDdataAll","", nBins, bins, xmin, xmax);
736 SetUpHist(fhPIDdataAll, binsPt, deltaPrimeBins, binsCent, binsJetPt);
737 fOutputContainer->Add(fhPIDdataAll);
740 // Generated histograms (so far, bins are the same as for primary THnSparse)
741 const Int_t nGenBins = fStoreAdditionalJetInformation ? nBinsJets : nBinsNoJets;
742 // MC PID, SelectSpecies, Pt, deltaPrimeSpecies, jet pT, z = track_pT/jet_pT, xi = log(1/z)
744 Int_t *genBins = fStoreAdditionalJetInformation ? &binsJets[0] : &binsNoJets[0];
745 Double_t *genXmin = fStoreAdditionalJetInformation? &xminJets[0] : &xminNoJets[0];
746 Double_t *genXmax = fStoreAdditionalJetInformation? &xmaxJets[0] : &xmaxNoJets[0];
749 fhGenEl = new THnSparseD("hGenEl", "", nGenBins, genBins, genXmin, genXmax);
750 SetUpGenHist(fhGenEl, binsPt, deltaPrimeBins, binsCent, binsJetPt);
751 fOutputContainer->Add(fhGenEl);
753 fhGenKa = new THnSparseD("hGenKa", "", nGenBins, genBins, genXmin, genXmax);
754 SetUpGenHist(fhGenKa, binsPt, deltaPrimeBins, binsCent, binsJetPt);
755 fOutputContainer->Add(fhGenKa);
757 fhGenPi = new THnSparseD("hGenPi", "", nGenBins, genBins, genXmin, genXmax);
758 SetUpGenHist(fhGenPi, binsPt, deltaPrimeBins, binsCent, binsJetPt);
759 fOutputContainer->Add(fhGenPi);
761 if (fTakeIntoAccountMuons) {
762 fhGenMu = new THnSparseD("hGenMu", "", nGenBins, genBins, genXmin, genXmax);
763 SetUpGenHist(fhGenMu, binsPt, deltaPrimeBins, binsCent, binsJetPt);
764 fOutputContainer->Add(fhGenMu);
767 fhGenPr = new THnSparseD("hGenPr", "", nGenBins, genBins, genXmin, genXmax);
768 SetUpGenHist(fhGenPr, binsPt, deltaPrimeBins, binsCent, binsJetPt);
769 fOutputContainer->Add(fhGenPr);
773 fhEventsProcessed = new TH1D("fhEventsProcessed",
774 "Number of events passing trigger selection, vtx and zvtx cuts and pile-up rejection;Centrality Percentile",
775 nCentBins, binsCent);
776 fhEventsProcessed->Sumw2();
777 fOutputContainer->Add(fhEventsProcessed);
779 fhEventsTriggerSelVtxCut = new TH1D("fhEventsTriggerSelVtxCut",
780 "Number of events passing trigger selection and vtx cut;Centrality Percentile",
781 nCentBins, binsCent);
782 fhEventsTriggerSelVtxCut->Sumw2();
783 fOutputContainer->Add(fhEventsTriggerSelVtxCut);
785 fhEventsTriggerSel = new TH1D("fhEventsTriggerSel",
786 "Number of events passing trigger selection;Centrality Percentile",
787 nCentBins, binsCent);
788 fOutputContainer->Add(fhEventsTriggerSel);
789 fhEventsTriggerSel->Sumw2();
792 fhEventsProcessedNoPileUpRejection = new TH1D("fhEventsProcessedNoPileUpRejection",
793 "Number of events passing trigger selection, vtx and zvtx cuts;Centrality Percentile",
794 nCentBins, binsCent);
795 fOutputContainer->Add(fhEventsProcessedNoPileUpRejection);
796 fhEventsProcessedNoPileUpRejection->Sumw2();
799 // Generated yields within acceptance
800 const Int_t nBinsGenYields = fStoreAdditionalJetInformation ? kGenYieldNumAxes : kGenYieldNumAxes - 3;
801 Int_t genYieldsBins[kGenYieldNumAxes] = { nMCPIDbins, nPtBins, nCentBins, nJetPtBins, nZBins, nXiBins,
803 genYieldsBins[GetIndexOfChargeAxisGenYield()] = nChargeBins;
804 Double_t genYieldsXmin[kGenYieldNumAxes] = { mcPIDmin, binsPt[0], binsCent[0], binsJetPt[0], zMin, xiMin,
806 genYieldsXmin[GetIndexOfChargeAxisGenYield()] = binsCharge[0];
807 Double_t genYieldsXmax[kGenYieldNumAxes] = { mcPIDmax, binsPt[nPtBins], binsCent[nCentBins], binsJetPt[nJetPtBins], zMax, xiMax,
808 binsCharge[nChargeBins] };
809 genYieldsXmax[GetIndexOfChargeAxisGenYield()] = binsCharge[nChargeBins];
812 fhMCgeneratedYieldsPrimaries = new THnSparseD("fhMCgeneratedYieldsPrimaries",
813 "Generated yields w/o reco and cuts inside acceptance (physical primaries)",
814 nBinsGenYields, genYieldsBins, genYieldsXmin, genYieldsXmax);
815 SetUpGenYieldHist(fhMCgeneratedYieldsPrimaries, binsPt, binsCent, binsJetPt);
816 fOutputContainer->Add(fhMCgeneratedYieldsPrimaries);
819 // Container with several process steps (generated and reconstructed level with some variations)
824 printf("File: %s, Line: %d: UserCreateOutputObjects -> OpenFile(2) successful\n", (char*)__FILE__, __LINE__);
826 // Array for the number of bins in each dimension
827 // Dimensions: MC-ID, trackPt, trackEta, trackCharge, cenrality percentile, jetPt, z, xi TODO phi???
828 const Int_t nEffDims = fStoreAdditionalJetInformation ? kEffNumAxes : kEffNumAxes - 3; // Number of dimensions for the efficiency
830 const Int_t nMCIDbins = AliPID::kSPECIES;
831 Double_t binsMCID[nMCIDbins + 1];
833 for(Int_t i = 0; i <= nMCIDbins; i++) {
837 const Int_t nEtaBins = 18;
838 const Double_t binsEta[nEtaBins+1] = {-0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1,
839 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 };
841 const Int_t nEffBins[kEffNumAxes] = { nMCIDbins, nPtBins, nEtaBins, nChargeBins, nCentBins, nJetPtBins, nZBins, nXiBins };
843 fContainerEff = new AliCFContainer("containerEff", "Reconstruction Efficiency x Acceptance x Resolution and Secondary Correction",
844 kNumSteps, nEffDims, nEffBins);
846 // Setting the bin limits
847 fContainerEff->SetBinLimits(kEffMCID, binsMCID);
848 fContainerEff->SetBinLimits(kEffTrackPt, binsPt);
849 fContainerEff->SetBinLimits(kEffTrackEta, binsEta);
850 fContainerEff->SetBinLimits(kEffTrackCharge, binsCharge);
851 fContainerEff->SetBinLimits(kEffCentrality, binsCent);
852 if (fStoreAdditionalJetInformation) {
853 fContainerEff->SetBinLimits(kEffJetPt, binsJetPt);
854 fContainerEff->SetBinLimits(kEffZ, zMin, zMax);
855 fContainerEff->SetBinLimits(kEffXi, xiMin, xiMax);
858 fContainerEff->SetVarTitle(kEffMCID,"MC ID");
859 fContainerEff->SetVarTitle(kEffTrackPt,"p_{T} (GeV/c)");
860 fContainerEff->SetVarTitle(kEffTrackEta,"#eta");
861 fContainerEff->SetVarTitle(kEffTrackCharge,"Charge (e_{0})");
862 fContainerEff->SetVarTitle(kEffCentrality, "Centrality Percentile");
863 if (fStoreAdditionalJetInformation) {
864 fContainerEff->SetVarTitle(kEffJetPt, "p_{T}^{jet} (GeV/c)");
865 fContainerEff->SetVarTitle(kEffZ, "z = p_{T}^{track} / p_{T}^{jet}");
866 fContainerEff->SetVarTitle(kEffXi, "#xi = ln(p_{T}^{jet} / p_{T}^{track})");
869 // Define clean MC sample
870 fContainerEff->SetStepTitle(kStepGenWithGenCuts, "Particle level, cuts on particle level");
871 // For Acceptance x Efficiency correction of primaries
872 fContainerEff->SetStepTitle(kStepRecWithGenCuts, "Detector level (rec) with cuts on particle level");
873 // For (pT) resolution correction
874 fContainerEff->SetStepTitle(kStepRecWithGenCutsMeasuredObs,
875 "Detector level (rec) with cuts on particle level with measured observables");
876 // For secondary correction
877 fContainerEff->SetStepTitle(kStepRecWithRecCutsMeasuredObs,
878 "Detector level, all cuts on detector level with measured observables");
879 fContainerEff->SetStepTitle(kStepRecWithRecCutsPrimaries,
880 "Detector level, all cuts on detector level, only MC primaries");
881 fContainerEff->SetStepTitle(kStepRecWithRecCutsMeasuredObsPrimaries,
882 "Detector level, all cuts on detector level with measured observables, only MC primaries");
883 fContainerEff->SetStepTitle(kStepRecWithRecCutsMeasuredObsStrangenessScaled,
884 "Detector level (strangeness scaled), all cuts on detector level with measured observables");
887 if (fDoPID || fDoEfficiency) {
889 fh2FFJetPtRec = new TH2D("fh2FFJetPtRec", "Number of reconstructed jets;Centrality Percentile;p_{T}^{jet} (GeV/c)",
890 nCentBins, binsCent, nJetPtBins, binsJetPt);
891 fh2FFJetPtRec->Sumw2();
892 fOutputContainer->Add(fh2FFJetPtRec);
893 fh2FFJetPtGen = new TH2D("fh2FFJetPtGen", "Number of generated jets;Centrality Percentile;p_{T}^{jet} (GeV/c)",
894 nCentBins, binsCent, nJetPtBins, binsJetPt);
895 fh2FFJetPtGen->Sumw2();
896 fOutputContainer->Add(fh2FFJetPtGen);
899 // Pythia information
900 fh1Xsec = new TProfile("fh1Xsec", "xsec from pyxsec.root", 1, 0, 1);
902 fh1Xsec->GetXaxis()->SetBinLabel(1, "<#sigma>");
903 fh1Trials = new TH1D("fh1Trials", "trials from pyxsec.root", 1, 0, 1);
905 fh1Trials->GetXaxis()->SetBinLabel(1, "#sum{ntrials}");
907 fOutputContainer->Add(fh1Xsec);
908 fOutputContainer->Add(fh1Trials);
910 if (fDoDeDxCheck || fDoPtResolution) {
912 fQAContainer = new TObjArray(1);
913 fQAContainer->SetName(Form("%s_QA", GetName()));
914 fQAContainer->SetOwner(kTRUE);
917 printf("File: %s, Line: %d: UserCreateOutputObjects -> OpenFile(3) successful\n", (char*)__FILE__, __LINE__);
920 if (fDoPtResolution) {
921 const Int_t nPtBinsRes = 100;
922 Double_t pTbinsRes[nPtBinsRes + 1];
924 const Double_t fromLowPtRes = 0.15;
925 const Double_t toHighPtRes = 50.;
926 const Double_t factorPtRes = TMath::Power(toHighPtRes/fromLowPtRes, 1./nPtBinsRes);
927 // Log binning for whole pT range
928 pTbinsRes[0] = fromLowPtRes;
929 for (Int_t i = 0 + 1; i <= nPtBinsRes; i++) {
930 pTbinsRes[i] = factorPtRes * pTbinsRes[i - 1];
933 const Int_t nBinsPtResolution = kPtResNumAxes;
934 Int_t ptResolutionBins[kPtResNumAxes] = { nJetPtBins, nPtBinsRes, nPtBinsRes,
935 nChargeBins, nCentBins };
936 Double_t ptResolutionXmin[kPtResNumAxes] = { binsJetPt[0], pTbinsRes[0], pTbinsRes[0],
937 binsCharge[0], binsCent[0] };
938 Double_t ptResolutionXmax[kPtResNumAxes] = { binsJetPt[nJetPtBins], pTbinsRes[nPtBinsRes], pTbinsRes[nPtBinsRes],
939 binsCharge[nChargeBins], binsCent[nCentBins] };
941 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
942 fPtResolution[i] = new THnSparseD(Form("fPtResolution_%s", AliPID::ParticleShortName(i)),
943 Form("Pt resolution for primaries, %s", AliPID::ParticleLatexName(i)),
944 nBinsPtResolution, ptResolutionBins, ptResolutionXmin, ptResolutionXmax);
945 SetUpPtResHist(fPtResolution[i], pTbinsRes, binsJetPt, binsCent);
946 fQAContainer->Add(fPtResolution[i]);
950 // Besides the pT resolution, also perform check on shared clusters
951 const Int_t nBinsQASharedCls = kQASharedClsNumAxes;
952 Int_t qaSharedClsBins[kQASharedClsNumAxes] = { nJetPtBins, nPtBinsRes, 160, 160 };
953 Double_t qaSharedClsXmin[kQASharedClsNumAxes] = { binsJetPt[0], pTbinsRes[0], 0, -1 };
954 Double_t qaSharedClsXmax[kQASharedClsNumAxes] = { binsJetPt[nJetPtBins], pTbinsRes[nPtBinsRes], 160, 159 };
956 fQASharedCls = new THnSparseD("fQASharedCls", "QA shared clusters", nBinsQASharedCls, qaSharedClsBins, qaSharedClsXmin, qaSharedClsXmax);
958 SetUpSharedClsHist(fQASharedCls, pTbinsRes, binsJetPt);
959 fQAContainer->Add(fQASharedCls);
965 const Int_t nEtaAbsBins = 9;
966 const Double_t binsEtaAbs[nEtaAbsBins+1] = { 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 };
968 const Double_t dEdxMin = 20;
969 const Double_t dEdxMax = 110;
970 const Int_t nDeDxBins = (Int_t) ((dEdxMax - dEdxMin) / 0.02);
971 const Int_t nBinsDeDxCheck= kDeDxCheckNumAxes;
972 Int_t dEdxCheckBins[kDeDxCheckNumAxes] = { nSelSpeciesBins, nPtBins, nJetPtBins, nEtaAbsBins, nDeDxBins };
973 Double_t dEdxCheckXmin[kDeDxCheckNumAxes] = { selSpeciesMin, binsPt[0], binsJetPt[0], binsEtaAbs[0], dEdxMin };
974 Double_t dEdxCheckXmax[kDeDxCheckNumAxes] = { selSpeciesMax, binsPt[nPtBins], binsJetPt[nJetPtBins], binsEtaAbs[nEtaAbsBins], dEdxMax };
976 fDeDxCheck = new THnSparseD("fDeDxCheck", "dEdx check", nBinsDeDxCheck, dEdxCheckBins, dEdxCheckXmin, dEdxCheckXmax);
977 SetUpDeDxCheckHist(fDeDxCheck, binsPt, binsJetPt, binsEtaAbs);
978 fQAContainer->Add(fDeDxCheck);
981 if (fDoBinZeroStudy) {
982 const Double_t etaLow = -0.9;
983 const Double_t etaUp = 0.9;
984 const Int_t nEtaBins = 18;
986 const Int_t nBinsBinZeroStudy = kBinZeroStudyNumAxes;
987 Int_t binZeroStudyBins[nBinsBinZeroStudy] = { nCentBins, nPtBins, nEtaBins };
988 Double_t binZeroStudyXmin[nBinsBinZeroStudy] = { binsCent[0], binsPt[0], etaLow };
989 Double_t binZeroStudyXmax[nBinsBinZeroStudy] = { binsCent[nCentBins], binsPt[nPtBins], etaUp };
991 fChargedGenPrimariesTriggerSel = new THnSparseD("fChargedGenPrimariesTriggerSel", "Trigger sel.", nBinsBinZeroStudy, binZeroStudyBins,
992 binZeroStudyXmin, binZeroStudyXmax);
993 SetUpBinZeroStudyHist(fChargedGenPrimariesTriggerSel, binsCent, binsPt);
994 fOutputContainer->Add(fChargedGenPrimariesTriggerSel);
996 fChargedGenPrimariesTriggerSelVtxCut = new THnSparseD("fChargedGenPrimariesTriggerSelVtxCut", "Vertex cut", nBinsBinZeroStudy,
997 binZeroStudyBins, binZeroStudyXmin, binZeroStudyXmax);
998 SetUpBinZeroStudyHist(fChargedGenPrimariesTriggerSelVtxCut, binsCent, binsPt);
999 fOutputContainer->Add(fChargedGenPrimariesTriggerSelVtxCut);
1001 fChargedGenPrimariesTriggerSelVtxCutZ = new THnSparseD("fChargedGenPrimariesTriggerSelVtxCutZ", "Vertex #it{z} cut", nBinsBinZeroStudy,
1002 binZeroStudyBins, binZeroStudyXmin, binZeroStudyXmax);
1003 SetUpBinZeroStudyHist(fChargedGenPrimariesTriggerSelVtxCutZ, binsCent, binsPt);
1004 fOutputContainer->Add(fChargedGenPrimariesTriggerSelVtxCutZ);
1006 fChargedGenPrimariesTriggerSelVtxCutZPileUpRej = new THnSparseD("fChargedGenPrimariesTriggerSelVtxCutZPileUpRej", "Vertex #it{z} cut",
1007 nBinsBinZeroStudy, binZeroStudyBins, binZeroStudyXmin, binZeroStudyXmax);
1008 SetUpBinZeroStudyHist(fChargedGenPrimariesTriggerSelVtxCutZPileUpRej, binsCent, binsPt);
1009 fOutputContainer->Add(fChargedGenPrimariesTriggerSelVtxCutZPileUpRej);
1013 printf("File: %s, Line: %d: UserCreateOutputObjects -> Posting output data\n", (char*)__FILE__, __LINE__);
1015 PostData(1, fOutputContainer);
1017 PostData(2, fContainerEff);
1018 if (fDoDeDxCheck || fDoPtResolution)
1019 PostData(3, fQAContainer);
1022 printf("File: %s, Line: %d: UserCreateOutputObjects -> Done\n", (char*)__FILE__, __LINE__);
1026 //________________________________________________________________________
1027 void AliAnalysisTaskPID::UserExec(Option_t *)
1030 // Called for each event
1033 printf("File: %s, Line: %d: UserExec\n", (char*)__FILE__, __LINE__);
1035 // No processing of event, if input is fed in directly from another task
1036 if (fInputFromOtherTask)
1040 printf("File: %s, Line: %d: UserExec -> Processing started\n", (char*)__FILE__, __LINE__);
1042 fEvent = dynamic_cast<AliVEvent*>(InputEvent());
1044 Printf("ERROR: fEvent not available");
1048 ConfigureTaskForCurrentEvent(fEvent);
1050 fMC = dynamic_cast<AliMCEvent*>(MCEvent());
1052 if (!fPIDResponse || !fPIDcombined)
1055 Double_t centralityPercentile = -1;
1056 if (fStoreCentralityPercentile) {
1057 if (fCentralityEstimator.CompareTo("ITSTPCtracklets", TString::kIgnoreCase) == 0) {
1058 // Special pp centrality estimator
1059 AliESDEvent* esdEvent = dynamic_cast<AliESDEvent*>(fEvent);
1061 AliError("Not esd event -> Cannot use tracklet multiplicity estimator!");
1062 centralityPercentile = -1;
1065 centralityPercentile = AliESDtrackCuts::GetReferenceMultiplicity(esdEvent, AliESDtrackCuts::kTrackletsITSTPC, fEtaAbsCutUp);
1067 else if (fCentralityEstimator.Contains("ppMult", TString::kIgnoreCase)) {
1068 // Another special pp centrality estimator
1069 centralityPercentile = fAnaUtils->GetMultiplicityPercentile(fEvent, GetPPCentralityEstimator().Data());
1072 // Ordinary centrality estimator
1073 centralityPercentile = fEvent->GetCentrality()->GetCentralityPercentile(fCentralityEstimator.Data());
1077 // Check if vertex is ok, but don't apply cut on z position
1078 const Bool_t passedVertexSelection = GetVertexIsOk(fEvent, kFALSE);
1079 // Now check again, but also require z position to be in desired range
1080 const Bool_t passedVertexZSelection = GetVertexIsOk(fEvent, kTRUE);
1082 const Bool_t isPileUp = GetIsPileUp(fEvent, fPileUpRejectionType);
1086 if (fDoBinZeroStudy && fMC) {
1087 for (Int_t iPart = 0; iPart < fMC->GetNumberOfTracks(); iPart++) {
1088 AliMCParticle *mcPart = dynamic_cast<AliMCParticle*>(fMC->GetTrack(iPart));
1093 if (!fMC->IsPhysicalPrimary(iPart))
1096 const Double_t etaGen = mcPart->Eta();
1097 const Double_t ptGen = mcPart->Pt();
1099 Double_t values[kBinZeroStudyNumAxes] = { 0. };
1100 values[kBinZeroStudyCentrality] = centralityPercentile;
1101 values[kBinZeroStudyGenPt] = ptGen;
1102 values[kBinZeroStudyGenEta] = etaGen;
1104 fChargedGenPrimariesTriggerSel->Fill(values);
1105 if (passedVertexSelection) {
1106 fChargedGenPrimariesTriggerSelVtxCut->Fill(values);
1107 if (passedVertexZSelection) {
1108 fChargedGenPrimariesTriggerSelVtxCutZ->Fill(values);
1110 fChargedGenPrimariesTriggerSelVtxCutZPileUpRej->Fill(values);
1119 // Event counters for trigger selection, vertex cuts and pile-up rejection
1120 IncrementEventCounter(centralityPercentile, kTriggerSel);
1122 if (!passedVertexSelection)
1125 IncrementEventCounter(centralityPercentile, kTriggerSelAndVtxCut);
1127 if (!passedVertexZSelection)
1130 IncrementEventCounter(centralityPercentile, kTriggerSelAndVtxCutAndZvtxCutNoPileUpRejection);
1131 // ATTENTION: Is this the right place for the pile-up rejection? Important to have still the proper bin-0 correction,
1132 // which is done solely with sel and selVtx, since the zvtx selection does ~not change the spectra. The question is whether the pile-up
1133 // rejection changes the spectra. If not, then it is perfectly fine to put it here and keep the usual histo for the normalisation to number
1134 // of events. But if it does change the spectra, this must somehow be corrected for.
1135 // NOTE: multiplicity >= 0 usually implies a properly reconstructed vertex. Hence, the bin-0 correction cannot be done in multiplicity bins.
1136 // Furthermore, there seems to be no MC simulation with pile-up rejection, so the bin-0 correction cannot be extracted with it. Pile-up
1137 // rejection has only a minor impact, so maybe there is no need to dig further.
1141 IncrementEventCounter(centralityPercentile, kTriggerSelAndVtxCutAndZvtxCut);
1143 Double_t magField = fEvent->GetMagneticField();
1146 if (fDoPID || fDoEfficiency) {
1147 for (Int_t iPart = 0; iPart < fMC->GetNumberOfTracks(); iPart++) {
1148 AliMCParticle *mcPart = dynamic_cast<AliMCParticle*>(fMC->GetTrack(iPart));
1153 // Define clean MC sample with corresponding particle level track cuts:
1154 // - MC-track must be in desired eta range
1155 // - MC-track must be physical primary
1156 // - Species must be one of those in question (everything else goes to the overflow bin of mcID)
1158 // Geometrie should be the same as on the reconstructed level -> By definition analysis within this eta interval
1159 if (!IsInAcceptedEtaRange(TMath::Abs(mcPart->Eta()))) continue;
1161 Int_t mcID = PDGtoMCID(mcPart->PdgCode());
1163 // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
1164 Double_t chargeMC = mcPart->Charge() / 3.;
1166 if (TMath::Abs(chargeMC) < 0.01)
1167 continue; // Reject neutral particles (only relevant, if mcID is not used)
1169 if (!fMC->IsPhysicalPrimary(iPart))
1173 Double_t valuesGenYield[kGenYieldNumAxes] = { static_cast<Double_t>(mcID), mcPart->Pt(), centralityPercentile, -1, -1, -1, -1 };
1174 valuesGenYield[GetIndexOfChargeAxisGenYield()] = chargeMC;
1176 fhMCgeneratedYieldsPrimaries->Fill(valuesGenYield);
1180 if (fDoEfficiency) {
1181 Double_t valueEff[kEffNumAxes] = { static_cast<Double_t>(mcID), mcPart->Pt(), mcPart->Eta(), chargeMC, centralityPercentile,
1184 fContainerEff->Fill(valueEff, kStepGenWithGenCuts);
1190 // Track loop to fill a Train spectrum
1191 for (Int_t iTracks = 0; iTracks < fEvent->GetNumberOfTracks(); iTracks++) {
1192 AliVTrack* track = dynamic_cast<AliVTrack*>(fEvent->GetTrack(iTracks));
1194 Printf("ERROR: Could not retrieve track %d", iTracks);
1199 // Apply detector level track cuts
1200 const Bool_t tuneOnDataTPC = fPIDResponse->IsTunedOnData() &&
1201 ((fPIDResponse->GetTunedOnDataMask() & AliPIDResponse::kDetTPC) == AliPIDResponse::kDetTPC);
1202 Double_t dEdxTPC = tuneOnDataTPC ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
1206 if(fTrackFilter && !fTrackFilter->IsSelected(track))
1209 if (GetUseTPCCutMIGeo()) {
1210 if (!TPCCutMIGeo(track, fEvent))
1213 else if (GetUseTPCnclCut()) {
1214 if (!TPCnclCut(track))
1219 if (!PhiPrimeCut(track, magField))
1220 continue; // reject track
1223 Int_t pdg = 0; // = 0 indicates data for the moment
1224 AliMCParticle* mcTrack = 0x0;
1225 Int_t mcID = AliPID::kUnknown;
1229 label = track->GetLabel();
1234 mcTrack = dynamic_cast<AliMCParticle*>(fMC->GetTrack(TMath::Abs(label)));
1236 Printf("ERROR: Could not retrieve mcTrack with label %d for track %d", label, iTracks);
1240 pdg = mcTrack->PdgCode();
1241 mcID = PDGtoMCID(mcTrack->PdgCode());
1243 if (fDoEfficiency) {
1244 // For efficiency: Reconstructed track has survived all cuts on the detector level (excluding acceptance)
1245 // and has an associated MC track which is a physical primary and was generated inside the acceptance
1246 if (fMC->IsPhysicalPrimary(TMath::Abs(label)) &&
1247 IsInAcceptedEtaRange(TMath::Abs(mcTrack->Eta()))) {
1249 // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
1250 Double_t value[kEffNumAxes] = { static_cast<Double_t>(mcID), mcTrack->Pt(), mcTrack->Eta(), mcTrack->Charge() / 3., centralityPercentile,
1252 fContainerEff->Fill(value, kStepRecWithGenCuts);
1254 Double_t valueMeas[kEffNumAxes] = { static_cast<Double_t>(mcID), track->Pt(), track->Eta(), static_cast<Double_t>(track->Charge()), centralityPercentile,
1256 fContainerEff->Fill(valueMeas, kStepRecWithGenCutsMeasuredObs);
1261 // Only process tracks inside the desired eta window
1262 if (!IsInAcceptedEtaRange(TMath::Abs(track->Eta()))) continue;
1264 if (fDoPID || fDoDeDxCheck || fDoPtResolution)
1265 ProcessTrack(track, pdg, centralityPercentile, -1); // No jet information in this case -> Set jet pT to -1
1267 if (fDoPtResolution) {
1268 if (mcTrack && fMC->IsPhysicalPrimary(TMath::Abs(label))) {
1269 // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
1270 Double_t valuePtRes[kPtResNumAxes] = { -1, mcTrack->Pt(), track->Pt(), mcTrack->Charge() / 3., centralityPercentile };
1271 FillPtResolution(mcID, valuePtRes);
1275 if (fDoEfficiency) {
1277 Double_t valueRecAllCuts[kEffNumAxes] = { static_cast<Double_t>(mcID), track->Pt(), track->Eta(), static_cast<Double_t>(track->Charge()), centralityPercentile,
1279 fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObs);
1281 Double_t weight = IsSecondaryWithStrangeMotherMC(fMC, TMath::Abs(label)) ?
1282 GetMCStrangenessFactorCMS(fMC, mcTrack) : 1.0;
1283 fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObsStrangenessScaled, weight);
1285 // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
1286 Double_t valueGenAllCuts[kEffNumAxes] = { static_cast<Double_t>(mcID), mcTrack->Pt(), mcTrack->Eta(), mcTrack->Charge() / 3.,
1287 centralityPercentile, -1, -1, -1 };
1288 if (fMC->IsPhysicalPrimary(TMath::Abs(label))) {
1289 fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObsPrimaries);
1290 fContainerEff->Fill(valueGenAllCuts, kStepRecWithRecCutsPrimaries);
1297 printf("File: %s, Line: %d: UserExec -> Processing done\n", (char*)__FILE__, __LINE__);
1302 printf("File: %s, Line: %d: UserExec -> Done\n", (char*)__FILE__, __LINE__);
1305 //________________________________________________________________________
1306 void AliAnalysisTaskPID::Terminate(const Option_t *)
1308 // Draw result to the screen
1309 // Called once at the end of the query
1313 //_____________________________________________________________________________
1314 void AliAnalysisTaskPID::CheckDoAnyStematicStudiesOnTheExpectedSignal()
1316 // Check whether at least one scale factor indicates the ussage of systematic studies
1317 // and set internal flag accordingly.
1319 fDoAnySystematicStudiesOnTheExpectedSignal = kFALSE;
1322 if ((TMath::Abs(fSystematicScalingSplinesBelowThreshold - 1.0) > fgkEpsilon) ||
1323 (TMath::Abs(fSystematicScalingSplinesAboveThreshold - 1.0) > fgkEpsilon)) {
1324 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1328 if ((TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon) ||
1329 (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon)) {
1330 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1334 if ((TMath::Abs(fSystematicScalingEtaSigmaParaBelowThreshold - 1.0) > fgkEpsilon) ||
1335 (TMath::Abs(fSystematicScalingEtaSigmaParaAboveThreshold - 1.0) > fgkEpsilon)) {
1336 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1340 if (TMath::Abs(fSystematicScalingMultCorrection - 1.0) > fgkEpsilon) {
1341 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1347 //_____________________________________________________________________________
1348 void AliAnalysisTaskPID::ConfigureTaskForCurrentEvent(AliVEvent* event)
1350 // Configure the task for the current event. In particular, this is needed if the run number changes
1353 AliError("Could not set up task: no event!");
1357 Int_t run = event->GetRunNumber();
1360 // If systematics on eta is investigated, need to calculate the maxEtaVariationMap
1361 if ((TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon) ||
1362 (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon)) {
1363 if (!CalculateMaxEtaVariationMapFromPIDResponse())
1364 AliFatal("Systematics on eta correction requested, but failed to calculate max eta varation map!");
1372 //_____________________________________________________________________________
1373 Int_t AliAnalysisTaskPID::PDGtoMCID(Int_t pdg)
1375 // Returns the corresponding AliPID index to the given pdg code.
1376 // Returns AliPID::kUnkown if pdg belongs to a not considered species.
1378 Int_t absPDGcode = TMath::Abs(pdg);
1379 if (absPDGcode == 211) {//Pion
1380 return AliPID::kPion;
1382 else if (absPDGcode == 321) {//Kaon
1383 return AliPID::kKaon;
1385 else if (absPDGcode == 2212) {//Proton
1386 return AliPID::kProton;
1388 else if (absPDGcode == 11) {//Electron
1389 return AliPID::kElectron;
1391 else if (absPDGcode == 13) {//Muon
1392 return AliPID::kMuon;
1395 return AliPID::kUnknown;
1399 //_____________________________________________________________________________
1400 void AliAnalysisTaskPID::GetJetTrackObservables(Double_t trackPt, Double_t jetPt, Double_t& z, Double_t& xi)
1402 // Uses trackPt and jetPt to obtain z and xi.
1404 z = (jetPt > 0 && trackPt >= 0) ? (trackPt / jetPt) : -1;
1405 xi = (z > 0) ? TMath::Log(1. / z) : -1;
1407 if(trackPt > (1. - 1e-06) * jetPt && trackPt < (1. + 1e-06) * jetPt) { // case z=1 : move entry to last histo bin <1
1414 //_____________________________________________________________________________
1415 void AliAnalysisTaskPID::CleanupParticleFractionHistos()
1417 // Delete histos with particle fractions
1419 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1420 delete fFractionHists[i];
1421 fFractionHists[i] = 0x0;
1423 delete fFractionSysErrorHists[i];
1424 fFractionSysErrorHists[i] = 0x0;
1429 //_____________________________________________________________________________
1430 Double_t AliAnalysisTaskPID::ConvolutedGaus(const Double_t* xx, const Double_t* par) const
1432 // Convolutes gauss with an exponential tail which describes dEdx-response better than pure gaussian
1434 const Double_t mean = par[0];
1435 const Double_t sigma = par[1];
1436 const Double_t lambda = par[2];
1439 printf("File: %s, Line: %d: ConvolutedGaus: mean %e, sigma %e, lambda %e\n", (char*)__FILE__, __LINE__, mean, sigma, lambda);
1441 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);
1445 //_____________________________________________________________________________
1446 inline Double_t AliAnalysisTaskPID::FastGaus(Double_t x, Double_t mean, Double_t sigma) const
1448 // Calculate an unnormalised gaussian function with mean and sigma.
1450 if (sigma < fgkEpsilon)
1453 const Double_t arg = (x - mean) / sigma;
1454 return exp(-0.5 * arg * arg);
1458 //_____________________________________________________________________________
1459 inline Double_t AliAnalysisTaskPID::FastNormalisedGaus(Double_t x, Double_t mean, Double_t sigma) const
1461 // Calculate a normalised (divided by sqrt(2*Pi)*sigma) gaussian function with mean and sigma.
1463 if (sigma < fgkEpsilon)
1466 const Double_t arg = (x - mean) / sigma;
1467 const Double_t res = exp(-0.5 * arg * arg);
1468 return res / (2.50662827463100024 * sigma); //sqrt(2*Pi)=2.50662827463100024
1472 //_____________________________________________________________________________
1473 Int_t AliAnalysisTaskPID::FindBinWithinRange(TAxis* axis, Double_t value) const
1475 // Find the corresponding bin of the axis. Values outside the range (also under and overflow) will be set to the first/last
1478 Int_t bin = axis->FindFixBin(value);
1482 if (bin > axis->GetNbins())
1483 bin = axis->GetNbins();
1489 //_____________________________________________________________________________
1490 Int_t AliAnalysisTaskPID::FindFirstBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yBin,
1493 // Kind of projects a TH3 to 1 bin combination in y and z
1494 // and looks for the first x bin above a threshold for this projection.
1495 // If no such bin is found, -1 is returned.
1500 Int_t nBinsX = hist->GetNbinsX();
1501 for (Int_t xBin = 1; xBin <= nBinsX; xBin++) {
1502 if (hist->GetBinContent(xBin, yBin, zBin) > threshold)
1510 //_____________________________________________________________________________
1511 Int_t AliAnalysisTaskPID::FindLastBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yBin,
1514 // Kind of projects a TH3 to 1 bin combination in y and z
1515 // and looks for the last x bin above a threshold for this projection.
1516 // If no such bin is found, -1 is returned.
1521 Int_t nBinsX = hist->GetNbinsX();
1522 for (Int_t xBin = nBinsX; xBin >= 1; xBin--) {
1523 if (hist->GetBinContent(xBin, yBin, zBin) > threshold)
1531 //_____________________________________________________________________________
1532 Bool_t AliAnalysisTaskPID::GetParticleFraction(Double_t trackPt, Double_t jetPt, Double_t centralityPercentile,
1533 AliPID::EParticleType species,
1534 Double_t& fraction, Double_t& fractionErrorStat, Double_t& fractionErrorSys) const
1536 // Computes the particle fraction for the corresponding species for the given trackPt, jetPt and centrality.
1537 // Use jetPt = -1 for inclusive spectra and centralityPercentile = -1 for pp.
1538 // On success (return value kTRUE), fraction contains the particle fraction, fractionErrorStat(Sys) the sigma of its
1539 // statistical (systematic) error
1542 fractionErrorStat = 999.;
1543 fractionErrorSys = 999.;
1545 if (species > AliPID::kProton || species < AliPID::kElectron) {
1546 AliError(Form("Only fractions for species index %d to %d availabe, but not for the requested one: %d", 0, AliPID::kProton, species));
1550 if (!fFractionHists[species]) {
1551 AliError(Form("Histo with particle fractions for species %d not loaded!", species));
1556 Int_t jetPtBin = FindBinWithinRange(fFractionHists[species]->GetYaxis(), jetPt);
1557 Int_t centBin = FindBinWithinRange(fFractionHists[species]->GetZaxis(), centralityPercentile);
1559 // The following interpolation takes the bin content of the first/last available bin,
1560 // if requested point lies beyond bin center of first/last bin.
1561 // The interpolation is only done for the x-axis (track pT), i.e. jetPtBin and centBin are fix,
1562 // because the analysis will anyhow be run in bins of jetPt and centrality and
1563 // it is not desired to correlate different jetPt bins via interpolation.
1565 // The same procedure is used for the error of the fraction
1566 TAxis* xAxis = fFractionHists[species]->GetXaxis();
1568 // No interpolation to values beyond the centers of the first/last bins (we don't know exactly where the spectra start or stop,
1569 // thus, search for the first and last bin above 0.0 to constrain the range
1570 Int_t firstBin = TMath::Max(1, FindFirstBinAboveIn3dSubset(fFractionHists[species], 0.0, jetPtBin, centBin));
1571 Int_t lastBin = TMath::Min(fFractionHists[species]->GetNbinsX(),
1572 FindLastBinAboveIn3dSubset(fFractionHists[species], 0.0, jetPtBin, centBin));
1574 if (trackPt <= xAxis->GetBinCenter(firstBin)) {
1575 fraction = fFractionHists[species]->GetBinContent(firstBin, jetPtBin, centBin);
1576 fractionErrorStat = fFractionHists[species]->GetBinError(firstBin, jetPtBin, centBin);
1577 fractionErrorSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(firstBin, jetPtBin, centBin) : 0.;
1579 else if (trackPt >= xAxis->GetBinCenter(lastBin)) {
1580 fraction = fFractionHists[species]->GetBinContent(lastBin, jetPtBin, centBin);
1581 fractionErrorStat = fFractionHists[species]->GetBinError(lastBin, jetPtBin, centBin);
1582 fractionErrorSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(lastBin, jetPtBin, centBin) : 0.;
1585 Double_t x0 = 0., x1 = 0., y0 = 0., y1 = 0.;
1586 Double_t y0errStat = 0., y1errStat = 0., y0errSys = 0., y1errSys = 0.;
1587 Int_t trackPtBin = xAxis->FindFixBin(trackPt);
1589 // Linear interpolation between nearest neighbours in trackPt
1590 if (trackPt <= xAxis->GetBinCenter(trackPtBin)) {
1591 y0 = fFractionHists[species]->GetBinContent(trackPtBin - 1, jetPtBin, centBin);
1592 y0errStat = fFractionHists[species]->GetBinError(trackPtBin - 1, jetPtBin, centBin);
1593 y0errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin - 1, jetPtBin, centBin)
1595 x0 = xAxis->GetBinCenter(trackPtBin - 1);
1596 y1 = fFractionHists[species]->GetBinContent(trackPtBin, jetPtBin, centBin);
1597 y1errStat = fFractionHists[species]->GetBinError(trackPtBin, jetPtBin, centBin);
1598 y1errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin, jetPtBin, centBin)
1600 x1 = xAxis->GetBinCenter(trackPtBin);
1603 y0 = fFractionHists[species]->GetBinContent(trackPtBin, jetPtBin, centBin);
1604 y0errStat = fFractionHists[species]->GetBinError(trackPtBin, jetPtBin, centBin);
1605 y0errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin, jetPtBin, centBin)
1607 x0 = xAxis->GetBinCenter(trackPtBin);
1608 y1 = fFractionHists[species]->GetBinContent(trackPtBin + 1, jetPtBin, centBin);
1609 y1errStat = fFractionHists[species]->GetBinError(trackPtBin + 1, jetPtBin, centBin);
1610 y1errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin + 1, jetPtBin, centBin)
1612 x1 = xAxis->GetBinCenter(trackPtBin + 1);
1615 // Per construction: x0 < trackPt < x1
1616 fraction = y0 + (trackPt - x0) * ((y1 - y0) / (x1 - x0));
1617 fractionErrorStat = y0errStat + (trackPt - x0) * ((y1errStat - y0errStat) / (x1 - x0));
1618 fractionErrorSys = fFractionSysErrorHists[species] ? (y0errSys + (trackPt - x0) * ((y1errSys - y0errSys) / (x1 - x0))) : 0.;
1625 //_____________________________________________________________________________
1626 Bool_t AliAnalysisTaskPID::GetParticleFractions(Double_t trackPt, Double_t jetPt, Double_t centralityPercentile,
1627 Double_t* prob, Int_t smearSpeciesByError,
1628 Int_t takeIntoAccountSpeciesSysError, Bool_t uniformSystematicError) const
1630 // Fills the particle fractions for the given trackPt, jetPt and centrality into "prob".
1631 // Use jetPt = -1 for inclusive spectra and centralityPercentile = -1 for pp.
1632 // If smearSpeciesByError is >= 0 && < AliPID::kSPECIES, the returned fractions will be a random number distributed
1633 // with a gauss with mean being the corresponding particle fraction and sigma it's error for the considered species
1634 // "smearSpeciesByError".
1635 // Note that in this case the fractions for all species will NOT sum up to 1!
1636 // Thus, all other species fractions will be re-scaled weighted with their corresponding statistical error.
1637 // A similar procedure is used for "takeIntoAccountSpeciesSysError": The systematic error of the corresponding species
1638 // is used to generate a random number with uniform distribution in [mean - sysError, mean + sysError] for the new mean
1639 // (in cace of uniformSystematicError = kTRUE, otherwise it will be a gaus(mean, sysError)),
1640 // then the other species will be re-scaled according to their systematic errors.
1641 // First, the systematic error uncertainty procedure will be performed (that is including re-scaling), then the statistical
1642 // uncertainty procedure.
1643 // On success, kTRUE is returned.
1645 if (!prob || smearSpeciesByError >= AliPID::kSPECIES || takeIntoAccountSpeciesSysError >= AliPID::kSPECIES)
1648 Double_t probTemp[AliPID::kSPECIES];
1649 Double_t probErrorStat[AliPID::kSPECIES];
1650 Double_t probErrorSys[AliPID::kSPECIES];
1652 Bool_t success = kTRUE;
1653 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kElectron,
1654 probTemp[AliPID::kElectron], probErrorStat[AliPID::kElectron],
1655 probErrorSys[AliPID::kElectron]);
1656 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kMuon,
1657 probTemp[AliPID::kMuon], probErrorStat[AliPID::kMuon], probErrorSys[AliPID::kMuon]);
1658 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kPion,
1659 probTemp[AliPID::kPion], probErrorStat[AliPID::kPion], probErrorSys[AliPID::kPion]);
1660 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kKaon,
1661 probTemp[AliPID::kKaon], probErrorStat[AliPID::kKaon], probErrorSys[AliPID::kKaon]);
1662 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kProton,
1663 probTemp[AliPID::kProton], probErrorStat[AliPID::kProton], probErrorSys[AliPID::kProton]);
1668 // If desired, take into account the systematic error of the corresponding species and re-generate probTemp accordingly
1669 if (takeIntoAccountSpeciesSysError >= 0) {
1670 // Generate random fraction of the considered species "smearSpeciesByError" according to mean and sigma
1671 Double_t generatedFraction = uniformSystematicError
1672 ? fRandom->Rndm() * 2. * probErrorSys[takeIntoAccountSpeciesSysError]
1673 - probErrorSys[takeIntoAccountSpeciesSysError]
1674 + probTemp[takeIntoAccountSpeciesSysError]
1675 : fRandom->Gaus(probTemp[takeIntoAccountSpeciesSysError],
1676 probErrorSys[takeIntoAccountSpeciesSysError]);
1678 // Catch cases with invalid fraction (can happen for large errors), i.e. fraction < 0 or > 1
1679 if (generatedFraction < 0.)
1680 generatedFraction = 0.;
1681 else if (generatedFraction > 1.)
1682 generatedFraction = 1.;
1684 // Calculate difference from original fraction (original fractions sum up to 1!)
1685 Double_t deltaFraction = generatedFraction - probTemp[takeIntoAccountSpeciesSysError];
1687 // Fractions must (including errors) lie inside [0,1] -> Adapt weights accordingly by setting the errors
1688 if (deltaFraction > 0) {
1689 // Some part will be SUBTRACTED from the other fractions
1690 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1691 if (probTemp[i] - probErrorSys[i] < 0)
1692 probErrorSys[i] = probTemp[i];
1696 // Some part will be ADDED to the other fractions
1697 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1698 if (probTemp[i] + probErrorSys[i] > 1)
1699 probErrorSys[i] = 1. - probTemp[i];
1703 // Compute summed weight of all fractions except for the considered one
1704 Double_t summedWeight = 0.;
1705 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1706 if (i != takeIntoAccountSpeciesSysError)
1707 summedWeight += probErrorSys[i];
1710 // Compute the weight for the other species
1712 if (summedWeight <= 1e-13) {
1713 // If this happens for some reason (it should not!), just assume flat weight
1714 printf("Error: summedWeight (sys error) ~ 0 for trackPt %f, jetPt %f, centralityPercentile %f. Setting flat weight!\n",
1715 trackPt, jetPt, centralityPercentile);
1718 Double_t weight[AliPID::kSPECIES];
1719 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1720 if (i != takeIntoAccountSpeciesSysError) {
1721 if (summedWeight > 1e-13)
1722 weight[i] = probErrorSys[i] / summedWeight;
1724 weight[i] = probErrorSys[i] / (AliPID::kSPECIES - 1);
1728 // For the final generated fractions, set the generated value for the considered species
1729 // and the generated value minus delta times statistical weight
1730 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1731 if (i != takeIntoAccountSpeciesSysError)
1732 probTemp[i] = probTemp[i] - weight[i] * deltaFraction;
1734 probTemp[i] = generatedFraction;
1738 // Using the values of probTemp (either the original ones or those after taking into account the systematic error),
1739 // calculate the final fractions - if the statistical error is to be taken into account, smear the corresponding
1740 // fraction. If not, just write probTemp to the final result array.
1741 if (smearSpeciesByError >= 0) {
1742 // Generate random fraction of the considered species "smearSpeciesByError" according to mean and sigma
1743 Double_t generatedFraction = fRandom->Gaus(probTemp[smearSpeciesByError], probErrorStat[smearSpeciesByError]);
1745 // Catch cases with invalid fraction (can happen for large errors), i.e. fraction < 0 or > 1
1746 if (generatedFraction < 0.)
1747 generatedFraction = 0.;
1748 else if (generatedFraction > 1.)
1749 generatedFraction = 1.;
1751 // Calculate difference from original fraction (original fractions sum up to 1!)
1752 Double_t deltaFraction = generatedFraction - probTemp[smearSpeciesByError];
1754 // Fractions must (including errors) lie inside [0,1] -> Adapt weights accordingly by setting the errors
1755 if (deltaFraction > 0) {
1756 // Some part will be SUBTRACTED from the other fractions
1757 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1758 if (probTemp[i] - probErrorStat[i] < 0)
1759 probErrorStat[i] = probTemp[i];
1763 // Some part will be ADDED to the other fractions
1764 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1765 if (probTemp[i] + probErrorStat[i] > 1)
1766 probErrorStat[i] = 1. - probTemp[i];
1770 // Compute summed weight of all fractions except for the considered one
1771 Double_t summedWeight = 0.;
1772 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1773 if (i != smearSpeciesByError)
1774 summedWeight += probErrorStat[i];
1777 // Compute the weight for the other species
1779 if (summedWeight <= 1e-13) {
1780 // If this happens for some reason (it should not!), just assume flat weight
1781 printf("Error: summedWeight (stat error) ~ 0 for trackPt %f, jetPt %f, centralityPercentile %f. Setting flat weight!\n",
1782 trackPt, jetPt, centralityPercentile);
1785 Double_t weight[AliPID::kSPECIES];
1786 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1787 if (i != smearSpeciesByError) {
1788 if (summedWeight > 1e-13)
1789 weight[i] = probErrorStat[i] / summedWeight;
1791 weight[i] = probErrorStat[i] / (AliPID::kSPECIES - 1);
1795 // For the final generated fractions, set the generated value for the considered species
1796 // and the generated value minus delta times statistical weight
1797 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1798 if (i != smearSpeciesByError)
1799 prob[i] = probTemp[i] - weight[i] * deltaFraction;
1801 prob[i] = generatedFraction;
1805 // Just take the generated values
1806 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
1807 prob[i] = probTemp[i];
1811 // Should already be normalised, but make sure that it really is:
1812 Double_t probSum = 0.;
1813 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1820 if (TMath::Abs(probSum - 1.0) > 1e-4) {
1821 printf("Warning: Re-normalising sum of fractions: Sum is %e\n", probSum);
1822 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1831 //_____________________________________________________________________________
1832 const TH3D* AliAnalysisTaskPID::GetParticleFractionHisto(Int_t species, Bool_t sysError) const
1834 if (species < AliPID::kElectron || species > AliPID::kProton)
1837 return sysError ? fFractionSysErrorHists[species] : fFractionHists[species];
1841 //_____________________________________________________________________________
1842 Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(Int_t motherPDG, Double_t motherGenPt)
1844 // Strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
1845 // -> Based on function in PWGJE/AliAnalysisTaskFragmentationFunction, which uses
1846 // the following data from CMS pp @ 7 TeV inclusive (JHEP 05 (2011) 064)
1850 const Int_t absMotherPDG = TMath::Abs(motherPDG);
1852 if (absMotherPDG == 310 || absMotherPDG == 321) { // K0s / K+ / K-
1853 if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.768049;
1854 else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.732933;
1855 else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.650298;
1856 else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.571332;
1857 else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.518734;
1858 else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.492543;
1859 else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.482704;
1860 else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.488056;
1861 else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.488861;
1862 else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.492862;
1863 else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.504332;
1864 else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.501858;
1865 else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.512970;
1866 else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.524131;
1867 else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.539130;
1868 else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.554101;
1869 else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.560348;
1870 else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.568869;
1871 else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.583310;
1872 else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.604818;
1873 else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.632630;
1874 else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.710070;
1875 else if(6.00 <= motherGenPt && motherGenPt < 8.00) fac = 0.736365;
1876 else if(8.00 <= motherGenPt && motherGenPt < 10.00) fac = 0.835865;
1879 if (absMotherPDG == 3122) { // Lambda
1880 //if (absMotherPDG == 3122 || absMotherPDG == 3112 || absMotherPDG == 3222) { // Lambda / Sigma- / Sigma+
1881 if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.645162;
1882 else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.627431;
1883 else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.457136;
1884 else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.384369;
1885 else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.330597;
1886 else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.309571;
1887 else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.293620;
1888 else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.283709;
1889 else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.282047;
1890 else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.277261;
1891 else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.275772;
1892 else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.280726;
1893 else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.288540;
1894 else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.288315;
1895 else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.296619;
1896 else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.302993;
1897 else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.338121;
1898 else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.349800;
1899 else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.356802;
1900 else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.391202;
1901 else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.422573;
1902 else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.573815;
1903 else if(6.00 <= motherGenPt && motherGenPt < 8.00) fac = 0.786984;
1904 else if(8.00 <= motherGenPt && motherGenPt < 10.00) fac = 1.020021;
1907 if (absMotherPDG == 3312 || absMotherPDG == 3322) { // xi
1908 if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.666620;
1909 else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.575908;
1910 else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.433198;
1911 else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.340901;
1912 else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.290896;
1913 else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.236074;
1914 else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.218681;
1915 else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.207763;
1916 else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.222848;
1917 else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.208806;
1918 else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.197275;
1919 else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.183645;
1920 else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.188788;
1921 else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.188282;
1922 else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.207442;
1923 else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.240388;
1924 else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.241916;
1925 else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.208276;
1926 else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.234550;
1927 else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.251689;
1928 else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.310204;
1929 else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.343492;
1932 const Double_t weight = 1. / fac;
1938 //_____________________________________________________________________________
1939 Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(AliMCEvent* mcEvent, AliMCParticle* daughter)
1941 // Strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
1942 // -> Based on function in PWGJE/AliAnalysisTaskFragmentationFunction
1947 AliMCParticle* currentMother = daughter;
1948 AliMCParticle* currentDaughter = daughter;
1951 // find first primary mother K0s, Lambda or Xi
1953 Int_t daughterPDG = currentDaughter->PdgCode();
1955 Int_t motherLabel = currentDaughter->GetMother();
1956 if(motherLabel >= mcEvent->GetNumberOfTracks()){ // protection
1957 currentMother = currentDaughter;
1961 currentMother = (AliMCParticle*)mcEvent->GetTrack(motherLabel);
1963 if (!currentMother) {
1964 currentMother = currentDaughter;
1968 Int_t motherPDG = currentMother->PdgCode();
1970 // phys. primary found ?
1971 if (mcEvent->IsPhysicalPrimary(motherLabel))
1974 if (TMath::Abs(daughterPDG) == 321) {
1975 // K+/K- e.g. from phi (ref data not feeddown corrected)
1976 currentMother = currentDaughter;
1979 if (TMath::Abs(motherPDG) == 310) {
1980 // K0s e.g. from phi (ref data not feeddown corrected)
1983 if (TMath::Abs(motherPDG) == 3212 && TMath::Abs(daughterPDG) == 3122) {
1984 // Mother Sigma0, daughter Lambda (this case not included in feeddown corr.)
1985 currentMother = currentDaughter;
1989 currentDaughter = currentMother;
1993 Int_t motherPDG = currentMother->PdgCode();
1994 Double_t motherGenPt = currentMother->Pt();
1996 return GetMCStrangenessFactorCMS(motherPDG, motherGenPt);
2000 // _________________________________________________________________________________
2001 AliAnalysisTaskPID::TOFpidInfo AliAnalysisTaskPID::GetTOFType(const AliVTrack* track, Int_t tofMode) const
2003 // Get the (locally defined) particle type judged by TOF
2005 if (!fPIDResponse) {
2006 Printf("ERROR: fPIDResponse not available -> Cannot determine TOF type!");
2010 /*TODO still needs some further thinking how to implement it....
2011 // TOF PID status (kTOFout, kTIME) automatically checked by return value of ComputeTPCProbability;
2012 // also, probability array will be set there (no need to initialise here)
2013 Double_t p[AliPID::kSPECIES];
2014 const AliPIDResponse::EDetPidStatus tofStatus = fPIDResponse->ComputeTPCProbability(track, AliPID::kSPECIES, p);
2015 if (tofStatus != AliPIDResponse::kDetPidOk)
2018 // Do not consider muons
2019 p[AliPID::kMuon] = 0.;
2021 // Probabilities are not normalised yet
2023 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2029 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2032 Double_t probThreshold = -999.;
2034 // If there is only one distribution, the threshold corresponds to...
2038 else if (tofMode == 1) { // default
2039 probThreshold = 0.9973; // a 3-sigma inclusion cut
2041 else if (tofMode == 2) {
2046 Printf("ERROR: Bad TOF mode: %d!", tofMode);
2052 ///* OLD: cut with n-sigmas (ATTENTION: Does not take into account properly the TOF tail!)
2053 // Check kTOFout, kTIME, mismatch
2054 const AliPIDResponse::EDetPidStatus tofStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, track);
2055 if (tofStatus != AliPIDResponse::kDetPidOk)
2058 Double_t nsigma[kNumTOFspecies + 1] = { -999., -999., -999., -999. };
2059 nsigma[kTOFpion] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion);
2060 nsigma[kTOFkaon] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon);
2061 nsigma[kTOFproton] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton);
2063 Double_t inclusion = -999;
2064 Double_t exclusion = -999;
2070 else if (tofMode == 1) { // default
2074 else if (tofMode == 2) {
2079 Printf("ERROR: Bad TOF mode: %d!", tofMode);
2083 // Do not cut on nSigma electron because this would also remove pions in a large pT range.
2084 // The electron contamination of the pion sample can be treated by TPC dEdx fits afterwards.
2085 if (TMath::Abs(nsigma[kTOFpion]) < inclusion && TMath::Abs(nsigma[kTOFkaon]) > exclusion && TMath::Abs(nsigma[kTOFproton]) > exclusion)
2087 if (TMath::Abs(nsigma[kTOFpion]) > exclusion && TMath::Abs(nsigma[kTOFkaon]) < inclusion && TMath::Abs(nsigma[kTOFproton]) > exclusion)
2089 if (TMath::Abs(nsigma[kTOFpion]) > exclusion && TMath::Abs(nsigma[kTOFkaon]) > exclusion && TMath::Abs(nsigma[kTOFproton]) < inclusion)
2093 // There are no TOF electrons selected because the purity is rather bad, even for small momenta
2094 // (also a small mismatch probability significantly affects electrons because their fraction
2095 // is small). There is no need for TOF electrons anyway, since the dEdx distribution of kaons and
2096 // protons in a given pT bin (also at the dEdx crossings) is very different from that of electrons.
2097 // This is due to the steeply falling dEdx of p and K with momentum, whereas the electron dEdx stays
2099 // As a result, the TPC fit yields a more accurate electron fraction than the TOF selection can do.
2105 // _________________________________________________________________________________
2106 Bool_t AliAnalysisTaskPID::IsSecondaryWithStrangeMotherMC(AliMCEvent* mcEvent, Int_t partLabel)
2108 // Check whether particle is a secondary with strange mother, i.e. returns kTRUE if a strange mother is found
2109 // and the particle is NOT a physical primary. In all other cases kFALSE is returned
2111 if (!mcEvent || partLabel < 0)
2114 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(partLabel);
2119 if (mcEvent->IsPhysicalPrimary(partLabel))
2122 Int_t iMother = part->GetMother();
2127 AliMCParticle* partM = (AliMCParticle*)mcEvent->GetTrack(iMother);
2131 Int_t codeM = TMath::Abs(partM->PdgCode());
2132 Int_t mfl = Int_t(codeM / TMath::Power(10, Int_t(TMath::Log10(codeM))));
2133 if (mfl == 3 && codeM != 3) // codeM = 3 is for s quark
2140 //_____________________________________________________________________________
2141 Bool_t AliAnalysisTaskPID::SetParticleFractionHisto(const TH3D* hist, Int_t species, Bool_t sysError)
2143 // Store a clone of hist (containing the particle fractions of the corresponding species with statistical error (sysError = kFALSE)
2144 // or systematic error (sysError = kTRUE), respectively), internally
2146 if (species < AliPID::kElectron || species > AliPID::kProton) {
2147 AliError(Form("Only fractions for species index %d to %d can be set, but not for the requested one: %d", 0,
2148 AliPID::kProton, species));
2153 delete fFractionSysErrorHists[species];
2155 fFractionSysErrorHists[species] = new TH3D(*hist);
2158 delete fFractionHists[species];
2160 fFractionHists[species] = new TH3D(*hist);
2167 //_____________________________________________________________________________
2168 Bool_t AliAnalysisTaskPID::SetParticleFractionHistosFromFile(const TString filePathName, Bool_t sysError)
2170 // Loads particle fractions for all species from the desired file and returns kTRUE on success.
2171 // The maps are assumed to be of Type TH3D, to sit in the main directory and to have names
2172 // Form("hFraction_%e", AliPID::ParticleName(i)) for sysError = kFALSE and
2173 // Form("hFractionSysError_%e", AliPID::ParticleName(i)) for sysError = kTRUE.
2175 TFile* f = TFile::Open(filePathName.Data());
2177 std::cout << "Failed to open file with particle fractions \"" << filePathName.Data() << "\"!" << std::endl;
2182 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
2183 TString histName = Form("hFraction%s_%s", sysError ? "SysError" : "", AliPID::ParticleName(i));
2184 hist = dynamic_cast<TH3D*>(f->Get(histName.Data()));
2186 std::cout << "Failed to load particle fractions for " << histName.Data() << "!";
2187 std::cout << std::endl << "Cleaning up particle fraction histos!" << std::endl;
2188 CleanupParticleFractionHistos();
2192 if (!SetParticleFractionHisto(hist, i, sysError)) {
2193 std::cout << "Failed to load particle fractions for " << histName.Data() << "!";
2194 std::cout << std::endl << "Cleaning up particle fraction histos!" << std::endl;
2195 CleanupParticleFractionHistos();
2207 //_____________________________________________________________________________
2208 Double_t AliAnalysisTaskPID::GetMaxEtaVariation(Double_t dEdxSplines)
2210 // Returns the maximum eta variation (i.e. deviation of eta correction factor from unity) for the
2211 // given (spline) dEdx
2213 if (dEdxSplines < 1. || !fhMaxEtaVariation) {
2214 Printf("Error GetMaxEtaVariation: No map or invalid dEdxSplines (%f)!", dEdxSplines);
2218 Int_t bin = fhMaxEtaVariation->GetXaxis()->FindFixBin(1. / dEdxSplines);
2222 if (bin > fhMaxEtaVariation->GetXaxis()->GetNbins())
2223 bin = fhMaxEtaVariation->GetXaxis()->GetNbins();
2225 return fhMaxEtaVariation->GetBinContent(bin);
2229 //_____________________________________________________________________________
2230 Int_t AliAnalysisTaskPID::GetRandomParticleTypeAccordingToParticleFractions(Double_t trackPt, Double_t jetPt,
2231 Double_t centralityPercentile,
2232 Bool_t smearByError,
2233 Bool_t takeIntoAccountSysError) const
2235 // Uses the stored histograms with the particle fractions to generate a random particle type according to these fractions.
2236 // In case of problems (e.g. histo missing), AliPID::kUnknown is returned.
2237 // If smearByError is kTRUE, the used fractions will be random numbers distributed with a gauss with mean
2238 // being the corresponding particle fraction and sigma it's error.
2239 // Note that in this case only the fraction of a random species is varied in this way. The other fractions
2240 // will be re-normalised according their statistical errors.
2241 // The same holds for the systematic error of species "takeIntoAccountSpeciesSysError", but the random number will be
2242 // uniformly distributed within [mean - sys, mean + sys] and the re-normalisation will be weighted with the systematic errors.
2243 // Note that the fractions will be calculated first with only the systematic error taken into account (if desired), including
2244 // re-normalisation. Then, the resulting fractions will be used to calculate the final fractions - either with statistical error
2245 // or without. The species, for which the error will be used for smearing, is the same for sys and stat error.
2247 Double_t prob[AliPID::kSPECIES];
2248 Int_t randomSpecies = (smearByError || takeIntoAccountSysError) ? (Int_t)(fRandom->Rndm() * AliPID::kSPECIES) : -1;
2249 Bool_t success = GetParticleFractions(trackPt, jetPt, centralityPercentile, prob, randomSpecies, randomSpecies);
2252 return AliPID::kUnknown;
2254 Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
2256 if (rnd <= prob[AliPID::kPion])
2257 return AliPID::kPion;
2258 else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon])
2259 return AliPID::kKaon;
2260 else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon] + prob[AliPID::kProton])
2261 return AliPID::kProton;
2262 else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon] + prob[AliPID::kProton] + prob[AliPID::kElectron])
2263 return AliPID::kElectron;
2265 return AliPID::kMuon; //else it must be a muon (only species left)
2269 //_____________________________________________________________________________
2270 AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::GenerateDetectorResponse(AliAnalysisTaskPID::ErrorCode errCode,
2271 Double_t mean, Double_t sigma,
2272 Double_t* responses, Int_t nResponses,
2275 // Generate detector response. If a previous generation was not successful or there is something wrong with this signal generation,
2276 // the function will return kFALSE
2280 // Reset response array
2281 for (Int_t i = 0; i < nResponses; i++)
2282 responses[i] = -999;
2284 if (errCode == kError)
2287 ErrorCode ownErrCode = kNoErrors;
2289 if (fUseConvolutedGaus && !usePureGaus) {
2290 // In case of convoluted gauss, calculate the probability density only once to save a lot of time!
2292 TH1* hProbDensity = 0x0;
2293 ownErrCode = SetParamsForConvolutedGaus(mean, sigma);
2294 if (ownErrCode == kError)
2297 hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
2299 for (Int_t i = 0; i < nResponses; i++) {
2300 responses[i] = hProbDensity->GetRandom();
2301 //responses[i] fConvolutedGausDeltaPrime->GetRandom(); // MUCH slower than using the binned version via the histogram
2305 for (Int_t i = 0; i < nResponses; i++) {
2306 responses[i] = fRandom->Gaus(mean, sigma);
2310 // If forwarded error code was a warning (error case has been handled before), return a warning
2311 if (errCode == kWarning)
2314 return ownErrCode; // Forward success/warning
2318 //_____________________________________________________________________________
2319 void AliAnalysisTaskPID::PrintSettings(Bool_t printSystematicsSettings) const
2321 // Print current settings.
2323 printf("\n\nSettings for task %s:\n", GetName());
2324 printf("Is pPb/Pbp: %d -> %s\n", GetIsPbpOrpPb(), GetIsPbpOrpPb() ? "Adapting vertex cuts" : "Using standard vertex cuts");
2325 printf("Track cuts: %s\n", fTrackFilter ? fTrackFilter->GetTitle() : "-");
2326 printf("Eta cut: %.2f <= |eta| <= %.2f\n", GetEtaAbsCutLow(), GetEtaAbsCutUp());
2327 printf("Phi' cut: %d\n", GetUsePhiCut());
2328 printf("TPCCutMIGeo: %d\n", GetUseTPCCutMIGeo());
2329 if (GetUseTPCCutMIGeo()) {
2330 printf("GetCutGeo: %f\n", GetCutGeo());
2331 printf("GetCutNcr: %f\n", GetCutNcr());
2332 printf("GetCutNcl: %f\n", GetCutNcl());
2334 printf("TPCnclCut: %d\n", GetUseTPCnclCut());
2335 if (GetUseTPCnclCut()) {
2336 printf("GetCutPureNcl: %d\n", GetCutPureNcl());
2341 printf("Centrality estimator: %s\n", GetCentralityEstimator().Data());
2345 printf("Pile up rejection type: %d\n", (Int_t)fPileUpRejectionType);
2349 printf("Use MC-ID for signal generation: %d\n", GetUseMCidForGeneration());
2350 printf("Use ITS: %d\n", GetUseITS());
2351 printf("Use TOF: %d\n", GetUseTOF());
2352 printf("Use priors: %d\n", GetUsePriors());
2353 printf("Use TPC default priors: %d\n", GetUseTPCDefaultPriors());
2354 printf("Use convoluted Gauss: %d\n", GetUseConvolutedGaus());
2355 printf("Accuracy of non-Gaussian tail: %e\n", GetAccuracyNonGaussianTail());
2356 printf("Take into account muons: %d\n", GetTakeIntoAccountMuons());
2357 printf("TOF mode: %d\n", GetTOFmode());
2358 printf("\nParams for transition from gauss to asymmetric shape:\n");
2359 printf("[0]: %e\n", GetConvolutedGaussTransitionPar(0));
2360 printf("[1]: %e\n", GetConvolutedGaussTransitionPar(1));
2361 printf("[2]: %e\n", GetConvolutedGaussTransitionPar(2));
2365 printf("Do PID: %d\n", fDoPID);
2366 printf("Do Efficiency: %d\n", fDoEfficiency);
2367 printf("Do PtResolution: %d\n", fDoPtResolution);
2368 printf("Do dEdxCheck: %d\n", fDoDeDxCheck);
2369 printf("Do binZeroStudy: %d\n", fDoBinZeroStudy);
2373 printf("Input from other task: %d\n", GetInputFromOtherTask());
2374 printf("Store additional jet information: %d\n", GetStoreAdditionalJetInformation());
2375 printf("Store centrality percentile: %d", GetStoreCentralityPercentile());
2377 if (printSystematicsSettings)
2378 PrintSystematicsSettings();
2384 //_____________________________________________________________________________
2385 void AliAnalysisTaskPID::PrintSystematicsSettings() const
2387 // Print current settings for systematic studies.
2389 printf("\n\nSettings for systematics for task %s:\n", GetName());
2390 printf("SplinesThr:\t%f\n", GetSystematicScalingSplinesThreshold());
2391 printf("SplinesBelowThr:\t%f\n", GetSystematicScalingSplinesBelowThreshold());
2392 printf("SplinesAboveThr:\t%f\n", GetSystematicScalingSplinesAboveThreshold());
2393 printf("EtaCorrMomThr:\t%f\n", GetSystematicScalingEtaCorrectionMomentumThr());
2394 printf("EtaCorrLowP:\t%f\n", GetSystematicScalingEtaCorrectionLowMomenta());
2395 printf("EtaCorrHighP:\t%f\n", GetSystematicScalingEtaCorrectionHighMomenta());
2396 printf("SigmaParaThr:\t%f\n", GetSystematicScalingEtaSigmaParaThreshold());
2397 printf("SigmaParaBelowThr:\t%f\n", GetSystematicScalingEtaSigmaParaBelowThreshold());
2398 printf("SigmaParaAboveThr:\t%f\n", GetSystematicScalingEtaSigmaParaAboveThreshold());
2399 printf("MultCorr:\t%f\n", GetSystematicScalingMultCorrection());
2400 printf("TOF mode: %d\n", GetTOFmode());
2406 //_____________________________________________________________________________
2407 Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePDGcode, Double_t centralityPercentile,
2410 // Process the track (generate expected response, fill histos, etc.).
2411 // particlePDGcode == 0 means data. Otherwise, the corresponding MC ID will be assumed.
2413 //Printf("Debug: Task %s is starting to process track: dEdx %f, pTPC %f, eta %f, ncl %d\n", GetName(), track->GetTPCsignal(), track->GetTPCmomentum(),
2414 // track->Eta(), track->GetTPCsignalN());
2417 printf("File: %s, Line: %d: ProcessTrack\n", (char*)__FILE__, __LINE__);
2419 if (!fDoPID && !fDoDeDxCheck && !fDoPtResolution)
2423 printf("File: %s, Line: %d: ProcessTrack -> Processing started\n", (char*)__FILE__, __LINE__);
2425 const Bool_t isMC = (particlePDGcode == 0) ? kFALSE : kTRUE;
2430 if (TMath::Abs(particlePDGcode) == 211) {//Pion
2433 else if (TMath::Abs(particlePDGcode) == 321) {//Kaon
2436 else if (TMath::Abs(particlePDGcode) == 2212) {//Proton
2439 else if (TMath::Abs(particlePDGcode) == 11) {//Electron
2442 else if (TMath::Abs(particlePDGcode) == 13) {//Muon
2445 else // In MC-ID case, set to underflow bin such that the response from this track is only used for unidentified signal generation
2446 // or signal generation with PID response and the track is still there (as in data) - e.g. relevant w.r.t. deuterons.
2447 // This is important to be as much as possible consistent with data. And the tracks can still be removed by disabling the
2448 // underflow bin for the projections
2453 //Double_t p = track->GetP();
2454 const Double_t pTPC = track->GetTPCmomentum();
2455 Double_t pT = track->Pt();
2457 Double_t z = -1, xi = -1;
2458 GetJetTrackObservables(pT, jetPt, z, xi);
2461 Double_t trackCharge = track->Charge();
2464 const Bool_t tuneOnDataTPC = fPIDResponse->IsTunedOnData() &&
2465 ((fPIDResponse->GetTunedOnDataMask() & AliPIDResponse::kDetTPC) == AliPIDResponse::kDetTPC);
2466 Double_t dEdxTPC = tuneOnDataTPC ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
2470 Printf("Skipping track with strange dEdx value: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), pTPC,
2471 track->Eta(), track->GetTPCsignalN());
2475 // Completely remove tracks from light nuclei defined by el and pr <dEdx> as function of pT.
2476 // A very loose cut to be sure not to throw away electrons and/or protons
2477 Double_t nSigmaPr = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton);
2478 Double_t nSigmaEl = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
2480 if ((pTPC >= 0.3 && (nSigmaPr > 10 && nSigmaEl > 10)) ||
2481 (pTPC < 0.3 && (nSigmaPr > 15 && nSigmaEl > 15))) {
2483 Printf("Skipping track which seems to be a light nucleus: dEdx %f, pTPC %f, pT %f, eta %f, ncl %d, nSigmaPr %f, nSigmaEl %f\n",
2484 track->GetTPCsignal(), pTPC, pT, track->Eta(), track->GetTPCsignalN(), nSigmaPr, nSigmaEl);
2489 Double_t dEdxEl, dEdxKa, dEdxPi, dEdxMu, dEdxPr;
2490 Double_t sigmaEl, sigmaKa, sigmaPi, sigmaMu, sigmaPr;
2492 if (fDoAnySystematicStudiesOnTheExpectedSignal) {
2493 // Get the uncorrected signal first and the corresponding correction factors.
2494 // Then modify the correction factors and properly recalculate the corrected dEdx
2496 // Get pure spline values for dEdx_expected, without any correction
2497 dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2498 dEdxKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2499 dEdxPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2500 dEdxMu = !fTakeIntoAccountMuons ? -1 :
2501 fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2502 dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2504 // Scale splines, if desired
2505 if ((TMath::Abs(fSystematicScalingSplinesBelowThreshold - 1.0) > fgkEpsilon) ||
2506 (TMath::Abs(fSystematicScalingSplinesAboveThreshold - 1.0) > fgkEpsilon)) {
2508 // Tune turn-on of correction for pions (not so relevant for the other species, since at very large momenta)
2510 Double_t scaleFactor = 1.;
2511 Double_t usedSystematicScalingSplines = 1.;
2513 bg = pTPC / AliPID::ParticleMass(AliPID::kElectron);
2514 scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2515 usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2516 + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2517 dEdxEl *= usedSystematicScalingSplines;
2519 bg = pTPC / AliPID::ParticleMass(AliPID::kKaon);
2520 scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2521 usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2522 + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2523 dEdxKa *= usedSystematicScalingSplines;
2525 bg = pTPC / AliPID::ParticleMass(AliPID::kPion);
2526 scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2527 usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2528 + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2529 dEdxPi *= usedSystematicScalingSplines;
2531 if (fTakeIntoAccountMuons) {
2532 bg = pTPC / AliPID::ParticleMass(AliPID::kMuon);
2533 scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2534 usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2535 + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2536 dEdxMu *= usedSystematicScalingSplines;
2540 bg = pTPC / AliPID::ParticleMass(AliPID::kProton);
2541 scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2542 usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2543 + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2544 dEdxPr *= usedSystematicScalingSplines;
2547 // Get the eta correction factors for the (modified) expected dEdx
2548 Double_t etaCorrEl = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxEl) : 1.;
2549 Double_t etaCorrKa = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxKa) : 1.;
2550 Double_t etaCorrPi = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPi) : 1.;
2551 Double_t etaCorrMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCEtaCorrection() ?
2552 fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxMu) : 1.;
2553 Double_t etaCorrPr = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPr) : 1.;
2555 // Scale eta correction factors, if desired (and eta correction maps are to be used, otherwise it is not possible!)
2556 if (fPIDResponse->UseTPCEtaCorrection() &&
2557 (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon ||
2558 TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon)) {
2559 // Since we do not want to scale the splines with this, but only the eta variation, only scale the deviation of the correction factor!
2560 // 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!
2563 // Due to additional azimuthal effects, there is an additional eta dependence for low momenta which is not corrected successfully so far.
2564 // One can assign a different (higher) systematic scale factor for this low-p region and a threshold which separates low- and high-p.
2565 // 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
2566 Double_t usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionHighMomenta;
2568 if (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - fSystematicScalingEtaCorrectionLowMomenta) > fgkEpsilon) {
2569 const Double_t fractionHighMomentumScaleFactor = 0.5 * (1. + TMath::Erf((pTPC - fSystematicScalingEtaCorrectionMomentumThr) / 0.1));
2570 usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionLowMomenta * (1 - fractionHighMomentumScaleFactor)
2571 + fSystematicScalingEtaCorrectionHighMomenta * fractionHighMomentumScaleFactor;
2574 Double_t maxEtaVariationEl = GetMaxEtaVariation(dEdxEl);
2575 etaCorrEl = etaCorrEl * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrEl - 1.0) / maxEtaVariationEl);
2577 Double_t maxEtaVariationKa = GetMaxEtaVariation(dEdxKa);
2578 etaCorrKa = etaCorrKa * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrKa - 1.0) / maxEtaVariationKa);
2580 Double_t maxEtaVariationPi = GetMaxEtaVariation(dEdxPi);
2581 etaCorrPi = etaCorrPi * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrPi - 1.0) / maxEtaVariationPi);
2583 if (fTakeIntoAccountMuons) {
2584 Double_t maxEtaVariationMu = GetMaxEtaVariation(dEdxMu);
2585 etaCorrMu = etaCorrMu * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrMu - 1.0) / maxEtaVariationMu);
2590 Double_t maxEtaVariationPr = GetMaxEtaVariation(dEdxPr);
2591 etaCorrPr = etaCorrPr * (1.0 + (usedSystematicScalingEtaCorrection - 1.) * (etaCorrPr - 1.0) / maxEtaVariationPr);
2595 etaCorrEl = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrEl - 1.0);
2596 etaCorrKa = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrKa - 1.0);
2597 etaCorrPi = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPi - 1.0);
2598 etaCorrMu = fTakeIntoAccountMuons ? (1.0 + usedSystematicScalingEtaCorrection * (etaCorrMu - 1.0)) : 1.0;
2599 etaCorrPr = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPr - 1.0);
2603 // Get the multiplicity correction factors for the (modified) expected dEdx
2604 const Int_t currEvtMultiplicity = fPIDResponse->GetTPCResponse().GetCurrentEventMultiplicity();
2606 Double_t multiplicityCorrEl = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2607 dEdxEl * etaCorrEl, currEvtMultiplicity) : 1.;
2608 Double_t multiplicityCorrKa = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2609 dEdxKa * etaCorrKa, currEvtMultiplicity) : 1.;
2610 Double_t multiplicityCorrPi = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2611 dEdxPi * etaCorrPi, currEvtMultiplicity) : 1.;
2612 Double_t multiplicityCorrMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCMultiplicityCorrection() ?
2613 fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track, dEdxMu * etaCorrMu, currEvtMultiplicity) : 1.;
2614 Double_t multiplicityCorrPr = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2615 dEdxPr * etaCorrPr, currEvtMultiplicity) : 1.;
2617 Double_t multiplicityCorrSigmaEl = fPIDResponse->UseTPCMultiplicityCorrection() ?
2618 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxEl * etaCorrEl, currEvtMultiplicity) : 1.;
2619 Double_t multiplicityCorrSigmaKa = fPIDResponse->UseTPCMultiplicityCorrection() ?
2620 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxKa * etaCorrKa, currEvtMultiplicity) : 1.;
2621 Double_t multiplicityCorrSigmaPi = fPIDResponse->UseTPCMultiplicityCorrection() ?
2622 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxPi * etaCorrPi, currEvtMultiplicity) : 1.;
2623 Double_t multiplicityCorrSigmaMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCMultiplicityCorrection() ?
2624 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxMu * etaCorrMu, currEvtMultiplicity) : 1.;
2625 Double_t multiplicityCorrSigmaPr = fPIDResponse->UseTPCMultiplicityCorrection() ?
2626 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxPr * etaCorrPr, currEvtMultiplicity) : 1.;
2628 // Scale multiplicity correction factors, if desired (and multiplicity correction functions are to be used, otherwise it is not possible!)
2629 if (fPIDResponse->UseTPCMultiplicityCorrection() && TMath::Abs(fSystematicScalingMultCorrection - 1.0) > fgkEpsilon) {
2630 // Since we do not want to scale the splines with this, but only the multiplicity variation, only scale the deviation of the correction factor!
2631 // 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!
2633 multiplicityCorrEl = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrEl - 1.0);
2634 multiplicityCorrKa = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrKa - 1.0);
2635 multiplicityCorrPi = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrPi - 1.0);
2636 multiplicityCorrMu = fTakeIntoAccountMuons ? (1.0 + fSystematicScalingMultCorrection * (multiplicityCorrMu - 1.0)) : 1.0;
2637 multiplicityCorrPr = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrPr - 1.0);
2640 // eta correction must be enabled in order to use the new sigma parametrisation maps. Since this is the absolute sigma
2641 // for a track calculated with the unscaled paramaters, we have to devide out dEdxExpectedEtaCorrected and then need
2642 // to scale with the multiplicitySigmaCorrFactor * fSystematicScalingEtaSigmaPara. In the end, one has to scale with the
2643 // (modified) dEdx to get the absolute sigma
2644 // This means there is no extra parameter for the multiplicitySigmaCorrFactor, but only for the sigma map itself.
2645 // This is valid, since it appears only as a product. One has to assume a larger systematic shift in case of additional
2646 // multiplicity dependence....
2648 Bool_t doSigmaSystematics = (TMath::Abs(fSystematicScalingEtaSigmaParaBelowThreshold - 1.0) > fgkEpsilon) ||
2649 (TMath::Abs(fSystematicScalingEtaSigmaParaAboveThreshold - 1.0) > fgkEpsilon);
2652 Double_t dEdxElExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2654 Double_t systematicScalingEtaSigmaParaEl = 1.;
2655 if (doSigmaSystematics) {
2656 Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxElExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2657 systematicScalingEtaSigmaParaEl = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2658 + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2660 Double_t sigmaRelEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2662 / dEdxElExpected * systematicScalingEtaSigmaParaEl * multiplicityCorrSigmaEl;
2665 Double_t dEdxKaExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2667 Double_t systematicScalingEtaSigmaParaKa = 1.;
2668 if (doSigmaSystematics) {
2669 Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxKaExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2670 systematicScalingEtaSigmaParaKa = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2671 + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2673 Double_t sigmaRelKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2675 / dEdxKaExpected * systematicScalingEtaSigmaParaKa * multiplicityCorrSigmaKa;
2678 Double_t dEdxPiExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2680 Double_t systematicScalingEtaSigmaParaPi = 1.;
2681 if (doSigmaSystematics) {
2682 Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxPiExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2683 systematicScalingEtaSigmaParaPi = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2684 + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2686 Double_t sigmaRelPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2688 / dEdxPiExpected * systematicScalingEtaSigmaParaPi * multiplicityCorrSigmaPi;
2691 Double_t sigmaRelMu = 999.;
2692 if (fTakeIntoAccountMuons) {
2693 Double_t dEdxMuExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2695 Double_t systematicScalingEtaSigmaParaMu = 1.;
2696 if (doSigmaSystematics) {
2697 Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxMuExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2698 systematicScalingEtaSigmaParaMu = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2699 + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2701 sigmaRelMu = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2702 / dEdxMuExpected * systematicScalingEtaSigmaParaMu * multiplicityCorrSigmaMu;
2706 Double_t dEdxPrExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2708 Double_t systematicScalingEtaSigmaParaPr = 1.;
2709 if (doSigmaSystematics) {
2710 Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxPrExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2711 systematicScalingEtaSigmaParaPr = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2712 + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2714 Double_t sigmaRelPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2716 / dEdxPrExpected * systematicScalingEtaSigmaParaPr * multiplicityCorrSigmaPr;
2718 // Now scale the (possibly modified) spline values with the (possibly modified) correction factors
2719 dEdxEl *= etaCorrEl * multiplicityCorrEl;
2720 dEdxKa *= etaCorrKa * multiplicityCorrKa;
2721 dEdxPi *= etaCorrPi * multiplicityCorrPi;
2722 dEdxMu *= etaCorrMu * multiplicityCorrMu;
2723 dEdxPr *= etaCorrPr * multiplicityCorrPr;
2725 // Finally, get the absolute sigma
2726 sigmaEl = sigmaRelEl * dEdxEl;
2727 sigmaKa = sigmaRelKa * dEdxKa;
2728 sigmaPi = sigmaRelPi * dEdxPi;
2729 sigmaMu = sigmaRelMu * dEdxMu;
2730 sigmaPr = sigmaRelPr * dEdxPr;
2734 // No systematic studies on expected signal - just take it as it comes from the TPCPIDResponse
2735 dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2736 fPIDResponse->UseTPCEtaCorrection(),
2737 fPIDResponse->UseTPCMultiplicityCorrection());
2738 dEdxKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2739 fPIDResponse->UseTPCEtaCorrection(),
2740 fPIDResponse->UseTPCMultiplicityCorrection());
2741 dEdxPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2742 fPIDResponse->UseTPCEtaCorrection(),
2743 fPIDResponse->UseTPCMultiplicityCorrection());
2744 dEdxMu = !fTakeIntoAccountMuons ? -1 :
2745 fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2746 fPIDResponse->UseTPCEtaCorrection(),
2747 fPIDResponse->UseTPCMultiplicityCorrection());
2748 dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2749 fPIDResponse->UseTPCEtaCorrection(),
2750 fPIDResponse->UseTPCMultiplicityCorrection());
2752 sigmaEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2753 fPIDResponse->UseTPCEtaCorrection(),
2754 fPIDResponse->UseTPCMultiplicityCorrection());
2755 sigmaKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2756 fPIDResponse->UseTPCEtaCorrection(),
2757 fPIDResponse->UseTPCMultiplicityCorrection());
2758 sigmaPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2759 fPIDResponse->UseTPCEtaCorrection(),
2760 fPIDResponse->UseTPCMultiplicityCorrection());
2761 sigmaMu = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2762 fPIDResponse->UseTPCEtaCorrection(),
2763 fPIDResponse->UseTPCMultiplicityCorrection());
2764 sigmaPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2765 fPIDResponse->UseTPCEtaCorrection(),
2766 fPIDResponse->UseTPCMultiplicityCorrection());
2769 Double_t deltaPrimeElectron = (dEdxEl > 0) ? dEdxTPC / dEdxEl : -1;
2771 Printf("Error: Expected TPC signal <= 0 for electron hypothesis");
2775 Double_t deltaPrimeKaon = (dEdxKa > 0) ? dEdxTPC / dEdxKa : -1;
2777 Printf("Error: Expected TPC signal <= 0 for kaon hypothesis");
2781 Double_t deltaPrimePion = (dEdxPi > 0) ? dEdxTPC / dEdxPi : -1;
2783 Printf("Error: Expected TPC signal <= 0 for pion hypothesis");
2787 Double_t deltaPrimeProton = (dEdxPr > 0) ? dEdxTPC / dEdxPr : -1;
2789 Printf("Error: Expected TPC signal <= 0 for proton hypothesis");
2794 printf("File: %s, Line: %d: ProcessTrack -> Compute probabilities\n", (char*)__FILE__, __LINE__);
2797 // Use probabilities to weigh the response generation for the different species.
2798 // Also determine the most probable particle type.
2799 Double_t prob[AliPID::kSPECIESC];
2800 for (Int_t i = 0; i < AliPID::kSPECIESC; i++)
2803 fPIDcombined->ComputeProbabilities(track, fPIDResponse, prob);
2805 // If muons are not to be taken into account, just set their probability to zero and normalise the remaining probabilities later
2806 if (!fTakeIntoAccountMuons)
2807 prob[AliPID::kMuon] = 0;
2809 // Priors might cause problems in case of mismatches, i.e. mismatch particles will be identified as pions.
2810 // Can be easily caught for low p_TPC by just setting pion (and muon) probs to zero, if dEdx is larger than
2811 // expected dEdx of electrons and kaons
2812 if (pTPC < 1. && dEdxTPC > dEdxEl && dEdxTPC > dEdxKa) {
2813 prob[AliPID::kMuon] = 0;
2814 prob[AliPID::kPion] = 0;
2818 // Bug: One can set the number of species for PIDcombined, but PIDcombined will call PIDresponse, which writes without testing
2819 // the probs for kSPECIESC (including light nuclei) into the array.
2820 // In this case, when only kSPECIES are considered, the probabilities have to be rescaled!
2823 // 1. If the most probable PID is a light nucleus and above expected dEdx + 5 sigma of el to be sure not to mix up things in the
2824 // high-pT region (also contribution there completely negligable!)
2825 // 2. If dEdx is larger than the expected dEdx + 5 sigma of el and pr, it is most likely a light nucleus
2826 // (if dEdx is more than 5 sigma away from expectation, flat TPC probs are set... and the light nuclei splines are not
2829 // The highest expected dEdx is either for pr or for el. Assign 100% probability to the species with highest expected dEdx then.
2830 // In all other cases: Throw away light nuclei probs and rescale other probs to 100%.
2832 // Find most probable species for the ID
2833 Int_t maxIndex = TMath::LocMax(AliPID::kSPECIESC, prob);
2835 if ((prob[maxIndex] > 0 && maxIndex >= AliPID::kSPECIES && dEdxTPC > dEdxEl + 5. * sigmaEl) ||
2836 (dEdxTPC > dEdxEl + 5. * sigmaEl && dEdxTPC > dEdxPr + 5. * sigmaPr)) {
2837 for (Int_t i = 0; i < AliPID::kSPECIESC; i++)
2840 if (dEdxEl > dEdxPr) {
2841 prob[AliPID::kElectron] = 1.;
2842 maxIndex = AliPID::kElectron;
2845 prob[AliPID::kProton] = 1.;
2846 maxIndex = AliPID::kProton;
2850 for (Int_t i = AliPID::kSPECIES; i < AliPID::kSPECIESC; i++)
2853 Double_t probSum = 0.;
2854 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2858 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2862 // If all probabilities are equal (should only happen if no priors are used), set most probable PID to pion
2863 // because LocMax returns just the first maximal value (i.e. AliPID::kElectron)
2864 if (TMath::Abs(prob[maxIndex] - 1. / AliPID::kSPECIES) < 1e-5)
2865 maxIndex = AliPID::kPion;
2867 // In all other cases, the maximum remains untouched from the scaling (and is < AliPID::kSPECIES by construction)
2871 // Generate the expected responses in DeltaPrime and translate to real dEdx. Only take responses for exactly that species
2872 // (i.e. with pre-PID)
2874 Int_t numGenEntries = fgkMaxNumGenEntries; // fgkMaxNumGenEntries = 500
2876 ErrorCode errCode = kNoErrors;
2879 printf("File: %s, Line: %d: ProcessTrack -> Generate Responses for dEdx check\n", (char*)__FILE__, __LINE__);
2882 errCode = GenerateDetectorResponse(errCode, 1., sigmaEl / dEdxEl, fGenRespElDeltaPrimeEl, numGenEntries);
2885 errCode = GenerateDetectorResponse(errCode, 1., sigmaKa / dEdxKa, fGenRespKaDeltaPrimeKa, numGenEntries);
2888 errCode = GenerateDetectorResponse(errCode, 1., sigmaPi / dEdxPi, fGenRespPiDeltaPrimePi, numGenEntries);
2891 errCode = GenerateDetectorResponse(errCode, 1., sigmaPr / dEdxPr, fGenRespPrDeltaPrimePr, numGenEntries);
2893 if (errCode == kNoErrors) {
2894 Double_t genEntry[kDeDxCheckNumAxes];
2895 genEntry[kDeDxCheckJetPt] = jetPt;
2896 genEntry[kDeDxCheckEtaAbs] = TMath::Abs(track->Eta());
2897 genEntry[kDeDxCheckP] = pTPC;
2900 for (Int_t n = 0; n < numGenEntries; n++) {
2901 // If no MC info is available or shall not be used, use weighting with priors to generate entries for the different species
2902 Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
2904 // Consider generated response as originating from...
2905 if (rnd <= prob[AliPID::kElectron]) {
2906 genEntry[kDeDxCheckPID] = 0; // ... an electron
2907 genEntry[kDeDxCheckDeDx] = fGenRespElDeltaPrimeEl[n] * dEdxEl;
2909 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon]) {
2910 genEntry[kDeDxCheckPID] = 1; // ... a kaon
2911 genEntry[kDeDxCheckDeDx] = fGenRespKaDeltaPrimeKa[n] * dEdxKa;
2913 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon] + prob[AliPID::kPion]) {
2914 genEntry[kDeDxCheckPID] = 2; // ... a pion (or a muon)
2915 genEntry[kDeDxCheckDeDx] = fGenRespPiDeltaPrimePi[n] * dEdxPi;
2918 genEntry[kDeDxCheckPID] = 3; // ... a proton
2919 genEntry[kDeDxCheckDeDx] = fGenRespPrDeltaPrimePr[n] * dEdxPr;
2922 fDeDxCheck->Fill(genEntry);
2927 printf("File: %s, Line: %d: ProcessTrack -> Generate Responses for dEdx check done\n", (char*)__FILE__, __LINE__);
2930 if (fDoPtResolution) {
2931 // Check shared clusters, which is done together with the pT resolution
2932 Double_t qaEntry[kQASharedClsNumAxes];
2933 qaEntry[kQASharedClsJetPt] = jetPt;
2934 qaEntry[kQASharedClsPt] = pT;
2935 qaEntry[kDeDxCheckP] = pTPC;
2936 qaEntry[kQASharedClsNumSharedCls] = track->GetTPCSharedMapPtr()->CountBits();
2939 // iRowInd == -1 for "all rows w/o multiple counting"
2940 qaEntry[kQASharedClsPadRow] = iRowInd;
2941 fQASharedCls->Fill(qaEntry);
2943 // Fill hist for every pad row with shared cluster
2944 for (iRowInd = 0; iRowInd < 159; iRowInd++) {
2945 if (track->GetTPCSharedMapPtr()->TestBitNumber(iRowInd)) {
2946 qaEntry[kQASharedClsPadRow] = iRowInd;
2947 fQASharedCls->Fill(qaEntry);
2956 // Clean up the most probable PID for low momenta (not an issue for the probabilities, since fractions small,
2957 // but to get clean (and nice looking) templates (if most probable PID is used instead of generated responses), it should be done).
2958 // Problem: If more than 5 sigma away (and since priors also included), most probable PID for dEdx around protons could be pions.
2959 // Idea: Only clean selection for K and p possible. So, just take for the most probable PID the probs from TPC only calculated
2960 // by hand without restriction to 5 sigma and without priors. This will not affect the template generation, but afterwards one
2961 // can replace the K and p templates with the most probable PID distributions, so all other species templates remain the same.
2962 // I.e. it doesn't hurt if the most probable PID is distributed incorrectly between el, pi, mu.
2963 Bool_t maxIndexSetManually = kFALSE;
2966 Double_t probManualTPC[AliPID::kSPECIES];
2967 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2968 probManualTPC[i] = 0;
2970 probManualTPC[AliPID::kElectron] = TMath::Exp(-0.5*(dEdxTPC-dEdxEl)*(dEdxTPC-dEdxEl)/(sigmaEl*sigmaEl))/sigmaEl;
2971 // Muons are not important here, just ignore and save processing time
2972 probManualTPC[AliPID::kMuon] = 0.;//TMath::Exp(-0.5*(dEdxTPC-dEdxMu)*(dEdxTPC-dEdxMu)/(sigmaMu*sigmaMu))/sigmaMu;
2973 probManualTPC[AliPID::kPion] = TMath::Exp(-0.5*(dEdxTPC-dEdxPi)*(dEdxTPC-dEdxPi)/(sigmaPi*sigmaPi))/sigmaPi;
2974 probManualTPC[AliPID::kKaon] = TMath::Exp(-0.5*(dEdxTPC-dEdxKa)*(dEdxTPC-dEdxKa)/(sigmaKa*sigmaKa))/sigmaKa;
2975 probManualTPC[AliPID::kProton] = TMath::Exp(-0.5*(dEdxTPC-dEdxPr)*(dEdxTPC-dEdxPr)/(sigmaPr*sigmaPr))/sigmaPr;
2977 const Int_t maxIndexManualTPC = TMath::LocMax(AliPID::kSPECIES, probManualTPC);
2978 // Should always be > 0, but in case the dEdx is far off such that the machine accuracy sets every prob to zero,
2979 // better take the "old" result
2980 if (probManualTPC[maxIndexManualTPC] > 0.)
2981 maxIndex = maxIndexManualTPC;
2983 maxIndexSetManually = kTRUE;
2987 // Translate from AliPID numbering to numbering of this class
2988 if (prob[maxIndex] > 0 || maxIndexSetManually) {
2989 if (maxIndex == AliPID::kElectron)
2991 else if (maxIndex == AliPID::kKaon)
2993 else if (maxIndex == AliPID::kMuon)
2995 else if (maxIndex == AliPID::kPion)
2997 else if (maxIndex == AliPID::kProton)
3003 // Only take track into account for expectation values, if valid pid response is available.. Otherwise: Set to underflow bin.
3009 //For testing: Swap p<->pT to analyse pure p-dependence => Needs to be removed later
3015 TOFpidInfo tofPIDinfo = GetTOFType(track, fTOFmode);
3017 Double_t entry[fStoreAdditionalJetInformation ? kDataNumAxes : kDataNumAxes - fgkNumJetAxes];
3018 entry[kDataMCID] = binMC;
3019 entry[kDataSelectSpecies] = 0;
3020 entry[kDataPt] = pT;
3021 entry[kDataDeltaPrimeSpecies] = deltaPrimeElectron;
3022 entry[kDataCentrality] = centralityPercentile;
3024 if (fStoreAdditionalJetInformation) {
3025 entry[kDataJetPt] = jetPt;
3027 entry[kDataXi] = xi;
3030 entry[GetIndexOfChargeAxisData()] = trackCharge;
3031 entry[GetIndexOfTOFpidInfoAxisData()] = tofPIDinfo;
3033 fhPIDdataAll->Fill(entry);
3035 entry[kDataSelectSpecies] = 1;
3036 entry[kDataDeltaPrimeSpecies] = deltaPrimeKaon;
3037 fhPIDdataAll->Fill(entry);
3039 entry[kDataSelectSpecies] = 2;
3040 entry[kDataDeltaPrimeSpecies] = deltaPrimePion;
3041 fhPIDdataAll->Fill(entry);
3043 entry[kDataSelectSpecies] = 3;
3044 entry[kDataDeltaPrimeSpecies] = deltaPrimeProton;
3045 fhPIDdataAll->Fill(entry);
3048 // Construct the expected shape for the transition p -> pT
3050 Double_t genEntry[fStoreAdditionalJetInformation ? kGenNumAxes : kGenNumAxes - fgkNumJetAxes];
3051 genEntry[kGenMCID] = binMC;
3052 genEntry[kGenSelectSpecies] = 0;
3053 genEntry[kGenPt] = pT;
3054 genEntry[kGenDeltaPrimeSpecies] = -999;
3055 genEntry[kGenCentrality] = centralityPercentile;
3057 if (fStoreAdditionalJetInformation) {
3058 genEntry[kGenJetPt] = jetPt;
3059 genEntry[kGenZ] = z;
3060 genEntry[kGenXi] = xi;
3063 genEntry[GetIndexOfChargeAxisGen()] = trackCharge;
3064 genEntry[GetIndexOfTOFpidInfoAxisGen()] = tofPIDinfo;
3066 // Generate numGenEntries "responses" with fluctuations around the expected values.
3067 // fgkMaxNumGenEntries = 500 turned out to give reasonable templates even for highest track pT in 15-20 GeV/c jet pT bin.
3068 Int_t numGenEntries = fgkMaxNumGenEntries; // fgkMaxNumGenEntries = 500
3070 /*OLD: Different number of responses depending on pT and jet pT for fgkMaxNumGenEntries = 1000
3071 * => Problem: If threshold to higher number of responses inside a bin (or after rebinning), then the template
3072 * is biased to the higher pT.
3073 // Generate numGenEntries "responses" with fluctuations around the expected values.
3074 // The higher the (transverse) momentum, the more "responses" will be generated in order not to run out of statistics too fast.
3075 Int_t numGenEntries = 10;
3077 // Jets have even worse statistics, therefore, scale numGenEntries further
3079 numGenEntries *= 20;
3080 else if (jetPt >= 20)
3081 numGenEntries *= 10;
3082 else if (jetPt >= 10)
3087 // Do not generate more entries than available in memory!
3088 if (numGenEntries > fgkMaxNumGenEntries)// fgkMaxNumGenEntries = 1000
3089 numGenEntries = fgkMaxNumGenEntries;
3093 ErrorCode errCode = kNoErrors;
3096 printf("File: %s, Line: %d: ProcessTrack -> Generate Responses\n", (char*)__FILE__, __LINE__);
3099 errCode = GenerateDetectorResponse(errCode, 1., sigmaEl / dEdxEl, fGenRespElDeltaPrimeEl, numGenEntries);
3100 errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxKa, sigmaEl / dEdxKa, fGenRespElDeltaPrimeKa, numGenEntries);
3101 errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxPi, sigmaEl / dEdxPi, fGenRespElDeltaPrimePi, numGenEntries);
3102 errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxPr, sigmaEl / dEdxPr, fGenRespElDeltaPrimePr, numGenEntries);
3105 errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxEl, sigmaKa / dEdxEl, fGenRespKaDeltaPrimeEl, numGenEntries);
3106 errCode = GenerateDetectorResponse(errCode, 1., sigmaKa / dEdxKa, fGenRespKaDeltaPrimeKa, numGenEntries);
3107 errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxPi, sigmaKa / dEdxPi, fGenRespKaDeltaPrimePi, numGenEntries);
3108 errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxPr, sigmaKa / dEdxPr, fGenRespKaDeltaPrimePr, numGenEntries);
3111 errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxEl, sigmaPi / dEdxEl, fGenRespPiDeltaPrimeEl, numGenEntries);
3112 errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxKa, sigmaPi / dEdxKa, fGenRespPiDeltaPrimeKa, numGenEntries);
3113 errCode = GenerateDetectorResponse(errCode, 1., sigmaPi / dEdxPi, fGenRespPiDeltaPrimePi, numGenEntries);
3114 errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxPr, sigmaPi / dEdxPr, fGenRespPiDeltaPrimePr, numGenEntries);
3116 // Muons, if desired
3117 if (fTakeIntoAccountMuons) {
3118 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxEl, sigmaMu / dEdxEl, fGenRespMuDeltaPrimeEl, numGenEntries);
3119 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxKa, sigmaMu / dEdxKa, fGenRespMuDeltaPrimeKa, numGenEntries);
3120 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxPi, sigmaMu / dEdxPi, fGenRespMuDeltaPrimePi, numGenEntries);
3121 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxPr, sigmaMu / dEdxPr, fGenRespMuDeltaPrimePr, numGenEntries);
3125 errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxEl, sigmaPr / dEdxEl, fGenRespPrDeltaPrimeEl, numGenEntries);
3126 errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxKa, sigmaPr / dEdxKa, fGenRespPrDeltaPrimeKa, numGenEntries);
3127 errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxPi, sigmaPr / dEdxPi, fGenRespPrDeltaPrimePi, numGenEntries);
3128 errCode = GenerateDetectorResponse(errCode, 1., sigmaPr / dEdxPr, fGenRespPrDeltaPrimePr, numGenEntries);
3130 if (errCode != kNoErrors) {
3131 if (errCode == kWarning && fDebug > 1) {
3132 Printf("Warning: Questionable detector response for track, most likely due to very low number of PID clusters! Debug output (dEdx_expected, sigma_expected):");
3135 Printf("Error: Failed to generate detector response for track - skipped! Debug output (dEdx_expected, sigma_expected):");
3138 Printf("Pr: %e, %e", dEdxPr, sigmaPr);
3139 Printf("Pi: %e, %e", dEdxPi, sigmaPi);
3140 Printf("El: %e, %e", dEdxEl, sigmaEl);
3141 Printf("Mu: %e, %e", dEdxMu, sigmaMu);
3142 Printf("Ka: %e, %e", dEdxKa, sigmaKa);
3143 Printf("track: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), track->GetTPCmomentum(), track->Eta(),
3144 track->GetTPCsignalN());
3147 if (errCode != kWarning) {
3148 return kFALSE;// Skip generated response in case of error
3152 for (Int_t n = 0; n < numGenEntries; n++) {
3153 if (!isMC || !fUseMCidForGeneration) {
3154 // If no MC info is available or shall not be used, use weighting with priors to generate entries for the different species
3155 Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
3157 // Consider generated response as originating from...
3158 if (rnd <= prob[AliPID::kElectron])
3159 genEntry[kGenMCID] = 0; // ... an electron
3160 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon])
3161 genEntry[kGenMCID] = 1; // ... a kaon
3162 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon])
3163 genEntry[kGenMCID] = 2; // ... a muon -> NOTE: prob[AliPID::kMuon] = 0 in case of fTakeIntoAccountMuons = kFALSE
3164 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon] + prob[AliPID::kPion])
3165 genEntry[kGenMCID] = 3; // ... a pion
3167 genEntry[kGenMCID] = 4; // ... a proton
3171 genEntry[kGenSelectSpecies] = 0;
3172 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeEl[n];
3173 fhGenEl->Fill(genEntry);
3175 genEntry[kGenSelectSpecies] = 1;
3176 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeKa[n];
3177 fhGenEl->Fill(genEntry);
3179 genEntry[kGenSelectSpecies] = 2;
3180 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePi[n];
3181 fhGenEl->Fill(genEntry);
3183 genEntry[kGenSelectSpecies] = 3;
3184 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePr[n];
3185 fhGenEl->Fill(genEntry);
3188 genEntry[kGenSelectSpecies] = 0;
3189 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeEl[n];
3190 fhGenKa->Fill(genEntry);
3192 genEntry[kGenSelectSpecies] = 1;
3193 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeKa[n];
3194 fhGenKa->Fill(genEntry);
3196 genEntry[kGenSelectSpecies] = 2;
3197 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePi[n];
3198 fhGenKa->Fill(genEntry);
3200 genEntry[kGenSelectSpecies] = 3;
3201 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePr[n];
3202 fhGenKa->Fill(genEntry);
3205 genEntry[kGenSelectSpecies] = 0;
3206 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeEl[n];
3207 fhGenPi->Fill(genEntry);
3209 genEntry[kGenSelectSpecies] = 1;
3210 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeKa[n];
3211 fhGenPi->Fill(genEntry);
3213 genEntry[kGenSelectSpecies] = 2;
3214 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePi[n];
3215 fhGenPi->Fill(genEntry);
3217 genEntry[kGenSelectSpecies] = 3;
3218 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePr[n];
3219 fhGenPi->Fill(genEntry);
3221 if (fTakeIntoAccountMuons) {
3222 // Muons, if desired
3223 genEntry[kGenSelectSpecies] = 0;
3224 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeEl[n];
3225 fhGenMu->Fill(genEntry);
3227 genEntry[kGenSelectSpecies] = 1;
3228 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeKa[n];
3229 fhGenMu->Fill(genEntry);
3231 genEntry[kGenSelectSpecies] = 2;
3232 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePi[n];
3233 fhGenMu->Fill(genEntry);
3235 genEntry[kGenSelectSpecies] = 3;
3236 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePr[n];
3237 fhGenMu->Fill(genEntry);
3241 genEntry[kGenSelectSpecies] = 0;
3242 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeEl[n];
3243 fhGenPr->Fill(genEntry);
3245 genEntry[kGenSelectSpecies] = 1;
3246 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeKa[n];
3247 fhGenPr->Fill(genEntry);
3249 genEntry[kGenSelectSpecies] = 2;
3250 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePi[n];
3251 fhGenPr->Fill(genEntry);
3253 genEntry[kGenSelectSpecies] = 3;
3254 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePr[n];
3255 fhGenPr->Fill(genEntry);
3259 printf("File: %s, Line: %d: ProcessTrack -> Done\n", (char*)__FILE__, __LINE__);
3265 //_____________________________________________________________________________
3266 Bool_t AliAnalysisTaskPID::SetConvolutedGaussLambdaParameter(Double_t lambda)
3268 // Set the lambda parameter of the convolution to the desired value. Automatically
3269 // calculates the parameters for the transition (restricted) gauss -> convoluted gauss.
3270 fConvolutedGaussTransitionPars[2] = lambda;
3272 // Save old parameters and settings of function to restore them later:
3273 Double_t* oldFuncParams = new Double_t[fkConvolutedGausNPar];
3274 for (Int_t i = 0; i < fkConvolutedGausNPar; i++)
3275 oldFuncParams[i] = fConvolutedGausDeltaPrime->GetParameter(i);
3276 Int_t oldFuncNpx = fConvolutedGausDeltaPrime->GetNpx();
3277 Double_t oldFuncRangeLow = 0, oldFuncRangeUp = 100;
3278 fConvolutedGausDeltaPrime->GetRange(oldFuncRangeLow, oldFuncRangeUp);
3280 // Choose some sufficiently large range
3281 const Double_t rangeStart = 0.5;
3282 const Double_t rangeEnd = 2.0;
3284 // To get the parameters for the transition, just choose arbitrarily input parameters for mu and sigma
3285 // (it makes sense to choose typical values). The ratio sigma_gauss / sigma_convolution is independent
3286 // of mu and as well the difference mu_gauss - mu_convolution.
3287 Double_t muInput = 1.0;
3288 Double_t sigmaInput = fgkSigmaReferenceForTransitionPars;
3291 // Step 1: Generate distribution with input parameters
3292 const Int_t nPar = 3;
3293 Double_t inputPar[nPar] = { muInput, sigmaInput, lambda };
3295 TH1D* hInput = new TH1D("hInput", "Input distribution", 2000, rangeStart, rangeEnd);
3297 fConvolutedGausDeltaPrime->SetParameters(inputPar);
3298 fConvolutedGausDeltaPrime->SetRange(rangeStart, rangeEnd);
3299 fConvolutedGausDeltaPrime->SetNpx(2000);
3302 // The resolution and mean of the AliPIDResponse are determined in finite intervals
3303 // of ncl (also second order effects due to finite dEdx and tanTheta).
3304 // Take this into account for the transition parameters via assuming an approximately flat
3305 // distritubtion in ncl in this interval.
3306 // NOTE: The ncl interval should be the same as the one used for the sigma map creation!
3307 const Int_t nclStart = 151;
3308 const Int_t nclEnd = 160;
3309 const Int_t nclSteps = (nclEnd - nclStart) + 1;
3310 for (Int_t ncl = nclStart; ncl <= nclEnd; ncl++) {
3311 // Resolution scales with 1/sqrt(ncl)
3312 fConvolutedGausDeltaPrime->SetParameter(1, inputPar[1] * sqrt(nclEnd) / sqrt(ncl));
3313 TH1* hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
3315 for (Int_t i = 0; i < 50000000 / nclSteps; i++)
3316 hInput->Fill(hProbDensity->GetRandom());
3320 TH1* hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
3322 for (Int_t i = 0; i < 50000000; i++)
3323 hInput->Fill(hProbDensity->GetRandom());
3325 // Step 2: Fit generated distribution with restricted gaussian
3326 Int_t maxBin = hInput->GetMaximumBin();
3327 Double_t max = hInput->GetBinContent(maxBin);
3329 UChar_t usedBins = 1;
3331 max += hInput->GetBinContent(maxBin - 1);
3334 if (maxBin < hInput->GetNbinsX()) {
3335 max += hInput->GetBinContent(maxBin + 1);
3340 // NOTE: The range (<-> fraction of maximum) should be the same
3341 // as for the spline and eta maps creation
3342 const Double_t lowThreshold = hInput->GetXaxis()->GetBinLowEdge(hInput->FindFirstBinAbove(0.1 * max));
3343 const Double_t highThreshold = hInput->GetXaxis()->GetBinUpEdge(hInput->FindLastBinAbove(0.1 * max));
3345 TFitResultPtr fitResGaussFirstStep = hInput->Fit("gaus", "QNRS", "", lowThreshold, highThreshold);
3347 TFitResultPtr fitResGauss;
3349 if ((Int_t)fitResGaussFirstStep == 0) {
3350 TF1 fGauss("fGauss", "[0]*TMath::Gaus(x, [1], [2], 1)", rangeStart, rangeEnd);
3351 fGauss.SetParameter(0, fitResGaussFirstStep->GetParams()[0]);
3352 fGauss.SetParError(0, fitResGaussFirstStep->GetErrors()[0]);
3353 fGauss.SetParameter(2, fitResGaussFirstStep->GetParams()[2]);
3354 fGauss.SetParError(2, fitResGaussFirstStep->GetErrors()[2]);
3356 fGauss.FixParameter(1, fitResGaussFirstStep->GetParams()[1]);
3357 fitResGauss = hInput->Fit(&fGauss, "QNS", "s", rangeStart, rangeEnd);
3360 fitResGauss = hInput->Fit("gaus", "QNRS", "same", rangeStart, rangeEnd);
3362 //OLD TFitResultPtr fitResGauss = hInput->Fit("gaus", "QNRS", "", hInput->GetXaxis()->GetBinLowEdge(hInput->FindFirstBinAbove(0.1 * max)),
3363 // hInput->GetXaxis()->GetBinUpEdge(hInput->FindLastBinAbove(0.1 * max)));
3366 // Step 3: Use parameters from gaussian fit to obtain parameters for the transition "restricted gauss" -> "convoluted gauss"
3368 // 3.1 The ratio sigmaInput / sigma_gaussFit ONLY depends on lambda (which is fixed per period) -> Calculate this first
3369 // for an arbitrary (here: typical) sigma. The ratio is then ~the same for ALL sigma for given lambda!
3370 if ((Int_t)fitResGauss != 0) {
3371 AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Gauss Fit failed!\n");
3372 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
3373 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
3374 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
3377 delete [] oldFuncParams;
3382 if (fitResGauss->GetParams()[2] <= 0) {
3383 AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Sigma of gauss fit <= 0!\n");
3384 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
3385 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
3386 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
3389 delete [] oldFuncParams;
3394 // sigma correction factor
3395 fConvolutedGaussTransitionPars[1] = sigmaInput / fitResGauss->GetParams()[2];
3397 // 3.2 Now that sigma und lambda are determined, one can calculate mu by shifting the maximum to the desired position,
3398 // i.e. the maximum of the input distribution should coincide with that of the re-calculated distribution,
3399 // which is achieved by getting the same mu for the same sigma.
3400 // NOTE: For fixed lambda, the shift is proportional to sigma and independent of mu!
3401 // So, one can calculate the shift for an arbitrary fixed (here: typical)
3402 // sigma and then simply use this shift for any other sigma by scaling it correspondingly!!!
3404 // Mu shift correction:
3405 // Shift in mu (difference between mean of gauss and mean of convolution) is proportional to sigma!
3406 // Thus, choose a reference sigma (typical -> 0.05), such that for arbitrary sigma one can simple scale
3407 // this shift correction with sigma / referenceSigma.
3408 fConvolutedGaussTransitionPars[0] = (fitResGauss->GetParams()[1] - muInput);
3411 /*Changed on 03.07.2013
3413 // Maximum of fConvolutedGausDeltaPrime should agree with maximum of input
3414 Double_t par[nPar] = { fitResGauss->GetParams()[1], // just as a guess of the maximum position
3418 fConvolutedGausDeltaPrime->SetParameters(par);
3420 Double_t maxXInput = fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001, muInput - 3. * sigmaInput),
3421 muInput + 10. * sigmaInput,
3424 // Maximum shifts proportional to sigma and is linear in mu (~mean of gauss)
3425 // Maximum should be typically located within [gaussMean, gaussMean + 3 gaussSigma].
3426 // -> Larger search range for safety reasons (also: sigma and/or mean might be not 100% accurate).
3427 Double_t maxXconvoluted = fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001,
3428 fitResGauss->GetParams()[1] - 3. * fitResGauss->GetParams()[2]),
3429 fitResGauss->GetParams()[1] + 10. * fitResGauss->GetParams()[2],
3431 if (maxXconvoluted <= 0) {
3432 AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Maximum of fConvolutedGausDeltaPrime <= 0!\n");
3434 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
3435 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
3436 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
3439 delete [] oldFuncParams;
3444 // maxX perfectly shifts as par[0] (scaled by sigma) -> Can shift maxX to input value.
3445 // Mu shift correction:
3446 fConvolutedGaussTransitionPars[0] = maxXconvoluted - maxXInput;
3451 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
3452 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
3453 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
3456 delete [] oldFuncParams;
3462 //_____________________________________________________________________________
3463 AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::SetParamsForConvolutedGaus(Double_t gausMean, Double_t gausSigma)
3465 // Set parameters for convoluted gauss using parameters for a pure gaussian.
3466 // If SetConvolutedGaussLambdaParameter has not been called before to initialise the translation parameters,
3467 // some default parameters will be used and an error will show up.
3470 printf("File: %s, Line: %d: SetParamsForConvolutedGaus: mean %e, sigma %e\n", (char*)__FILE__, __LINE__, gausMean, gausSigma);
3472 if (fConvolutedGaussTransitionPars[1] < -998) {
3473 AliError("Transition parameters not initialised! Default parameters will be used. Please call SetConvolutedGaussLambdaParameter(...) before any calculations!");
3474 SetConvolutedGaussLambdaParameter(2.0);
3475 AliError(Form("Parameters set to:\n[0]: %f\n[1]: %f\n[2]: %f\n", fConvolutedGaussTransitionPars[0],
3476 fConvolutedGaussTransitionPars[1], fConvolutedGaussTransitionPars[2]));
3479 Double_t par[fkConvolutedGausNPar];
3480 par[2] = fConvolutedGaussTransitionPars[2];
3481 par[1] = fConvolutedGaussTransitionPars[1] * gausSigma;
3482 // maxX perfectly shifts as par[0] (scaled by sigma) -> Can shift maxX so that it sits at the right place.
3483 par[0] = gausMean - fConvolutedGaussTransitionPars[0] * par[1] / fgkSigmaReferenceForTransitionPars;
3485 ErrorCode errCode = kNoErrors;
3486 fConvolutedGausDeltaPrime->SetParameters(par);
3489 printf("File: %s, Line: %d: SetParamsForConvolutedGaus -> Parameters set to: %e, %e, %e (transition pars: %e, %e, %e, %e)\n",
3490 (char*)__FILE__, __LINE__, par[0], par[1], par[2], fConvolutedGaussTransitionPars[0], fConvolutedGaussTransitionPars[1],
3491 fConvolutedGaussTransitionPars[2], fgkSigmaReferenceForTransitionPars);
3493 fConvolutedGausDeltaPrime->SetNpx(20); // Small value speeds up following algorithm (valid, since extrema far apart)
3495 // Accuracy of 10^-5 is enough to get 0.1% precise peak for MIPS w.r.t. to dEdx = 2000 of protons
3496 // (should boost up the algorithm, because 10^-10 is the default value!)
3497 Double_t maxX= fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001, gausMean - 2. * gausSigma),
3498 gausMean + 6. * gausSigma, 1.0E-5);
3500 const Double_t maximum = fConvolutedGausDeltaPrime->Eval(maxX);
3501 const Double_t maximumFraction = maximum * fAccuracyNonGaussianTail;
3503 // Estimate lower boundary for subsequent search:
3504 Double_t lowBoundSearchBoundLow = TMath::Max(1e-4, maxX - 5. * gausSigma);
3505 Double_t lowBoundSearchBoundUp = maxX;
3507 Bool_t lowerBoundaryFixedAtZero = kFALSE;
3509 while (fConvolutedGausDeltaPrime->Eval(lowBoundSearchBoundLow) >= maximumFraction) {
3510 if (lowBoundSearchBoundLow <= 0) {
3511 // This should only happen to low dEdx particles with very few clusters and therefore large sigma, such that the gauss goes below zero deltaPrime
3512 if (maximum <= 0) { // Something is weired
3513 printf("Error generating signal: maximum is <= 0!\n");
3517 const Double_t valueAtZero = fConvolutedGausDeltaPrime->Eval(0);
3518 if (valueAtZero / maximum > 0.05) {
3519 // Too large fraction below zero deltaPrime. Signal generation cannot be reliable in this case
3520 printf("Error generating signal: Too large fraction below zero deltaPrime: convGauss(0) / convGauss(max) = %e / %e = %e!\n",
3521 valueAtZero, maximum, valueAtZero / maximum);
3527 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",
3528 fConvolutedGausDeltaPrime->Eval(maxX), fAccuracyNonGaussianTail, maximumFraction, maxX, par[0], par[1], par[2], gausMean, gausSigma);
3531 lowerBoundaryFixedAtZero = kTRUE;
3533 if (errCode != kError)
3539 lowBoundSearchBoundUp -= gausSigma;
3540 lowBoundSearchBoundLow -= gausSigma;
3542 if (lowBoundSearchBoundLow < 0) {
3543 lowBoundSearchBoundLow = 0;
3544 lowBoundSearchBoundUp += gausSigma;
3548 // Determine lower boundary inside estimated range. For small values of the maximum: Need more precision, since finer binning!
3549 Double_t rangeStart = lowerBoundaryFixedAtZero ? 0 :
3550 fConvolutedGausDeltaPrime->GetX(maximumFraction, lowBoundSearchBoundLow, lowBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
3552 // .. and the same for the upper boundary
3553 Double_t rangeEnd = 0;
3554 // If distribution starts beyond upper boundary, everything ends up in the overflow bin. So, just reduce range and Npx to minimum
3555 if (rangeStart > fkDeltaPrimeUpLimit) {
3556 rangeEnd = rangeStart + 0.00001;
3557 fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
3558 fConvolutedGausDeltaPrime->SetNpx(4);
3561 // Estimate upper boundary for subsequent search:
3562 Double_t upBoundSearchBoundUp = maxX + 5 * gausSigma;
3563 Double_t upBoundSearchBoundLow = maxX;
3564 while (fConvolutedGausDeltaPrime->Eval(upBoundSearchBoundUp) >= maximumFraction) {
3565 upBoundSearchBoundUp += gausSigma;
3566 upBoundSearchBoundLow += gausSigma;
3569 // For small values of the maximum: Need more precision, since finer binning!
3570 rangeEnd = fConvolutedGausDeltaPrime->GetX(maximumFraction, upBoundSearchBoundLow, upBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
3572 fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
3573 fConvolutedGausDeltaPrime->SetNpx(fDeltaPrimeAxis->FindFixBin(rangeEnd) - fDeltaPrimeAxis->FindFixBin(rangeStart) + 1);
3575 fConvolutedGausDeltaPrime->SetNpx(fDeltaPrimeAxis->FindFixBin(rangeEnd) - fDeltaPrimeAxis->FindFixBin(rangeStart) + 1);
3576 //fConvolutedGausDeltaPrime->SetNpx(fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindFixBin(rangeEnd)
3577 // - fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindFixBin(rangeStart) + 1);
3578 //fConvolutedGausDeltaPrime->SetNpx((rangeEnd - rangeStart) / fDeltaPrimeBinWidth + 1);
3582 printf("File: %s, Line: %d: SetParamsForConvolutedGaus -> range %f - %f, error code %d\n", (char*)__FILE__, __LINE__,
3583 rangeStart, rangeEnd, errCode);
3589 //________________________________________________________________________
3590 void AliAnalysisTaskPID::SetUpGenHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const
3592 // Sets bin limits for axes which are not standard binned and the axes titles.
3594 hist->SetBinEdges(kGenPt, binsPt);
3595 hist->SetBinEdges(kGenDeltaPrimeSpecies, binsDeltaPrime);
3596 hist->SetBinEdges(kGenCentrality, binsCent);
3598 if (fStoreAdditionalJetInformation)
3599 hist->SetBinEdges(kGenJetPt, binsJetPt);
3602 hist->GetAxis(kGenMCID)->SetTitle("MC PID");
3603 hist->GetAxis(kGenMCID)->SetBinLabel(1, "e");
3604 hist->GetAxis(kGenMCID)->SetBinLabel(2, "K");
3605 hist->GetAxis(kGenMCID)->SetBinLabel(3, "#mu");
3606 hist->GetAxis(kGenMCID)->SetBinLabel(4, "#pi");
3607 hist->GetAxis(kGenMCID)->SetBinLabel(5, "p");
3609 hist->GetAxis(kGenSelectSpecies)->SetTitle("Select Species");
3610 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(1, "e");
3611 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(2, "K");
3612 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(3, "#pi");
3613 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(4, "p");
3615 hist->GetAxis(kGenPt)->SetTitle("p_{T} (GeV/c)");
3617 hist->GetAxis(kGenDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
3619 hist->GetAxis(kGenCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
3621 if (fStoreAdditionalJetInformation) {
3622 hist->GetAxis(kGenJetPt)->SetTitle("p_{T}^{jet} (GeV/c)");
3624 hist->GetAxis(kGenZ)->SetTitle("z = p_{T}^{track} / p_{T}^{jet}");
3626 hist->GetAxis(kGenXi)->SetTitle("#xi = ln(p_{T}^{jet} / p_{T}^{track})");
3629 hist->GetAxis(GetIndexOfChargeAxisGen())->SetTitle("Charge (e_{0})");
3631 hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetTitle("TOF PID Info");
3632 // Offset is (TMath::Abs(kNoTOFinfo) + 1), such that bin 1 (= first label) corresponds to kNoTOFinfo (< 0)
3633 hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kNoTOFinfo + (TMath::Abs(kNoTOFinfo) + 1), "No TOF Info");
3634 hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kNoTOFpid + (TMath::Abs(kNoTOFinfo) + 1), "No TOF PID");
3635 hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kTOFpion + (TMath::Abs(kNoTOFinfo) + 1), "#pi");
3636 hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kTOFkaon + (TMath::Abs(kNoTOFinfo) + 1), "K");
3637 hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kTOFproton + (TMath::Abs(kNoTOFinfo) + 1), "p");
3641 //________________________________________________________________________
3642 void AliAnalysisTaskPID::SetUpGenYieldHist(THnSparse* hist, Double_t* binsPt, Double_t* binsCent, Double_t* binsJetPt) const
3644 // Sets bin limits for axes which are not standard binned and the axes titles.
3646 hist->SetBinEdges(kGenYieldPt, binsPt);
3647 hist->SetBinEdges(kGenYieldCentrality, binsCent);
3648 if (fStoreAdditionalJetInformation)
3649 hist->SetBinEdges(kGenYieldJetPt, binsJetPt);
3651 for (Int_t i = 0; i < 5; i++)
3652 hist->GetAxis(kGenYieldMCID)->SetBinLabel(i + 1, AliPID::ParticleLatexName(i));
3655 hist->GetAxis(kGenYieldMCID)->SetTitle("MC PID");
3656 hist->GetAxis(kGenYieldPt)->SetTitle("p_{T}^{gen} (GeV/c)");
3657 hist->GetAxis(kGenYieldCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
3659 if (fStoreAdditionalJetInformation) {
3660 hist->GetAxis(kGenYieldJetPt)->SetTitle("p_{T}^{jet, gen} (GeV/c)");
3662 hist->GetAxis(kGenYieldZ)->SetTitle("z = p_{T}^{track} / p_{T}^{jet}");
3664 hist->GetAxis(kGenYieldXi)->SetTitle("#xi = ln(p_{T}^{jet} / p_{T}^{track})");
3667 hist->GetAxis(GetIndexOfChargeAxisGenYield())->SetTitle("Charge (e_{0})");
3671 //________________________________________________________________________
3672 void AliAnalysisTaskPID::SetUpHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const
3674 // Sets bin limits for axes which are not standard binned and the axes titles.
3676 hist->SetBinEdges(kDataPt, binsPt);
3677 hist->SetBinEdges(kDataDeltaPrimeSpecies, binsDeltaPrime);
3678 hist->SetBinEdges(kDataCentrality, binsCent);
3680 if (fStoreAdditionalJetInformation)
3681 hist->SetBinEdges(kDataJetPt, binsJetPt);
3684 hist->GetAxis(kDataMCID)->SetTitle("MC PID");
3685 hist->GetAxis(kDataMCID)->SetBinLabel(1, "e");
3686 hist->GetAxis(kDataMCID)->SetBinLabel(2, "K");
3687 hist->GetAxis(kDataMCID)->SetBinLabel(3, "#mu");
3688 hist->GetAxis(kDataMCID)->SetBinLabel(4, "#pi");
3689 hist->GetAxis(kDataMCID)->SetBinLabel(5, "p");
3691 hist->GetAxis(kDataSelectSpecies)->SetTitle("Select Species");
3692 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(1, "e");
3693 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(2, "K");
3694 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(3, "#pi");
3695 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(4, "p");
3697 hist->GetAxis(kDataPt)->SetTitle("p_{T} (GeV/c)");
3699 hist->GetAxis(kDataDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
3701 hist->GetAxis(kDataCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
3703 if (fStoreAdditionalJetInformation) {
3704 hist->GetAxis(kDataJetPt)->SetTitle("p_{T}^{jet} (GeV/c)");
3706 hist->GetAxis(kDataZ)->SetTitle("z = p_{T}^{track} / p_{T}^{jet}");
3708 hist->GetAxis(kDataXi)->SetTitle("#xi = ln(p_{T}^{jet} / p_{T}^{track})");
3711 hist->GetAxis(GetIndexOfChargeAxisData())->SetTitle("Charge (e_{0})");
3713 hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetTitle("TOF PID Info");
3714 // Offset is (TMath::Abs(kNoTOFinfo) + 1), such that bin 1 (= first label) corresponds to kNoTOFinfo (< 0)
3715 hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kNoTOFinfo + (TMath::Abs(kNoTOFinfo) + 1), "No TOF Info");
3716 hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kNoTOFpid + (TMath::Abs(kNoTOFinfo) + 1), "No TOF PID");
3717 hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kTOFpion + (TMath::Abs(kNoTOFinfo) + 1), "#pi");
3718 hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kTOFkaon + (TMath::Abs(kNoTOFinfo) + 1), "K");
3719 hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kTOFproton + (TMath::Abs(kNoTOFinfo) + 1), "p");
3723 //________________________________________________________________________
3724 void AliAnalysisTaskPID::SetUpPtResHist(THnSparse* hist, Double_t* binsPt, Double_t* binsJetPt, Double_t* binsCent) const
3726 // Sets bin limits for axes which are not standard binned and the axes titles.
3728 hist->SetBinEdges(kPtResJetPt, binsJetPt);
3729 hist->SetBinEdges(kPtResGenPt, binsPt);
3730 hist->SetBinEdges(kPtResRecPt, binsPt);
3731 hist->SetBinEdges(kPtResCentrality, binsCent);
3734 hist->GetAxis(kPtResJetPt)->SetTitle("p_{T}^{jet, rec} (GeV/c)");
3735 hist->GetAxis(kPtResGenPt)->SetTitle("p_{T}^{gen} (GeV/c)");
3736 hist->GetAxis(kPtResRecPt)->SetTitle("p_{T}^{rec} (GeV/c)");
3738 hist->GetAxis(kPtResCharge)->SetTitle("Charge (e_{0})");
3739 hist->GetAxis(kPtResCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
3743 //________________________________________________________________________
3744 void AliAnalysisTaskPID::SetUpSharedClsHist(THnSparse* hist, Double_t* binsPt, Double_t* binsJetPt) const
3746 // Sets bin limits for axes which are not standard binned and the axes titles.
3748 hist->SetBinEdges(kQASharedClsJetPt, binsJetPt);
3749 hist->SetBinEdges(kQASharedClsPt, binsPt);
3752 hist->GetAxis(kQASharedClsJetPt)->SetTitle("#it{p}_{T}^{jet} (GeV/#it{c})");
3753 hist->GetAxis(kQASharedClsPt)->SetTitle("#it{p}_{T} (GeV/#it{c})");
3754 hist->GetAxis(kQASharedClsNumSharedCls)->SetTitle("#it{N}_{shared}^{cls}");
3755 hist->GetAxis(kQASharedClsPadRow)->SetTitle("Pad row");
3760 //________________________________________________________________________
3761 void AliAnalysisTaskPID::SetUpDeDxCheckHist(THnSparse* hist, const Double_t* binsPt, const Double_t* binsJetPt, const Double_t* binsEtaAbs) const
3763 // Sets bin limits for axes which are not standard binned and the axes titles.
3764 hist->SetBinEdges(kDeDxCheckP, binsPt);
3765 hist->SetBinEdges(kDeDxCheckJetPt, binsJetPt);
3766 hist->SetBinEdges(kDeDxCheckEtaAbs, binsEtaAbs);
3770 hist->GetAxis(kDeDxCheckPID)->SetTitle("PID");
3771 hist->GetAxis(kDeDxCheckPID)->SetBinLabel(1, "e");
3772 hist->GetAxis(kDeDxCheckPID)->SetBinLabel(2, "K");
3773 hist->GetAxis(kDeDxCheckPID)->SetBinLabel(3, "#pi");
3774 hist->GetAxis(kDeDxCheckPID)->SetBinLabel(4, "p");
3777 hist->GetAxis(kDeDxCheckJetPt)->SetTitle("p_{T}^{jet} (GeV/c)");
3778 hist->GetAxis(kDeDxCheckEtaAbs)->SetTitle("|#eta|");
3779 hist->GetAxis(kDeDxCheckP)->SetTitle("p_{TPC} (GeV/c)");
3780 hist->GetAxis(kDeDxCheckDeDx)->SetTitle("TPC dE/dx (arb. unit)");
3784 //________________________________________________________________________
3785 void AliAnalysisTaskPID::SetUpBinZeroStudyHist(THnSparse* hist, const Double_t* binsCent, const Double_t* binsPt) const
3787 // Sets bin limits for axes which are not standard binned and the axes titles.
3788 hist->SetBinEdges(kBinZeroStudyCentrality, binsCent);
3789 hist->SetBinEdges(kBinZeroStudyGenPt, binsPt);
3792 hist->GetAxis(kBinZeroStudyCentrality)->SetTitle(Form("Centrality Percentile (%s)", GetPPCentralityEstimator().Data()));
3793 hist->GetAxis(kBinZeroStudyGenEta)->SetTitle("#it{#eta}^{gen}");
3794 hist->GetAxis(kBinZeroStudyGenPt)->SetTitle("#it{p}_{T}^{gen} (GeV/#it{c})");