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 "AliESDInputHandler.h"
18 #include "AliInputEventHandler.h"
20 #include "AliVVertex.h"
21 #include "AliAnalysisFilter.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()
51 , fPIDcombined(new AliPIDCombined())
52 , fInputFromOtherTask(kFALSE)
54 , fDoEfficiency(kTRUE)
55 , fDoPtResolution(kTRUE)
56 , fStoreCentralityPercentile(kFALSE)
57 , fStoreAdditionalJetInformation(kFALSE)
58 , fTakeIntoAccountMuons(kFALSE)
62 , fTPCDefaultPriors(kFALSE)
63 , fUseMCidForGeneration(kTRUE)
64 , fUseConvolutedGaus(kFALSE)
65 , fkConvolutedGausNPar(3)
66 , fAccuracyNonGaussianTail(1e-8)
67 , fkDeltaPrimeLowLimit(0.02)
68 , fkDeltaPrimeUpLimit(40.0)
69 , fConvolutedGausDeltaPrime(0x0)
73 , fDoAnySystematicStudiesOnTheExpectedSignal(kFALSE)
74 , fSystematicScalingSplinesThreshold(50.)
75 , fSystematicScalingSplinesBelowThreshold(1.0)
76 , fSystematicScalingSplinesAboveThreshold(1.0)
77 , fSystematicScalingEtaCorrectionMomentumThr(0.35)
78 , fSystematicScalingEtaCorrectionLowMomenta(1.0)
79 , fSystematicScalingEtaCorrectionHighMomenta(1.0)
80 , fSystematicScalingEtaSigmaParaThreshold(250.)
81 , fSystematicScalingEtaSigmaParaBelowThreshold(1.0)
82 , fSystematicScalingEtaSigmaParaAboveThreshold(1.0)
83 , fSystematicScalingMultCorrection(1.0)
84 , fCentralityEstimator("V0M")
91 , fGenRespElDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
92 , fGenRespElDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
93 , fGenRespElDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
94 , fGenRespElDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
95 , fGenRespKaDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
96 , fGenRespKaDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
97 , fGenRespKaDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
98 , fGenRespKaDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
99 , fGenRespPiDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
100 , fGenRespPiDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
101 , fGenRespPiDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
102 , fGenRespPiDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
103 , fGenRespMuDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
104 , fGenRespMuDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
105 , fGenRespMuDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
106 , fGenRespMuDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
107 , fGenRespPrDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
108 , fGenRespPrDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
109 , fGenRespPrDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
110 , fGenRespPrDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
112 , fGenRespElDeltaEl(new Double_t[fgkMaxNumGenEntries])
113 , fGenRespElDeltaKa(new Double_t[fgkMaxNumGenEntries])
114 , fGenRespElDeltaPi(new Double_t[fgkMaxNumGenEntries])
115 , fGenRespElDeltaPr(new Double_t[fgkMaxNumGenEntries])
116 , fGenRespKaDeltaEl(new Double_t[fgkMaxNumGenEntries])
117 , fGenRespKaDeltaKa(new Double_t[fgkMaxNumGenEntries])
118 , fGenRespKaDeltaPi(new Double_t[fgkMaxNumGenEntries])
119 , fGenRespKaDeltaPr(new Double_t[fgkMaxNumGenEntries])
120 , fGenRespPiDeltaEl(new Double_t[fgkMaxNumGenEntries])
121 , fGenRespPiDeltaKa(new Double_t[fgkMaxNumGenEntries])
122 , fGenRespPiDeltaPi(new Double_t[fgkMaxNumGenEntries])
123 , fGenRespPiDeltaPr(new Double_t[fgkMaxNumGenEntries])
124 , fGenRespMuDeltaEl(new Double_t[fgkMaxNumGenEntries])
125 , fGenRespMuDeltaKa(new Double_t[fgkMaxNumGenEntries])
126 , fGenRespMuDeltaPi(new Double_t[fgkMaxNumGenEntries])
127 , fGenRespMuDeltaPr(new Double_t[fgkMaxNumGenEntries])
128 , fGenRespPrDeltaEl(new Double_t[fgkMaxNumGenEntries])
129 , fGenRespPrDeltaKa(new Double_t[fgkMaxNumGenEntries])
130 , fGenRespPrDeltaPi(new Double_t[fgkMaxNumGenEntries])
131 , fGenRespPrDeltaPr(new Double_t[fgkMaxNumGenEntries])
133 , fhEventsProcessed(0x0)
134 , fhEventsTriggerSel(0x0)
135 , fhEventsTriggerSelVtxCut(0x0)
136 , fhSkippedTracksForSignalGeneration(0x0)
137 , fhMCgeneratedYieldsPrimaries(0x0)
143 , fOutputContainer(0x0)
144 , fPtResolutionContainer(0x0)
146 // default Constructor
148 AliLog::SetClassDebugLevel("AliAnalysisTaskPID", AliLog::kInfo);
150 fConvolutedGausDeltaPrime = new TF1("convolutedGausDeltaPrime", this, &AliAnalysisTaskPID::ConvolutedGaus,
151 fkDeltaPrimeLowLimit, fkDeltaPrimeUpLimit,
152 fkConvolutedGausNPar, "AliAnalysisTaskPID", "ConvolutedGaus");
154 // Set some arbitrary parameteres, such that the function call will not crash
155 // (although it should not be called with these parameters...)
156 fConvolutedGausDeltaPrime->SetParameter(0, 0);
157 fConvolutedGausDeltaPrime->SetParameter(1, 1);
158 fConvolutedGausDeltaPrime->SetParameter(2, 2);
161 // Initialisation of translation parameters is time consuming.
162 // Therefore, default values will only be initialised if they are really needed.
163 // Otherwise, it is left to the user to set the parameter properly.
164 fConvolutedGaussTransitionPars[0] = -999;
165 fConvolutedGaussTransitionPars[1] = -999;
166 fConvolutedGaussTransitionPars[2] = -999;
169 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
170 fFractionHists[i] = 0x0;
171 fFractionSysErrorHists[i] = 0x0;
173 fPtResolution[i] = 0x0;
177 //________________________________________________________________________
178 AliAnalysisTaskPID::AliAnalysisTaskPID(const char *name)
179 : AliAnalysisTaskPIDV0base(name)
180 , fPIDcombined(new AliPIDCombined())
181 , fInputFromOtherTask(kFALSE)
183 , fDoEfficiency(kTRUE)
184 , fDoPtResolution(kTRUE)
185 , fStoreCentralityPercentile(kFALSE)
186 , fStoreAdditionalJetInformation(kFALSE)
187 , fTakeIntoAccountMuons(kFALSE)
191 , fTPCDefaultPriors(kFALSE)
192 , fUseMCidForGeneration(kTRUE)
193 , fUseConvolutedGaus(kFALSE)
194 , fkConvolutedGausNPar(3)
195 , fAccuracyNonGaussianTail(1e-8)
196 , fkDeltaPrimeLowLimit(0.02)
197 , fkDeltaPrimeUpLimit(40.0)
198 , fConvolutedGausDeltaPrime(0x0)
202 , fDoAnySystematicStudiesOnTheExpectedSignal(kFALSE)
203 , fSystematicScalingSplinesThreshold(50.)
204 , fSystematicScalingSplinesBelowThreshold(1.0)
205 , fSystematicScalingSplinesAboveThreshold(1.0)
206 , fSystematicScalingEtaCorrectionMomentumThr(0.35)
207 , fSystematicScalingEtaCorrectionLowMomenta(1.0)
208 , fSystematicScalingEtaCorrectionHighMomenta(1.0)
209 , fSystematicScalingEtaSigmaParaThreshold(250.)
210 , fSystematicScalingEtaSigmaParaBelowThreshold(1.0)
211 , fSystematicScalingEtaSigmaParaAboveThreshold(1.0)
212 , fSystematicScalingMultCorrection(1.0)
213 , fCentralityEstimator("V0M")
220 , fGenRespElDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
221 , fGenRespElDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
222 , fGenRespElDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
223 , fGenRespElDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
224 , fGenRespKaDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
225 , fGenRespKaDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
226 , fGenRespKaDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
227 , fGenRespKaDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
228 , fGenRespPiDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
229 , fGenRespPiDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
230 , fGenRespPiDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
231 , fGenRespPiDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
232 , fGenRespMuDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
233 , fGenRespMuDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
234 , fGenRespMuDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
235 , fGenRespMuDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
236 , fGenRespPrDeltaPrimeEl(new Double_t[fgkMaxNumGenEntries])
237 , fGenRespPrDeltaPrimeKa(new Double_t[fgkMaxNumGenEntries])
238 , fGenRespPrDeltaPrimePi(new Double_t[fgkMaxNumGenEntries])
239 , fGenRespPrDeltaPrimePr(new Double_t[fgkMaxNumGenEntries])
241 , fGenRespElDeltaEl(new Double_t[fgkMaxNumGenEntries])
242 , fGenRespElDeltaKa(new Double_t[fgkMaxNumGenEntries])
243 , fGenRespElDeltaPi(new Double_t[fgkMaxNumGenEntries])
244 , fGenRespElDeltaPr(new Double_t[fgkMaxNumGenEntries])
245 , fGenRespKaDeltaEl(new Double_t[fgkMaxNumGenEntries])
246 , fGenRespKaDeltaKa(new Double_t[fgkMaxNumGenEntries])
247 , fGenRespKaDeltaPi(new Double_t[fgkMaxNumGenEntries])
248 , fGenRespKaDeltaPr(new Double_t[fgkMaxNumGenEntries])
249 , fGenRespPiDeltaEl(new Double_t[fgkMaxNumGenEntries])
250 , fGenRespPiDeltaKa(new Double_t[fgkMaxNumGenEntries])
251 , fGenRespPiDeltaPi(new Double_t[fgkMaxNumGenEntries])
252 , fGenRespPiDeltaPr(new Double_t[fgkMaxNumGenEntries])
253 , fGenRespMuDeltaEl(new Double_t[fgkMaxNumGenEntries])
254 , fGenRespMuDeltaKa(new Double_t[fgkMaxNumGenEntries])
255 , fGenRespMuDeltaPi(new Double_t[fgkMaxNumGenEntries])
256 , fGenRespMuDeltaPr(new Double_t[fgkMaxNumGenEntries])
257 , fGenRespPrDeltaEl(new Double_t[fgkMaxNumGenEntries])
258 , fGenRespPrDeltaKa(new Double_t[fgkMaxNumGenEntries])
259 , fGenRespPrDeltaPi(new Double_t[fgkMaxNumGenEntries])
260 , fGenRespPrDeltaPr(new Double_t[fgkMaxNumGenEntries])
262 , fhEventsProcessed(0x0)
263 , fhEventsTriggerSel(0x0)
264 , fhEventsTriggerSelVtxCut(0x0)
265 , fhSkippedTracksForSignalGeneration(0x0)
266 , fhMCgeneratedYieldsPrimaries(0x0)
272 , fOutputContainer(0x0)
273 , fPtResolutionContainer(0x0)
277 AliLog::SetClassDebugLevel("AliAnalysisTaskPID", AliLog::kInfo);
279 fConvolutedGausDeltaPrime = new TF1("convolutedGausDeltaPrime", this, &AliAnalysisTaskPID::ConvolutedGaus,
280 fkDeltaPrimeLowLimit, fkDeltaPrimeUpLimit,
281 fkConvolutedGausNPar, "AliAnalysisTaskPID", "ConvolutedGaus");
283 // Set some arbitrary parameteres, such that the function call will not crash
284 // (although it should not be called with these parameters...)
285 fConvolutedGausDeltaPrime->SetParameter(0, 0);
286 fConvolutedGausDeltaPrime->SetParameter(1, 1);
287 fConvolutedGausDeltaPrime->SetParameter(2, 2);
290 // Initialisation of translation parameters is time consuming.
291 // Therefore, default values will only be initialised if they are really needed.
292 // Otherwise, it is left to the user to set the parameter properly.
293 fConvolutedGaussTransitionPars[0] = -999;
294 fConvolutedGaussTransitionPars[1] = -999;
295 fConvolutedGaussTransitionPars[2] = -999;
298 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
299 fFractionHists[i] = 0x0;
300 fFractionSysErrorHists[i] = 0x0;
302 fPtResolution[i] = 0x0;
305 // Define input and output slots here
306 // Input slot #0 works with a TChain
307 DefineInput(0, TChain::Class());
309 DefineOutput(1, TObjArray::Class());
311 DefineOutput(2, AliCFContainer::Class());
313 DefineOutput(3, TObjArray::Class());
317 //________________________________________________________________________
318 AliAnalysisTaskPID::~AliAnalysisTaskPID()
322 CleanupParticleFractionHistos();
324 delete fOutputContainer;
325 fOutputContainer = 0x0;
327 delete fPtResolutionContainer;
328 fPtResolutionContainer = 0x0;
330 delete fConvolutedGausDeltaPrime;
331 fConvolutedGausDeltaPrime = 0x0;
333 delete [] fGenRespElDeltaPrimeEl;
334 delete [] fGenRespElDeltaPrimeKa;
335 delete [] fGenRespElDeltaPrimePi;
336 delete [] fGenRespElDeltaPrimePr;
338 fGenRespElDeltaPrimeEl = 0x0;
339 fGenRespElDeltaPrimeKa = 0x0;
340 fGenRespElDeltaPrimePi = 0x0;
341 fGenRespElDeltaPrimePr = 0x0;
343 delete [] fGenRespKaDeltaPrimeEl;
344 delete [] fGenRespKaDeltaPrimeKa;
345 delete [] fGenRespKaDeltaPrimePi;
346 delete [] fGenRespKaDeltaPrimePr;
348 fGenRespKaDeltaPrimeEl = 0x0;
349 fGenRespKaDeltaPrimeKa = 0x0;
350 fGenRespKaDeltaPrimePi = 0x0;
351 fGenRespKaDeltaPrimePr = 0x0;
353 delete [] fGenRespPiDeltaPrimeEl;
354 delete [] fGenRespPiDeltaPrimeKa;
355 delete [] fGenRespPiDeltaPrimePi;
356 delete [] fGenRespPiDeltaPrimePr;
358 fGenRespPiDeltaPrimeEl = 0x0;
359 fGenRespPiDeltaPrimeKa = 0x0;
360 fGenRespPiDeltaPrimePi = 0x0;
361 fGenRespPiDeltaPrimePr = 0x0;
363 delete [] fGenRespMuDeltaPrimeEl;
364 delete [] fGenRespMuDeltaPrimeKa;
365 delete [] fGenRespMuDeltaPrimePi;
366 delete [] fGenRespMuDeltaPrimePr;
368 fGenRespMuDeltaPrimeEl = 0x0;
369 fGenRespMuDeltaPrimeKa = 0x0;
370 fGenRespMuDeltaPrimePi = 0x0;
371 fGenRespMuDeltaPrimePr = 0x0;
373 delete [] fGenRespPrDeltaPrimeEl;
374 delete [] fGenRespPrDeltaPrimeKa;
375 delete [] fGenRespPrDeltaPrimePi;
376 delete [] fGenRespPrDeltaPrimePr;
378 fGenRespPrDeltaPrimeEl = 0x0;
379 fGenRespPrDeltaPrimeKa = 0x0;
380 fGenRespPrDeltaPrimePi = 0x0;
381 fGenRespPrDeltaPrimePr = 0x0;
383 /*OLD with deltaSpecies
384 delete [] fGenRespElDeltaEl;
385 delete [] fGenRespElDeltaKa;
386 delete [] fGenRespElDeltaPi;
387 delete [] fGenRespElDeltaPr;
389 fGenRespElDeltaEl = 0x0;
390 fGenRespElDeltaKa = 0x0;
391 fGenRespElDeltaPi = 0x0;
392 fGenRespElDeltaPr = 0x0;
394 delete [] fGenRespKaDeltaEl;
395 delete [] fGenRespKaDeltaKa;
396 delete [] fGenRespKaDeltaPi;
397 delete [] fGenRespKaDeltaPr;
399 fGenRespKaDeltaEl = 0x0;
400 fGenRespKaDeltaKa = 0x0;
401 fGenRespKaDeltaPi = 0x0;
402 fGenRespKaDeltaPr = 0x0;
404 delete [] fGenRespPiDeltaEl;
405 delete [] fGenRespPiDeltaKa;
406 delete [] fGenRespPiDeltaPi;
407 delete [] fGenRespPiDeltaPr;
409 fGenRespPiDeltaEl = 0x0;
410 fGenRespPiDeltaKa = 0x0;
411 fGenRespPiDeltaPi = 0x0;
412 fGenRespPiDeltaPr = 0x0;
414 delete [] fGenRespMuDeltaEl;
415 delete [] fGenRespMuDeltaKa;
416 delete [] fGenRespMuDeltaPi;
417 delete [] fGenRespMuDeltaPr;
419 fGenRespMuDeltaEl = 0x0;
420 fGenRespMuDeltaKa = 0x0;
421 fGenRespMuDeltaPi = 0x0;
422 fGenRespMuDeltaPr = 0x0;
424 delete [] fGenRespPrDeltaEl;
425 delete [] fGenRespPrDeltaKa;
426 delete [] fGenRespPrDeltaPi;
427 delete [] fGenRespPrDeltaPr;
429 fGenRespPrDeltaEl = 0x0;
430 fGenRespPrDeltaKa = 0x0;
431 fGenRespPrDeltaPi = 0x0;
432 fGenRespPrDeltaPr = 0x0;
437 //________________________________________________________________________
438 void AliAnalysisTaskPID::SetUpPIDcombined()
440 // Initialise the PIDcombined object
446 printf("File: %s, Line: %d: SetUpPIDcombined\n", (char*)__FILE__, __LINE__);
449 AliFatal("No PIDcombined object!\n");
453 fPIDcombined->SetSelectedSpecies(AliPID::kSPECIESC);
454 fPIDcombined->SetEnablePriors(fUsePriors);
456 if (fTPCDefaultPriors)
457 fPIDcombined->SetDefaultTPCPriors();
459 //TODO use individual priors...
461 // Change detector mask (TPC,TOF,ITS)
462 Int_t detectorMask = AliPIDResponse::kDetTPC;
464 // Reject mismatch mask - mismatch only relevant for TOF at the moment - other detectors do not use it
465 Int_t rejectMismatchMask = AliPIDResponse::kDetTPC;
469 detectorMask = detectorMask | AliPIDResponse::kDetITS;
470 rejectMismatchMask = rejectMismatchMask | AliPIDResponse::kDetITS;
473 detectorMask = detectorMask | AliPIDResponse::kDetTOF;
474 rejectMismatchMask = rejectMismatchMask | AliPIDResponse::kDetTOF;
477 fPIDcombined->SetDetectorMask(detectorMask);
478 fPIDcombined->SetRejectMismatchMask(rejectMismatchMask);
481 printf("File: %s, Line: %d: SetUpPIDcombined done\n", (char*)__FILE__, __LINE__);
485 //________________________________________________________________________
486 void AliAnalysisTaskPID::UserCreateOutputObjects()
492 printf("File: %s, Line: %d: UserCreateOutputObjects\n", (char*)__FILE__, __LINE__);
497 AliAnalysisManager* man = AliAnalysisManager::GetAnalysisManager();
498 AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
501 AliFatal("Input handler needed");
503 // PID response object
504 fPIDResponse = inputHandler->GetPIDResponse();
506 AliFatal("PIDResponse object was not created");
510 printf("File: %s, Line: %d: UserCreateOutputObjects -> Retrieved PIDresponse object\n", (char*)__FILE__, __LINE__);
515 printf("File: %s, Line: %d: UserCreateOutputObjects -> OpenFile(1) successful\n", (char*)__FILE__, __LINE__);
517 fOutputContainer = new TObjArray(1);
518 fOutputContainer->SetName(GetName());
519 fOutputContainer->SetOwner(kTRUE);
521 const Int_t nPtBins = 68;
522 Double_t binsPt[nPtBins+1] = {0. , 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45,
523 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95,
524 1.0, 1.1 , 1.2, 1.3 , 1.4, 1.5 , 1.6, 1.7 , 1.8, 1.9 ,
525 2.0, 2.2 , 2.4, 2.6 , 2.8, 3.0 , 3.2, 3.4 , 3.6, 3.8 ,
526 4.0, 4.5 , 5.0, 5.5 , 6.0, 6.5 , 7.0, 8.0 , 9.0, 10.0,
527 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0,
528 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 };
530 const Int_t nCentBins = 12;
531 //-1 for pp; 90-100 has huge electro-magnetic impurities
532 Double_t binsCent[nCentBins+1] = {-1, 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 };
534 const Int_t nJetPtBins = 11;
535 Double_t binsJetPt[nJetPtBins+1] = {0, 2, 5, 10, 15, 20, 30, 40, 60, 80, 120, 200};
537 const Int_t nChargeBins = 2;
538 const Double_t binsCharge[nChargeBins+1] = { -1.0 - 1e-4, 0.0, 1.0 + 1e-4 };
540 const Int_t nBinsJets = kDataNumAxes;
541 const Int_t nBinsNoJets = nBinsJets - fgkNumJetAxes;
543 const Int_t nBins = fStoreAdditionalJetInformation ? nBinsJets : nBinsNoJets;
545 // deltaPrimeSpecies binning
546 const Int_t deltaPrimeNBins = 600;
547 Double_t deltaPrimeBins[deltaPrimeNBins + 1];
549 const Double_t fromLow = fkDeltaPrimeLowLimit;
550 const Double_t toHigh = fkDeltaPrimeUpLimit;
551 const Double_t factor = TMath::Power(toHigh/fromLow, 1./deltaPrimeNBins);
553 // Log binning for whole deltaPrime range
554 deltaPrimeBins[0] = fromLow;
555 for (Int_t i = 0 + 1; i <= deltaPrimeNBins; i++) {
556 deltaPrimeBins[i] = factor * deltaPrimeBins[i - 1];
559 const Int_t nMCPIDbins = 5;
560 const Double_t mcPIDmin = 0.;
561 const Double_t mcPIDmax = 5.;
563 const Int_t nSelSpeciesBins = 4;
564 const Double_t selSpeciesMin = 0.;
565 const Double_t selSpeciesMax = 4.;
567 const Int_t nZBins = 20;
568 const Double_t zMin = 0.;
569 const Double_t zMax = 1.;
571 const Int_t nXiBins = 70;
572 const Double_t xiMin = 0.;
573 const Double_t xiMax = 7.;
575 const Int_t nTOFpidInfoBins = kNumTOFpidInfoBins;
576 const Double_t tofPIDinfoMin = kNoTOFinfo;
577 const Double_t tofPIDinfoMax = kNoTOFinfo + kNumTOFpidInfoBins;
579 // MC PID, SelectSpecies, pT, deltaPrimeSpecies, centrality percentile, jet pT, z = track_pT/jet_pT, xi = log(1/z)
580 Int_t binsNoJets[nBinsNoJets] = { nMCPIDbins,
588 Int_t binsJets[nBinsJets] = { nMCPIDbins,
599 Int_t *bins = fStoreAdditionalJetInformation ? &binsJets[0] : &binsNoJets[0];
601 Double_t xminNoJets[nBinsNoJets] = { mcPIDmin,
609 Double_t xminJets[nBinsJets] = { mcPIDmin,
620 Double_t *xmin = fStoreAdditionalJetInformation? &xminJets[0] : &xminNoJets[0];
622 Double_t xmaxNoJets[nBinsNoJets] = { mcPIDmax,
625 deltaPrimeBins[deltaPrimeNBins],
627 binsCharge[nChargeBins],
630 Double_t xmaxJets[nBinsJets] = { mcPIDmax,
633 deltaPrimeBins[deltaPrimeNBins],
635 binsJetPt[nJetPtBins],
638 binsCharge[nChargeBins],
641 Double_t *xmax = fStoreAdditionalJetInformation? &xmaxJets[0] : &xmaxNoJets[0];
643 fConvolutedGausDeltaPrime->SetNpx(deltaPrimeNBins);
646 fhPIDdataAll = new THnSparseD("hPIDdataAll","", nBins, bins, xmin, xmax);
647 SetUpHist(fhPIDdataAll, binsPt, deltaPrimeBins, binsCent, binsJetPt);
648 fOutputContainer->Add(fhPIDdataAll);
651 // Generated histograms (so far, bins are the same as for primary THnSparse)
652 const Int_t nGenBins = fStoreAdditionalJetInformation ? nBinsJets : nBinsNoJets;
653 // MC PID, SelectSpecies, Pt, deltaPrimeSpecies, jet pT, z = track_pT/jet_pT, xi = log(1/z)
655 Int_t *genBins = fStoreAdditionalJetInformation ? &binsJets[0] : &binsNoJets[0];
656 Double_t *genXmin = fStoreAdditionalJetInformation? &xminJets[0] : &xminNoJets[0];
657 Double_t *genXmax = fStoreAdditionalJetInformation? &xmaxJets[0] : &xmaxNoJets[0];
660 fhGenEl = new THnSparseD("hGenEl", "", nGenBins, genBins, genXmin, genXmax);
661 SetUpGenHist(fhGenEl, binsPt, deltaPrimeBins, binsCent, binsJetPt);
662 fOutputContainer->Add(fhGenEl);
664 fhGenKa = new THnSparseD("hGenKa", "", nGenBins, genBins, genXmin, genXmax);
665 SetUpGenHist(fhGenKa, binsPt, deltaPrimeBins, binsCent, binsJetPt);
666 fOutputContainer->Add(fhGenKa);
668 fhGenPi = new THnSparseD("hGenPi", "", nGenBins, genBins, genXmin, genXmax);
669 SetUpGenHist(fhGenPi, binsPt, deltaPrimeBins, binsCent, binsJetPt);
670 fOutputContainer->Add(fhGenPi);
672 if (fTakeIntoAccountMuons) {
673 fhGenMu = new THnSparseD("hGenMu", "", nGenBins, genBins, genXmin, genXmax);
674 SetUpGenHist(fhGenMu, binsPt, deltaPrimeBins, binsCent, binsJetPt);
675 fOutputContainer->Add(fhGenMu);
678 fhGenPr = new THnSparseD("hGenPr", "", nGenBins, genBins, genXmin, genXmax);
679 SetUpGenHist(fhGenPr, binsPt, deltaPrimeBins, binsCent, binsJetPt);
680 fOutputContainer->Add(fhGenPr);
682 fhSkippedTracksForSignalGeneration = new TH2D("fhSkippedTracksForSignalGeneration",
683 "Number of tracks skipped for the signal generation;p_{T}^{gen} (GeV/c);TPC signal N",
684 nPtBins, binsPt, 161, -0.5, 160.5);
685 fhSkippedTracksForSignalGeneration->Sumw2();
686 fOutputContainer->Add(fhSkippedTracksForSignalGeneration);
690 fhEventsProcessed = new TH1D("fhEventsProcessed",
691 "Number of events passing trigger selection, vtx and zvtx cuts;Centrality percentile",
692 nCentBins, binsCent);
693 fhEventsProcessed->Sumw2();
694 fOutputContainer->Add(fhEventsProcessed);
696 fhEventsTriggerSelVtxCut = new TH1D("fhEventsTriggerSelVtxCut",
697 "Number of events passing trigger selection and vtx cut;Centrality percentile",
698 nCentBins, binsCent);
699 fhEventsTriggerSelVtxCut->Sumw2();
700 fOutputContainer->Add(fhEventsTriggerSelVtxCut);
702 fhEventsTriggerSel = new TH1D("fhEventsTriggerSel",
703 "Number of events passing trigger selection;Centrality percentile",
704 nCentBins, binsCent);
705 fOutputContainer->Add(fhEventsTriggerSel);
706 fhEventsTriggerSel->Sumw2();
709 // Generated yields within acceptance
710 const Int_t nBinsGenYields = fStoreAdditionalJetInformation ? kGenYieldNumAxes : kGenYieldNumAxes - 3;
711 Int_t genYieldsBins[kGenYieldNumAxes] = { nMCPIDbins, nPtBins, nCentBins, nJetPtBins, nZBins, nXiBins,
713 genYieldsBins[GetIndexOfChargeAxisGenYield()] = nChargeBins;
714 Double_t genYieldsXmin[kGenYieldNumAxes] = { mcPIDmin, binsPt[0], binsCent[0], binsJetPt[0], zMin, xiMin,
716 genYieldsXmin[GetIndexOfChargeAxisGenYield()] = binsCharge[0];
717 Double_t genYieldsXmax[kGenYieldNumAxes] = { mcPIDmax, binsPt[nPtBins], binsCent[nCentBins], binsJetPt[nJetPtBins], zMax, xiMax,
718 binsCharge[nChargeBins] };
719 genYieldsXmax[GetIndexOfChargeAxisGenYield()] = binsCharge[nChargeBins];
722 fhMCgeneratedYieldsPrimaries = new THnSparseD("fhMCgeneratedYieldsPrimaries",
723 "Generated yields w/o reco and cuts inside acceptance (physical primaries)",
724 nBinsGenYields, genYieldsBins, genYieldsXmin, genYieldsXmax);
725 SetUpGenYieldHist(fhMCgeneratedYieldsPrimaries, binsPt, binsCent, binsJetPt);
726 fOutputContainer->Add(fhMCgeneratedYieldsPrimaries);
729 // Container with several process steps (generated and reconstructed level with some variations)
734 printf("File: %s, Line: %d: UserCreateOutputObjects -> OpenFile(2) successful\n", (char*)__FILE__, __LINE__);
736 // Array for the number of bins in each dimension
737 // Dimensions: MC-ID, trackPt, trackEta, trackCharge, cenrality percentile, jetPt, z, xi TODO phi???
738 const Int_t nEffDims = fStoreAdditionalJetInformation ? kEffNumAxes : kEffNumAxes - 3; // Number of dimensions for the efficiency
740 const Int_t nMCIDbins = AliPID::kSPECIES;
741 Double_t binsMCID[nMCIDbins + 1];
743 for(Int_t i = 0; i <= nMCIDbins; i++) {
747 const Int_t nEtaBins = 18;
748 const Double_t binsEta[nEtaBins+1] = {-0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1,
749 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 };
751 const Int_t nEffBins[kEffNumAxes] = { nMCIDbins, nPtBins, nEtaBins, nChargeBins, nCentBins, nJetPtBins, nZBins, nXiBins };
753 fContainerEff = new AliCFContainer("containerEff", "Reconstruction Efficiency x Acceptance x Resolution and Secondary Correction",
754 kNumSteps, nEffDims, nEffBins);
756 // Setting the bin limits
757 fContainerEff->SetBinLimits(kEffMCID, binsMCID);
758 fContainerEff->SetBinLimits(kEffTrackPt, binsPt);
759 fContainerEff->SetBinLimits(kEffTrackEta, binsEta);
760 fContainerEff->SetBinLimits(kEffTrackCharge, binsCharge);
761 fContainerEff->SetBinLimits(kEffCentrality, binsCent);
762 if (fStoreAdditionalJetInformation) {
763 fContainerEff->SetBinLimits(kEffJetPt, binsJetPt);
764 fContainerEff->SetBinLimits(kEffZ, zMin, zMax);
765 fContainerEff->SetBinLimits(kEffXi, xiMin, xiMax);
768 fContainerEff->SetVarTitle(kEffMCID,"MC ID");
769 fContainerEff->SetVarTitle(kEffTrackPt,"p_{T} (GeV/c)");
770 fContainerEff->SetVarTitle(kEffTrackEta,"#eta");
771 fContainerEff->SetVarTitle(kEffTrackCharge,"Charge (e_{0})");
772 fContainerEff->SetVarTitle(kEffCentrality, "Centrality Percentile");
773 if (fStoreAdditionalJetInformation) {
774 fContainerEff->SetVarTitle(kEffJetPt, "p_{T}^{jet} (GeV/c)");
775 fContainerEff->SetVarTitle(kEffZ, "z = p_{T}^{track} / p_{T}^{jet}");
776 fContainerEff->SetVarTitle(kEffXi, "#xi = ln(p_{T}^{jet} / p_{T}^{track})");
779 // Define clean MC sample
780 fContainerEff->SetStepTitle(kStepGenWithGenCuts, "Particle level, cuts on particle level");
781 // For Acceptance x Efficiency correction of primaries
782 fContainerEff->SetStepTitle(kStepRecWithGenCuts, "Detector level (rec) with cuts on particle level");
783 // For (pT) resolution correction
784 fContainerEff->SetStepTitle(kStepRecWithGenCutsMeasuredObs,
785 "Detector level (rec) with cuts on particle level with measured observables");
786 // For secondary correction
787 fContainerEff->SetStepTitle(kStepRecWithRecCutsMeasuredObs,
788 "Detector level, all cuts on detector level with measured observables");
789 fContainerEff->SetStepTitle(kStepRecWithRecCutsPrimaries,
790 "Detector level, all cuts on detector level, only MC primaries");
791 fContainerEff->SetStepTitle(kStepRecWithRecCutsMeasuredObsPrimaries,
792 "Detector level, all cuts on detector level with measured observables, only MC primaries");
793 fContainerEff->SetStepTitle(kStepRecWithRecCutsMeasuredObsStrangenessScaled,
794 "Detector level (strangeness scaled), all cuts on detector level with measured observables");
797 if (fDoPID || fDoEfficiency) {
799 fh2FFJetPtRec = new TH2D("fh2FFJetPtRec", "Number of reconstructed jets;Centrality Percentile;p_{T}^{jet} (GeV/c)",
800 nCentBins, binsCent, nJetPtBins, binsJetPt);
801 fh2FFJetPtRec->Sumw2();
802 fOutputContainer->Add(fh2FFJetPtRec);
803 fh2FFJetPtGen = new TH2D("fh2FFJetPtGen", "Number of generated jets;Centrality Percentile;p_{T}^{jet} (GeV/c)",
804 nCentBins, binsCent, nJetPtBins, binsJetPt);
805 fh2FFJetPtGen->Sumw2();
806 fOutputContainer->Add(fh2FFJetPtGen);
809 // Pythia information
810 fh1Xsec = new TProfile("fh1Xsec", "xsec from pyxsec.root", 1, 0, 1);
812 fh1Xsec->GetXaxis()->SetBinLabel(1, "<#sigma>");
813 fh1Trials = new TH1D("fh1Trials", "trials from pyxsec.root", 1, 0, 1);
815 fh1Trials->GetXaxis()->SetBinLabel(1, "#sum{ntrials}");
817 fOutputContainer->Add(fh1Xsec);
818 fOutputContainer->Add(fh1Trials);
820 if (fDoPtResolution) {
824 printf("File: %s, Line: %d: UserCreateOutputObjects -> OpenFile(3) successful\n", (char*)__FILE__, __LINE__);
826 fPtResolutionContainer = new TObjArray(1);
827 fPtResolutionContainer->SetName(Form("%s_PtResolution", GetName()));
828 fPtResolutionContainer->SetOwner(kTRUE);
830 const Int_t nPtBinsRes = 100;
831 Double_t pTbinsRes[nPtBinsRes + 1];
833 const Double_t fromLowPtRes = 0.15;
834 const Double_t toHighPtRes = 50.;
835 const Double_t factorPtRes = TMath::Power(toHighPtRes/fromLowPtRes, 1./nPtBinsRes);
836 // Log binning for whole pT range
837 pTbinsRes[0] = fromLowPtRes;
838 for (Int_t i = 0 + 1; i <= nPtBinsRes; i++) {
839 pTbinsRes[i] = factorPtRes * pTbinsRes[i - 1];
842 const Int_t nBinsPtResolution = kPtResNumAxes;
843 Int_t ptResolutionBins[kPtResNumAxes] = { nJetPtBins, nPtBinsRes, nPtBinsRes,
844 nChargeBins, nCentBins };
845 Double_t ptResolutionXmin[kPtResNumAxes] = { binsJetPt[0], pTbinsRes[0], pTbinsRes[0],
846 binsCharge[0], binsCent[0] };
847 Double_t ptResolutionXmax[kPtResNumAxes] = { binsJetPt[nJetPtBins], pTbinsRes[nPtBinsRes], pTbinsRes[nPtBinsRes],
848 binsCharge[nChargeBins], binsCent[nCentBins] };
850 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
851 fPtResolution[i] = new THnSparseD(Form("fPtResolution_%s", AliPID::ParticleShortName(i)),
852 Form("Pt resolution for primaries, %s", AliPID::ParticleLatexName(i)),
853 nBinsPtResolution, ptResolutionBins, ptResolutionXmin, ptResolutionXmax);
854 SetUpPtResHist(fPtResolution[i], pTbinsRes, binsJetPt, binsCent);
855 fPtResolutionContainer->Add(fPtResolution[i]);
860 printf("File: %s, Line: %d: UserCreateOutputObjects -> Posting output data\n", (char*)__FILE__, __LINE__);
862 PostData(1, fOutputContainer);
863 PostData(2, fContainerEff);
864 PostData(3, fPtResolutionContainer);
867 printf("File: %s, Line: %d: UserCreateOutputObjects -> Done\n", (char*)__FILE__, __LINE__);
871 //________________________________________________________________________
872 void AliAnalysisTaskPID::UserExec(Option_t *)
875 // Called for each event
878 printf("File: %s, Line: %d: UserExec\n", (char*)__FILE__, __LINE__);
880 // No processing of event, if input is fed in directly from another task
881 if (fInputFromOtherTask)
885 printf("File: %s, Line: %d: UserExec -> Processing started\n", (char*)__FILE__, __LINE__);
887 fEvent = dynamic_cast<AliVEvent*>(InputEvent());
889 Printf("ERROR: fEvent not available");
893 fMC = dynamic_cast<AliMCEvent*>(MCEvent());
895 if (!fPIDResponse || !fPIDcombined)
898 Double_t centralityPercentile = -1;
899 if (fStoreCentralityPercentile)
900 centralityPercentile = fEvent->GetCentrality()->GetCentralityPercentile(fCentralityEstimator.Data());
902 IncrementEventCounter(centralityPercentile, kTriggerSel);
904 // Check if vertex is ok, but don't apply cut on z position
905 if (!GetVertexIsOk(fEvent, kFALSE))
908 fESD = dynamic_cast<AliESDEvent*>(fEvent);
909 const AliVVertex* primaryVertex = fESD ? fESD->GetPrimaryVertexTracks() : fEvent->GetPrimaryVertex();
913 if(primaryVertex->GetNContributors() <= 0)
916 IncrementEventCounter(centralityPercentile, kTriggerSelAndVtxCut);
918 // Now check again, but also require z position to be in desired range
919 if (!GetVertexIsOk(fEvent, kTRUE))
922 IncrementEventCounter(centralityPercentile, kTriggerSelAndVtxCutAndZvtxCut);
924 Double_t magField = fEvent->GetMagneticField();
927 if (fDoPID || fDoEfficiency) {
928 for (Int_t iPart = 0; iPart < fMC->GetNumberOfTracks(); iPart++) {
929 AliMCParticle *mcPart = dynamic_cast<AliMCParticle*>(fMC->GetTrack(iPart));
934 // Define clean MC sample with corresponding particle level track cuts:
935 // - MC-track must be in desired eta range
936 // - MC-track must be physical primary
937 // - Species must be one of those in question (everything else goes to the overflow bin of mcID)
939 // Geometrie should be the same as on the reconstructed level -> By definition analysis within this eta interval
940 if (!IsInAcceptedEtaRange(TMath::Abs(mcPart->Eta()))) continue;
942 Int_t mcID = PDGtoMCID(mcPart->PdgCode());
944 // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
945 Double_t chargeMC = mcPart->Charge() / 3.;
947 if (TMath::Abs(chargeMC) < 0.01)
948 continue; // Reject neutral particles (only relevant, if mcID is not used)
950 if (!fMC->IsPhysicalPrimary(iPart))
954 Double_t valuesGenYield[kGenYieldNumAxes] = { mcID, mcPart->Pt(), centralityPercentile, -1, -1, -1, -1 };
955 valuesGenYield[GetIndexOfChargeAxisGenYield()] = chargeMC;
957 fhMCgeneratedYieldsPrimaries->Fill(valuesGenYield);
962 Double_t valueEff[kEffNumAxes] = { mcID, mcPart->Pt(), mcPart->Eta(), chargeMC, centralityPercentile,
965 fContainerEff->Fill(valueEff, kStepGenWithGenCuts);
971 // Track loop to fill a Train spectrum
972 for (Int_t iTracks = 0; iTracks < fEvent->GetNumberOfTracks(); iTracks++) {
973 AliVTrack* track = dynamic_cast<AliVTrack*>(fEvent->GetTrack(iTracks));
975 Printf("ERROR: Could not retrieve track %d", iTracks);
980 // Apply detector level track cuts
981 Double_t dEdxTPC = fPIDResponse->IsTunedOnData() ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
985 if(fTrackFilter && !fTrackFilter->IsSelected(track))
988 if (GetUseTPCCutMIGeo()) {
989 if (!TPCCutMIGeo(track, fEvent))
992 else if (GetUseTPCnclCut()) {
993 if (!TPCnclCut(track))
998 if (!PhiPrimeCut(track, magField))
999 continue; // reject track
1002 Int_t pdg = 0; // = 0 indicates data for the moment
1003 AliMCParticle* mcTrack = 0x0;
1004 Int_t mcID = AliPID::kUnknown;
1008 label = track->GetLabel();
1013 mcTrack = dynamic_cast<AliMCParticle*>(fMC->GetTrack(TMath::Abs(label)));
1015 Printf("ERROR: Could not retrieve mcTrack with label %d for track %d", label, iTracks);
1019 pdg = mcTrack->PdgCode();
1020 mcID = PDGtoMCID(mcTrack->PdgCode());
1022 if (fDoEfficiency) {
1023 // For efficiency: Reconstructed track has survived all cuts on the detector level (excluding acceptance)
1024 // and has an associated MC track which is a physical primary and was generated inside the acceptance
1025 if (fMC->IsPhysicalPrimary(TMath::Abs(label)) &&
1026 IsInAcceptedEtaRange(TMath::Abs(mcTrack->Eta()))) {
1028 // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
1029 Double_t value[kEffNumAxes] = { mcID, mcTrack->Pt(), mcTrack->Eta(), mcTrack->Charge() / 3., centralityPercentile,
1031 fContainerEff->Fill(value, kStepRecWithGenCuts);
1033 Double_t valueMeas[kEffNumAxes] = { mcID, track->Pt(), track->Eta(), track->Charge(), centralityPercentile,
1035 fContainerEff->Fill(valueMeas, kStepRecWithGenCutsMeasuredObs);
1040 // Only process tracks inside the desired eta window
1041 if (!IsInAcceptedEtaRange(TMath::Abs(track->Eta()))) continue;
1044 ProcessTrack(track, pdg, centralityPercentile, -1); // No jet information in this case -> Set jet pT to -1
1046 if (fDoPtResolution) {
1047 if (mcTrack && fMC->IsPhysicalPrimary(TMath::Abs(label))) {
1048 // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
1049 Double_t valuePtRes[kPtResNumAxes] = { -1, mcTrack->Pt(), track->Pt(), mcTrack->Charge() / 3., centralityPercentile };
1050 FillPtResolution(mcID, valuePtRes);
1054 if (fDoEfficiency) {
1056 Double_t valueRecAllCuts[kEffNumAxes] = { mcID, track->Pt(), track->Eta(), track->Charge(), centralityPercentile,
1058 fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObs);
1060 Double_t weight = IsSecondaryWithStrangeMotherMC(fMC, TMath::Abs(label)) ?
1061 GetMCStrangenessFactorCMS(fMC, mcTrack) : 1.0;
1062 fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObsStrangenessScaled, weight);
1064 // AliMCParticle->Charge() calls TParticlePDG->Charge(), which returns the charge in units of e0 / 3
1065 Double_t valueGenAllCuts[kEffNumAxes] = { mcID, mcTrack->Pt(), mcTrack->Eta(), mcTrack->Charge() / 3.,
1066 centralityPercentile, -1, -1, -1 };
1067 if (fMC->IsPhysicalPrimary(TMath::Abs(label))) {
1068 fContainerEff->Fill(valueRecAllCuts, kStepRecWithRecCutsMeasuredObsPrimaries);
1069 fContainerEff->Fill(valueGenAllCuts, kStepRecWithRecCutsPrimaries);
1076 printf("File: %s, Line: %d: UserExec -> Processing done\n", (char*)__FILE__, __LINE__);
1081 printf("File: %s, Line: %d: UserExec -> Done\n", (char*)__FILE__, __LINE__);
1084 //________________________________________________________________________
1085 void AliAnalysisTaskPID::Terminate(const Option_t *)
1087 // Draw result to the screen
1088 // Called once at the end of the query
1092 //_____________________________________________________________________________
1093 void AliAnalysisTaskPID::CheckDoAnyStematicStudiesOnTheExpectedSignal()
1095 // Check whether at least one scale factor indicates the ussage of systematic studies
1096 // and set internal flag accordingly.
1098 fDoAnySystematicStudiesOnTheExpectedSignal = kFALSE;
1101 if ((TMath::Abs(fSystematicScalingSplinesBelowThreshold - 1.0) > fgkEpsilon) ||
1102 (TMath::Abs(fSystematicScalingSplinesAboveThreshold - 1.0) > fgkEpsilon)) {
1103 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1107 if ((TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon) ||
1108 (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon)) {
1109 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1113 if ((TMath::Abs(fSystematicScalingEtaSigmaParaBelowThreshold - 1.0) > fgkEpsilon) ||
1114 (TMath::Abs(fSystematicScalingEtaSigmaParaAboveThreshold - 1.0) > fgkEpsilon)) {
1115 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1119 if (TMath::Abs(fSystematicScalingMultCorrection - 1.0) > fgkEpsilon) {
1120 fDoAnySystematicStudiesOnTheExpectedSignal = kTRUE;
1126 //_____________________________________________________________________________
1127 Int_t AliAnalysisTaskPID::PDGtoMCID(Int_t pdg)
1129 // Returns the corresponding AliPID index to the given pdg code.
1130 // Returns AliPID::kUnkown if pdg belongs to a not considered species.
1132 Int_t absPDGcode = TMath::Abs(pdg);
1133 if (absPDGcode == 211) {//Pion
1134 return AliPID::kPion;
1136 else if (absPDGcode == 321) {//Kaon
1137 return AliPID::kKaon;
1139 else if (absPDGcode == 2212) {//Proton
1140 return AliPID::kProton;
1142 else if (absPDGcode == 11) {//Electron
1143 return AliPID::kElectron;
1145 else if (absPDGcode == 13) {//Muon
1146 return AliPID::kMuon;
1149 return AliPID::kUnknown;
1153 //_____________________________________________________________________________
1154 void AliAnalysisTaskPID::GetJetTrackObservables(Double_t trackPt, Double_t jetPt, Double_t& z, Double_t& xi)
1156 // Uses trackPt and jetPt to obtain z and xi.
1158 z = (jetPt > 0 && trackPt >= 0) ? (trackPt / jetPt) : -1;
1159 xi = (z > 0) ? TMath::Log(1. / z) : -1;
1161 if(trackPt > (1. - 1e-06) * jetPt && trackPt < (1. + 1e-06) * jetPt) { // case z=1 : move entry to last histo bin <1
1168 //_____________________________________________________________________________
1169 void AliAnalysisTaskPID::CleanupParticleFractionHistos()
1171 // Delete histos with particle fractions
1173 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1174 delete fFractionHists[i];
1175 fFractionHists[i] = 0x0;
1177 delete fFractionSysErrorHists[i];
1178 fFractionSysErrorHists[i] = 0x0;
1183 //_____________________________________________________________________________
1184 Double_t AliAnalysisTaskPID::ConvolutedGaus(const Double_t* xx, const Double_t* par) const
1186 // Convolutes gauss with an exponential tail which describes dEdx-response better than pure gaussian
1188 const Double_t mean = par[0];
1189 const Double_t sigma = par[1];
1190 const Double_t lambda = par[2];
1193 printf("File: %s, Line: %d: ConvolutedGaus: mean %e, sigma %e, lambda %e\n", (char*)__FILE__, __LINE__, mean, sigma, lambda);
1195 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);
1199 //_____________________________________________________________________________
1200 inline Double_t AliAnalysisTaskPID::FastGaus(Double_t x, Double_t mean, Double_t sigma) const
1202 // Calculate an unnormalised gaussian function with mean and sigma.
1204 if (sigma < fgkEpsilon)
1207 const Double_t arg = (x - mean) / sigma;
1208 return exp(-0.5 * arg * arg);
1212 //_____________________________________________________________________________
1213 inline Double_t AliAnalysisTaskPID::FastNormalisedGaus(Double_t x, Double_t mean, Double_t sigma) const
1215 // Calculate a normalised (divided by sqrt(2*Pi)*sigma) gaussian function with mean and sigma.
1217 if (sigma < fgkEpsilon)
1220 const Double_t arg = (x - mean) / sigma;
1221 const Double_t res = exp(-0.5 * arg * arg);
1222 return res / (2.50662827463100024 * sigma); //sqrt(2*Pi)=2.50662827463100024
1226 //_____________________________________________________________________________
1227 Int_t AliAnalysisTaskPID::FindBinWithinRange(TAxis* axis, Double_t value) const
1229 // Find the corresponding bin of the axis. Values outside the range (also under and overflow) will be set to the first/last
1232 Int_t bin = axis->FindFixBin(value);
1236 if (bin > axis->GetNbins())
1237 bin = axis->GetNbins();
1243 //_____________________________________________________________________________
1244 Int_t AliAnalysisTaskPID::FindFirstBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yBin,
1247 // Kind of projects a TH3 to 1 bin combination in y and z
1248 // and looks for the first x bin above a threshold for this projection.
1249 // If no such bin is found, -1 is returned.
1254 Int_t nBinsX = hist->GetNbinsX();
1255 for (Int_t xBin = 1; xBin <= nBinsX; xBin++) {
1256 if (hist->GetBinContent(xBin, yBin, zBin) > threshold)
1264 //_____________________________________________________________________________
1265 Int_t AliAnalysisTaskPID::FindLastBinAboveIn3dSubset(const TH3* hist, Double_t threshold, Int_t yBin,
1268 // Kind of projects a TH3 to 1 bin combination in y and z
1269 // and looks for the last x bin above a threshold for this projection.
1270 // If no such bin is found, -1 is returned.
1275 Int_t nBinsX = hist->GetNbinsX();
1276 for (Int_t xBin = nBinsX; xBin >= 1; xBin--) {
1277 if (hist->GetBinContent(xBin, yBin, zBin) > threshold)
1285 //_____________________________________________________________________________
1286 Bool_t AliAnalysisTaskPID::GetParticleFraction(Double_t trackPt, Double_t jetPt, Double_t centralityPercentile,
1287 AliPID::EParticleType species,
1288 Double_t& fraction, Double_t& fractionErrorStat, Double_t& fractionErrorSys) const
1290 // Computes the particle fraction for the corresponding species for the given trackPt, jetPt and centrality.
1291 // Use jetPt = -1 for inclusive spectra and centralityPercentile = -1 for pp.
1292 // On success (return value kTRUE), fraction contains the particle fraction, fractionErrorStat(Sys) the sigma of its
1293 // statistical (systematic) error
1296 fractionErrorStat = 999.;
1297 fractionErrorSys = 999.;
1299 if (species > AliPID::kProton || species < AliPID::kElectron) {
1300 AliError(Form("Only fractions for species index %d to %d availabe, but not for the requested one: %d", 0, AliPID::kProton, species));
1304 if (!fFractionHists[species]) {
1305 AliError(Form("Histo with particle fractions for species %d not loaded!", species));
1310 Int_t jetPtBin = FindBinWithinRange(fFractionHists[species]->GetYaxis(), jetPt);
1311 Int_t centBin = FindBinWithinRange(fFractionHists[species]->GetZaxis(), centralityPercentile);
1313 // The following interpolation takes the bin content of the first/last available bin,
1314 // if requested point lies beyond bin center of first/last bin.
1315 // The interpolation is only done for the x-axis (track pT), i.e. jetPtBin and centBin are fix,
1316 // because the analysis will anyhow be run in bins of jetPt and centrality and
1317 // it is not desired to correlate different jetPt bins via interpolation.
1319 // The same procedure is used for the error of the fraction
1320 TAxis* xAxis = fFractionHists[species]->GetXaxis();
1322 // No interpolation to values beyond the centers of the first/last bins (we don't know exactly where the spectra start or stop,
1323 // thus, search for the first and last bin above 0.0 to constrain the range
1324 Int_t firstBin = TMath::Max(1, FindFirstBinAboveIn3dSubset(fFractionHists[species], 0.0, jetPtBin, centBin));
1325 Int_t lastBin = TMath::Min(fFractionHists[species]->GetNbinsX(),
1326 FindLastBinAboveIn3dSubset(fFractionHists[species], 0.0, jetPtBin, centBin));
1328 if (trackPt <= xAxis->GetBinCenter(firstBin)) {
1329 fraction = fFractionHists[species]->GetBinContent(firstBin, jetPtBin, centBin);
1330 fractionErrorStat = fFractionHists[species]->GetBinError(firstBin, jetPtBin, centBin);
1331 fractionErrorSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(firstBin, jetPtBin, centBin) : 0.;
1333 else if (trackPt >= xAxis->GetBinCenter(lastBin)) {
1334 fraction = fFractionHists[species]->GetBinContent(lastBin, jetPtBin, centBin);
1335 fractionErrorStat = fFractionHists[species]->GetBinError(lastBin, jetPtBin, centBin);
1336 fractionErrorSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(lastBin, jetPtBin, centBin) : 0.;
1339 Double_t x0 = 0., x1 = 0., y0 = 0., y1 = 0.;
1340 Double_t y0errStat = 0., y1errStat = 0., y0errSys = 0., y1errSys = 0.;
1341 Int_t trackPtBin = xAxis->FindBin(trackPt);
1343 // Linear interpolation between nearest neighbours in trackPt
1344 if (trackPt <= xAxis->GetBinCenter(trackPtBin)) {
1345 y0 = fFractionHists[species]->GetBinContent(trackPtBin - 1, jetPtBin, centBin);
1346 y0errStat = fFractionHists[species]->GetBinError(trackPtBin - 1, jetPtBin, centBin);
1347 y0errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin - 1, jetPtBin, centBin)
1349 x0 = xAxis->GetBinCenter(trackPtBin - 1);
1350 y1 = fFractionHists[species]->GetBinContent(trackPtBin, jetPtBin, centBin);
1351 y1errStat = fFractionHists[species]->GetBinError(trackPtBin, jetPtBin, centBin);
1352 y1errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin, jetPtBin, centBin)
1354 x1 = xAxis->GetBinCenter(trackPtBin);
1357 y0 = fFractionHists[species]->GetBinContent(trackPtBin, jetPtBin, centBin);
1358 y0errStat = fFractionHists[species]->GetBinError(trackPtBin, jetPtBin, centBin);
1359 y0errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin, jetPtBin, centBin)
1361 x0 = xAxis->GetBinCenter(trackPtBin);
1362 y1 = fFractionHists[species]->GetBinContent(trackPtBin + 1, jetPtBin, centBin);
1363 y1errStat = fFractionHists[species]->GetBinError(trackPtBin + 1, jetPtBin, centBin);
1364 y1errSys = fFractionSysErrorHists[species] ? fFractionSysErrorHists[species]->GetBinError(trackPtBin + 1, jetPtBin, centBin)
1366 x1 = xAxis->GetBinCenter(trackPtBin + 1);
1369 // Per construction: x0 < trackPt < x1
1370 fraction = y0 + (trackPt - x0) * ((y1 - y0) / (x1 - x0));
1371 fractionErrorStat = y0errStat + (trackPt - x0) * ((y1errStat - y0errStat) / (x1 - x0));
1372 fractionErrorSys = fFractionSysErrorHists[species] ? (y0errSys + (trackPt - x0) * ((y1errSys - y0errSys) / (x1 - x0))) : 0.;
1379 //_____________________________________________________________________________
1380 Bool_t AliAnalysisTaskPID::GetParticleFractions(Double_t trackPt, Double_t jetPt, Double_t centralityPercentile,
1381 Double_t* prob, Int_t smearSpeciesByError,
1382 Int_t takeIntoAccountSpeciesSysError, Bool_t uniformSystematicError) const
1384 // Fills the particle fractions for the given trackPt, jetPt and centrality into "prob".
1385 // Use jetPt = -1 for inclusive spectra and centralityPercentile = -1 for pp.
1386 // If smearSpeciesByError is >= 0 && < AliPID::kSPECIES, the returned fractions will be a random number distributed
1387 // with a gauss with mean being the corresponding particle fraction and sigma it's error for the considered species
1388 // "smearSpeciesByError".
1389 // Note that in this case the fractions for all species will NOT sum up to 1!
1390 // Thus, all other species fractions will be re-scaled weighted with their corresponding statistical error.
1391 // A similar procedure is used for "takeIntoAccountSpeciesSysError": The systematic error of the corresponding species
1392 // is used to generate a random number with uniform distribution in [mean - sysError, mean + sysError] for the new mean
1393 // (in cace of uniformSystematicError = kTRUE, otherwise it will be a gaus(mean, sysError)),
1394 // then the other species will be re-scaled according to their systematic errors.
1395 // First, the systematic error uncertainty procedure will be performed (that is including re-scaling), then the statistical
1396 // uncertainty procedure.
1397 // On success, kTRUE is returned.
1399 if (!prob || smearSpeciesByError >= AliPID::kSPECIES || takeIntoAccountSpeciesSysError >= AliPID::kSPECIES)
1402 Double_t probTemp[AliPID::kSPECIES];
1403 Double_t probErrorStat[AliPID::kSPECIES];
1404 Double_t probErrorSys[AliPID::kSPECIES];
1406 Bool_t success = kTRUE;
1407 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kElectron,
1408 probTemp[AliPID::kElectron], probErrorStat[AliPID::kElectron],
1409 probErrorSys[AliPID::kElectron]);
1410 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kMuon,
1411 probTemp[AliPID::kMuon], probErrorStat[AliPID::kMuon], probErrorSys[AliPID::kMuon]);
1412 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kPion,
1413 probTemp[AliPID::kPion], probErrorStat[AliPID::kPion], probErrorSys[AliPID::kPion]);
1414 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kKaon,
1415 probTemp[AliPID::kKaon], probErrorStat[AliPID::kKaon], probErrorSys[AliPID::kKaon]);
1416 success = success && GetParticleFraction(trackPt, jetPt, centralityPercentile, AliPID::kProton,
1417 probTemp[AliPID::kProton], probErrorStat[AliPID::kProton], probErrorSys[AliPID::kProton]);
1422 // If desired, take into account the systematic error of the corresponding species and re-generate probTemp accordingly
1423 if (takeIntoAccountSpeciesSysError >= 0) {
1424 // Generate random fraction of the considered species "smearSpeciesByError" according to mean and sigma
1425 Double_t generatedFraction = uniformSystematicError
1426 ? fRandom->Rndm() * 2. * probErrorSys[takeIntoAccountSpeciesSysError]
1427 - probErrorSys[takeIntoAccountSpeciesSysError]
1428 + probTemp[takeIntoAccountSpeciesSysError]
1429 : fRandom->Gaus(probTemp[takeIntoAccountSpeciesSysError],
1430 probErrorSys[takeIntoAccountSpeciesSysError]);
1432 // Catch cases with invalid fraction (can happen for large errors), i.e. fraction < 0 or > 1
1433 if (generatedFraction < 0.)
1434 generatedFraction = 0.;
1435 else if (generatedFraction > 1.)
1436 generatedFraction = 1.;
1438 // Calculate difference from original fraction (original fractions sum up to 1!)
1439 Double_t deltaFraction = generatedFraction - probTemp[takeIntoAccountSpeciesSysError];
1441 // Fractions must (including errors) lie inside [0,1] -> Adapt weights accordingly by setting the errors
1442 if (deltaFraction > 0) {
1443 // Some part will be SUBTRACTED from the other fractions
1444 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1445 if (probTemp[i] - probErrorSys[i] < 0)
1446 probErrorSys[i] = probTemp[i];
1450 // Some part will be ADDED to the other fractions
1451 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1452 if (probTemp[i] + probErrorSys[i] > 1)
1453 probErrorSys[i] = 1. - probTemp[i];
1457 // Compute summed weight of all fractions except for the considered one
1458 Double_t summedWeight = 0.;
1459 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1460 if (i != takeIntoAccountSpeciesSysError)
1461 summedWeight += probErrorSys[i];
1464 // Compute the weight for the other species
1466 if (summedWeight <= 1e-13) {
1467 // If this happens for some reason (it should not!), just assume flat weight
1468 printf("Error: summedWeight (sys error) ~ 0 for trackPt %f, jetPt %f, centralityPercentile %f. Setting flat weight!\n",
1469 trackPt, jetPt, centralityPercentile);
1472 Double_t weight[AliPID::kSPECIES];
1473 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1474 if (i != takeIntoAccountSpeciesSysError) {
1475 if (summedWeight > 1e-13)
1476 weight[i] = probErrorSys[i] / summedWeight;
1478 weight[i] = probErrorSys[i] / (AliPID::kSPECIES - 1);
1482 // For the final generated fractions, set the generated value for the considered species
1483 // and the generated value minus delta times statistical weight
1484 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1485 if (i != takeIntoAccountSpeciesSysError)
1486 probTemp[i] = probTemp[i] - weight[i] * deltaFraction;
1488 probTemp[i] = generatedFraction;
1492 // Using the values of probTemp (either the original ones or those after taking into account the systematic error),
1493 // calculate the final fractions - if the statistical error is to be taken into account, smear the corresponding
1494 // fraction. If not, just write probTemp to the final result array.
1495 if (smearSpeciesByError >= 0) {
1496 // Generate random fraction of the considered species "smearSpeciesByError" according to mean and sigma
1497 Double_t generatedFraction = fRandom->Gaus(probTemp[smearSpeciesByError], probErrorStat[smearSpeciesByError]);
1499 // Catch cases with invalid fraction (can happen for large errors), i.e. fraction < 0 or > 1
1500 if (generatedFraction < 0.)
1501 generatedFraction = 0.;
1502 else if (generatedFraction > 1.)
1503 generatedFraction = 1.;
1505 // Calculate difference from original fraction (original fractions sum up to 1!)
1506 Double_t deltaFraction = generatedFraction - probTemp[smearSpeciesByError];
1508 // Fractions must (including errors) lie inside [0,1] -> Adapt weights accordingly by setting the errors
1509 if (deltaFraction > 0) {
1510 // Some part will be SUBTRACTED from the other fractions
1511 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1512 if (probTemp[i] - probErrorStat[i] < 0)
1513 probErrorStat[i] = probTemp[i];
1517 // Some part will be ADDED to the other fractions
1518 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1519 if (probTemp[i] + probErrorStat[i] > 1)
1520 probErrorStat[i] = 1. - probTemp[i];
1524 // Compute summed weight of all fractions except for the considered one
1525 Double_t summedWeight = 0.;
1526 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1527 if (i != smearSpeciesByError)
1528 summedWeight += probErrorStat[i];
1531 // Compute the weight for the other species
1533 if (summedWeight <= 1e-13) {
1534 // If this happens for some reason (it should not!), just assume flat weight
1535 printf("Error: summedWeight (stat error) ~ 0 for trackPt %f, jetPt %f, centralityPercentile %f. Setting flat weight!\n",
1536 trackPt, jetPt, centralityPercentile);
1539 Double_t weight[AliPID::kSPECIES];
1540 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1541 if (i != smearSpeciesByError) {
1542 if (summedWeight > 1e-13)
1543 weight[i] = probErrorStat[i] / summedWeight;
1545 weight[i] = probErrorStat[i] / (AliPID::kSPECIES - 1);
1549 // For the final generated fractions, set the generated value for the considered species
1550 // and the generated value minus delta times statistical weight
1551 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1552 if (i != smearSpeciesByError)
1553 prob[i] = probTemp[i] - weight[i] * deltaFraction;
1555 prob[i] = generatedFraction;
1559 // Just take the generated values
1560 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
1561 prob[i] = probTemp[i];
1565 // Should already be normalised, but make sure that it really is:
1566 Double_t probSum = 0.;
1567 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1574 if (TMath::Abs(probSum - 1.0) > 1e-4) {
1575 printf("Warning: Re-normalising sum of fractions: Sum is %e\n", probSum);
1576 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1585 //_____________________________________________________________________________
1586 const TH3D* AliAnalysisTaskPID::GetParticleFractionHisto(Int_t species, Bool_t sysError) const
1588 if (species < AliPID::kElectron || species > AliPID::kProton)
1591 return sysError ? fFractionSysErrorHists[species] : fFractionHists[species];
1595 //_____________________________________________________________________________
1596 Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(Int_t motherPDG, Double_t motherGenPt)
1598 // Strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
1599 // -> Based on function in PWGJE/AliAnalysisTaskFragmentationFunction, which uses
1600 // the following data from CMS pp @ 7 TeV inclusive (JHEP 05 (2011) 064)
1604 const Int_t absMotherPDG = TMath::Abs(motherPDG);
1606 if (absMotherPDG == 310 || absMotherPDG == 321) { // K0s / K+ / K-
1607 if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.768049;
1608 else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.732933;
1609 else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.650298;
1610 else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.571332;
1611 else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.518734;
1612 else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.492543;
1613 else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.482704;
1614 else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.488056;
1615 else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.488861;
1616 else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.492862;
1617 else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.504332;
1618 else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.501858;
1619 else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.512970;
1620 else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.524131;
1621 else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.539130;
1622 else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.554101;
1623 else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.560348;
1624 else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.568869;
1625 else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.583310;
1626 else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.604818;
1627 else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.632630;
1628 else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.710070;
1629 else if(6.00 <= motherGenPt && motherGenPt < 8.00) fac = 0.736365;
1630 else if(8.00 <= motherGenPt && motherGenPt < 10.00) fac = 0.835865;
1633 if (absMotherPDG == 3122) { // Lambda
1634 if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.645162;
1635 else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.627431;
1636 else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.457136;
1637 else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.384369;
1638 else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.330597;
1639 else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.309571;
1640 else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.293620;
1641 else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.283709;
1642 else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.282047;
1643 else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.277261;
1644 else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.275772;
1645 else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.280726;
1646 else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.288540;
1647 else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.288315;
1648 else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.296619;
1649 else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.302993;
1650 else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.338121;
1651 else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.349800;
1652 else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.356802;
1653 else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.391202;
1654 else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.422573;
1655 else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.573815;
1656 else if(6.00 <= motherGenPt && motherGenPt < 8.00) fac = 0.786984;
1657 else if(8.00 <= motherGenPt && motherGenPt < 10.00) fac = 1.020021;
1660 if (absMotherPDG == 3312 || absMotherPDG == 3322) { // xi
1661 if (0.00 <= motherGenPt && motherGenPt < 0.20) fac = 0.666620;
1662 else if(0.20 <= motherGenPt && motherGenPt < 0.40) fac = 0.575908;
1663 else if(0.40 <= motherGenPt && motherGenPt < 0.60) fac = 0.433198;
1664 else if(0.60 <= motherGenPt && motherGenPt < 0.80) fac = 0.340901;
1665 else if(0.80 <= motherGenPt && motherGenPt < 1.00) fac = 0.290896;
1666 else if(1.00 <= motherGenPt && motherGenPt < 1.20) fac = 0.236074;
1667 else if(1.20 <= motherGenPt && motherGenPt < 1.40) fac = 0.218681;
1668 else if(1.40 <= motherGenPt && motherGenPt < 1.60) fac = 0.207763;
1669 else if(1.60 <= motherGenPt && motherGenPt < 1.80) fac = 0.222848;
1670 else if(1.80 <= motherGenPt && motherGenPt < 2.00) fac = 0.208806;
1671 else if(2.00 <= motherGenPt && motherGenPt < 2.20) fac = 0.197275;
1672 else if(2.20 <= motherGenPt && motherGenPt < 2.40) fac = 0.183645;
1673 else if(2.40 <= motherGenPt && motherGenPt < 2.60) fac = 0.188788;
1674 else if(2.60 <= motherGenPt && motherGenPt < 2.80) fac = 0.188282;
1675 else if(2.80 <= motherGenPt && motherGenPt < 3.00) fac = 0.207442;
1676 else if(3.00 <= motherGenPt && motherGenPt < 3.20) fac = 0.240388;
1677 else if(3.20 <= motherGenPt && motherGenPt < 3.40) fac = 0.241916;
1678 else if(3.40 <= motherGenPt && motherGenPt < 3.60) fac = 0.208276;
1679 else if(3.60 <= motherGenPt && motherGenPt < 3.80) fac = 0.234550;
1680 else if(3.80 <= motherGenPt && motherGenPt < 4.00) fac = 0.251689;
1681 else if(4.00 <= motherGenPt && motherGenPt < 5.00) fac = 0.310204;
1682 else if(5.00 <= motherGenPt && motherGenPt < 6.00) fac = 0.343492;
1685 const Double_t weight = 1. / fac;
1691 //_____________________________________________________________________________
1692 Double_t AliAnalysisTaskPID::GetMCStrangenessFactorCMS(AliMCEvent* mcEvent, AliMCParticle* daughter)
1694 // Strangeness ratio MC/data as function of mother pt from CMS data in |eta|<2.0
1695 // -> Based on function in PWGJE/AliAnalysisTaskFragmentationFunction
1700 AliMCParticle* currentMother = daughter;
1701 AliMCParticle* currentDaughter = daughter;
1704 // find first primary mother K0s, Lambda or Xi
1706 Int_t daughterPDG = currentDaughter->PdgCode();
1708 Int_t motherLabel = currentDaughter->GetMother();
1709 if(motherLabel >= mcEvent->GetNumberOfTracks()){ // protection
1710 currentMother = currentDaughter;
1714 currentMother = (AliMCParticle*)mcEvent->GetTrack(motherLabel);
1716 if (!currentMother) {
1717 currentMother = currentDaughter;
1721 Int_t motherPDG = currentMother->PdgCode();
1723 // phys. primary found ?
1724 if (mcEvent->IsPhysicalPrimary(motherLabel))
1727 if (TMath::Abs(daughterPDG) == 321) {
1728 // K+/K- e.g. from phi (ref data not feeddown corrected)
1729 currentMother = currentDaughter;
1732 if (TMath::Abs(motherPDG) == 310) {
1733 // K0s e.g. from phi (ref data not feeddown corrected)
1736 if (TMath::Abs(motherPDG) == 3212 && TMath::Abs(daughterPDG) == 3122) {
1737 // Mother Sigma0, daughter Lambda (this case not included in feeddown corr.)
1738 currentMother = currentDaughter;
1742 currentDaughter = currentMother;
1746 Int_t motherPDG = currentMother->PdgCode();
1747 Double_t motherGenPt = currentMother->Pt();
1749 return GetMCStrangenessFactorCMS(motherPDG, motherGenPt);
1753 // _________________________________________________________________________________
1754 AliAnalysisTaskPID::TOFpidInfo AliAnalysisTaskPID::GetTOFType(const AliVTrack* track, Int_t tofMode) const
1756 // Get the (locally defined) particle type judged by TOF
1758 if (!fPIDResponse) {
1759 Printf("ERROR: fEvent not available -> Cannot determine TOF type!");
1763 // Check kTOFout, kTIME, mismatch
1764 const AliPIDResponse::EDetPidStatus tofStatus = fPIDResponse->CheckPIDStatus(AliPIDResponse::kTOF, track);
1765 if (tofStatus != AliPIDResponse::kDetPidOk)
1768 Double_t nsigma[kNumTOFspecies + 1] = { -999., -999., -999., -999. };
1769 nsigma[kTOFpion] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion);
1770 nsigma[kTOFkaon] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon);
1771 nsigma[kTOFproton] = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton);
1773 Double_t inclusion = -999;
1774 Double_t exclusion = -999;
1780 else if (tofMode == 1) { // default
1784 else if (tofMode == 2) {
1789 Printf("ERROR: Bad TOF mode: %d!", tofMode);
1793 // Do not cut on nSigma electron because this would also remove pions in a large pT range.
1794 // The electron contamination of the pion sample can be treated by TPC dEdx fits afterwards.
1795 if (TMath::Abs(nsigma[kTOFpion]) < inclusion && TMath::Abs(nsigma[kTOFkaon]) > exclusion && TMath::Abs(nsigma[kTOFproton]) > exclusion)
1797 if (TMath::Abs(nsigma[kTOFpion]) > exclusion && TMath::Abs(nsigma[kTOFkaon]) < inclusion && TMath::Abs(nsigma[kTOFproton]) > exclusion)
1799 if (TMath::Abs(nsigma[kTOFpion]) > exclusion && TMath::Abs(nsigma[kTOFkaon]) > exclusion && TMath::Abs(nsigma[kTOFproton]) < inclusion)
1802 // There are no TOF electrons selected because the purity is rather bad, even for small momenta
1803 // (also a small mismatch probability significantly affects electrons because their fraction
1804 // is small). There is no need for TOF electrons anyway, since the dEdx distribution of kaons and
1805 // protons in a given pT bin (also at the dEdx crossings) is very different from that of electrons.
1806 // This is due to the steeply falling dEdx of p and K with momentum, whereas the electron dEdx stays
1808 // As a result, the TPC fit yields a more accurate electron fraction than the TOF selection can do.
1814 // _________________________________________________________________________________
1815 Bool_t AliAnalysisTaskPID::IsSecondaryWithStrangeMotherMC(AliMCEvent* mcEvent, Int_t partLabel)
1817 // Check whether particle is a secondary with strange mother, i.e. returns kTRUE if a strange mother is found
1818 // and the particle is NOT a physical primary. In all other cases kFALSE is returned
1820 if (!mcEvent || partLabel < 0)
1823 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(partLabel);
1828 if (mcEvent->IsPhysicalPrimary(partLabel))
1831 Int_t iMother = part->GetMother();
1836 AliMCParticle* partM = (AliMCParticle*)mcEvent->GetTrack(iMother);
1840 Int_t codeM = TMath::Abs(partM->PdgCode());
1841 Int_t mfl = Int_t(codeM / TMath::Power(10, Int_t(TMath::Log10(codeM))));
1842 if (mfl == 3 && codeM != 3) // codeM = 3 is for s quark
1849 //_____________________________________________________________________________
1850 Bool_t AliAnalysisTaskPID::SetParticleFractionHisto(const TH3D* hist, Int_t species, Bool_t sysError)
1852 // Store a clone of hist (containing the particle fractions of the corresponding species with statistical error (sysError = kFALSE)
1853 // or systematic error (sysError = kTRUE), respectively), internally
1855 if (species < AliPID::kElectron || species > AliPID::kProton) {
1856 AliError(Form("Only fractions for species index %d to %d can be set, but not for the requested one: %d", 0,
1857 AliPID::kProton, species));
1862 delete fFractionSysErrorHists[species];
1864 fFractionSysErrorHists[species] = new TH3D(*hist);
1867 delete fFractionHists[species];
1869 fFractionHists[species] = new TH3D(*hist);
1876 //_____________________________________________________________________________
1877 Bool_t AliAnalysisTaskPID::SetParticleFractionHistosFromFile(const TString filePathName, Bool_t sysError)
1879 // Loads particle fractions for all species from the desired file and returns kTRUE on success.
1880 // The maps are assumed to be of Type TH3D, to sit in the main directory and to have names
1881 // Form("hFraction_%e", AliPID::ParticleName(i)) for sysError = kFALSE and
1882 // Form("hFractionSysError_%e", AliPID::ParticleName(i)) for sysError = kTRUE.
1884 TFile* f = TFile::Open(filePathName.Data());
1886 std::cout << "Failed to open file with particle fractions \"" << filePathName.Data() << "\"!" << std::endl;
1891 for (Int_t i = 0; i < AliPID::kSPECIES; i++) {
1892 TString histName = Form("hFraction%s_%s", sysError ? "SysError" : "", AliPID::ParticleName(i));
1893 hist = dynamic_cast<TH3D*>(f->Get(histName.Data()));
1895 std::cout << "Failed to load particle fractions for " << histName.Data() << "!";
1896 std::cout << std::endl << "Cleaning up particle fraction histos!" << std::endl;
1897 CleanupParticleFractionHistos();
1901 if (!SetParticleFractionHisto(hist, i, sysError)) {
1902 std::cout << "Failed to load particle fractions for " << histName.Data() << "!";
1903 std::cout << std::endl << "Cleaning up particle fraction histos!" << std::endl;
1904 CleanupParticleFractionHistos();
1916 //_____________________________________________________________________________
1917 Int_t AliAnalysisTaskPID::GetRandomParticleTypeAccordingToParticleFractions(Double_t trackPt, Double_t jetPt,
1918 Double_t centralityPercentile,
1919 Bool_t smearByError,
1920 Bool_t takeIntoAccountSysError) const
1922 // Uses the stored histograms with the particle fractions to generate a random particle type according to these fractions.
1923 // In case of problems (e.g. histo missing), AliPID::kUnknown is returned.
1924 // If smearByError is kTRUE, the used fractions will be random numbers distributed with a gauss with mean
1925 // being the corresponding particle fraction and sigma it's error.
1926 // Note that in this case only the fraction of a random species is varied in this way. The other fractions
1927 // will be re-normalised according their statistical errors.
1928 // The same holds for the systematic error of species "takeIntoAccountSpeciesSysError", but the random number will be
1929 // uniformly distributed within [mean - sys, mean + sys] and the re-normalisation will be weighted with the systematic errors.
1930 // Note that the fractions will be calculated first with only the systematic error taken into account (if desired), including
1931 // re-normalisation. Then, the resulting fractions will be used to calculate the final fractions - either with statistical error
1932 // or without. The species, for which the error will be used for smearing, is the same for sys and stat error.
1934 Double_t prob[AliPID::kSPECIES];
1935 Int_t randomSpecies = (smearByError || takeIntoAccountSysError) ? (Int_t)(fRandom->Rndm() * AliPID::kSPECIES) : -1;
1936 Bool_t success = GetParticleFractions(trackPt, jetPt, centralityPercentile, prob, randomSpecies, randomSpecies);
1939 return AliPID::kUnknown;
1941 Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
1943 if (rnd <= prob[AliPID::kPion])
1944 return AliPID::kPion;
1945 else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon])
1946 return AliPID::kKaon;
1947 else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon] + prob[AliPID::kProton])
1948 return AliPID::kProton;
1949 else if (rnd <= prob[AliPID::kPion] + prob[AliPID::kKaon] + prob[AliPID::kProton] + prob[AliPID::kElectron])
1950 return AliPID::kElectron;
1952 return AliPID::kMuon; //else it must be a muon (only species left)
1956 //_____________________________________________________________________________
1957 AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::GenerateDetectorResponse(AliAnalysisTaskPID::ErrorCode errCode,
1958 Double_t mean, Double_t sigma,
1959 Double_t* responses, Int_t nResponses,
1962 // Generate detector response. If a previous generation was not successful or there is something wrong with this signal generation,
1963 // the function will return kFALSE
1967 // Reset response array
1968 for (Int_t i = 0; i < nResponses; i++)
1969 responses[i] = -999;
1971 if (errCode == kError)
1974 ErrorCode ownErrCode = kNoErrors;
1976 if (fUseConvolutedGaus && !usePureGaus) {
1977 // In case of convoluted gauss, calculate the probability density only once to save a lot of time!
1979 TH1* hProbDensity = 0x0;
1980 ownErrCode = SetParamsForConvolutedGaus(mean, sigma);
1981 if (ownErrCode == kError)
1984 hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
1986 for (Int_t i = 0; i < nResponses; i++) {
1987 responses[i] = hProbDensity->GetRandom();
1988 //responses[i] fConvolutedGausDeltaPrime->GetRandom(); // MUCH slower than using the binned version via the histogram
1992 for (Int_t i = 0; i < nResponses; i++) {
1993 responses[i] = fRandom->Gaus(mean, sigma);
1997 // If forwarded error code was a warning (error case has been handled before), return a warning
1998 if (errCode == kWarning)
2001 return ownErrCode; // Forward success/warning
2005 //_____________________________________________________________________________
2006 void AliAnalysisTaskPID::PrintSettings(Bool_t printSystematicsSettings) const
2008 // Print current settings.
2010 printf("\n\nSettings for task %s:\n", GetName());
2011 printf("Is pPb/Pbp: %d -> %s\n", GetIsPbpOrpPb(), GetIsPbpOrpPb() ? "Adapting vertex cuts" : "Using standard vertex cuts");
2012 printf("Track cuts: %s\n", fTrackFilter ? fTrackFilter->GetTitle() : "-");
2013 printf("Eta cut: %.2f <= |eta| <= %.2f\n", GetEtaAbsCutLow(), GetEtaAbsCutUp());
2014 printf("Phi' cut: %d\n", GetUsePhiCut());
2015 printf("TPCCutMIGeo: %d\n", GetUseTPCCutMIGeo());
2016 if (GetUseTPCCutMIGeo()) {
2017 printf("GetCutGeo: %f\n", GetCutGeo());
2018 printf("GetCutNcr: %f\n", GetCutNcr());
2019 printf("GetCutNcl: %f\n", GetCutNcl());
2021 printf("TPCnclCut: %d\n", GetUseTPCnclCut());
2022 if (GetUseTPCnclCut()) {
2023 printf("GetCutPureNcl: %d\n", GetCutPureNcl());
2028 printf("Centrality estimator: %s\n", GetCentralityEstimator().Data());
2032 printf("Use MC-ID for signal generation: %d\n", GetUseMCidForGeneration());
2033 printf("Use ITS: %d\n", GetUseITS());
2034 printf("Use TOF: %d\n", GetUseTOF());
2035 printf("Use priors: %d\n", GetUsePriors());
2036 printf("Use TPC default priors: %d\n", GetUseTPCDefaultPriors());
2037 printf("Use convoluted Gauss: %d\n", GetUseConvolutedGaus());
2038 printf("Accuracy of non-Gaussian tail: %e\n", GetAccuracyNonGaussianTail());
2039 printf("Take into account muons: %d\n", GetTakeIntoAccountMuons());
2040 printf("TOF mode: %d\n", GetTOFmode());
2041 printf("\nParams for transition from gauss to asymmetric shape:\n");
2042 printf("[0]: %e\n", GetConvolutedGaussTransitionPar(0));
2043 printf("[1]: %e\n", GetConvolutedGaussTransitionPar(1));
2044 printf("[2]: %e\n", GetConvolutedGaussTransitionPar(2));
2048 printf("Do PID: %d\n", fDoPID);
2049 printf("Do Efficiency: %d\n", fDoEfficiency);
2050 printf("Do PtResolution: %d\n", fDoPtResolution);
2054 printf("Input from other task: %d\n", GetInputFromOtherTask());
2055 printf("Store additional jet information: %d\n", GetStoreAdditionalJetInformation());
2056 printf("Store centrality percentile: %d", GetStoreCentralityPercentile());
2058 if (printSystematicsSettings)
2059 PrintSystematicsSettings();
2065 //_____________________________________________________________________________
2066 void AliAnalysisTaskPID::PrintSystematicsSettings() const
2068 // Print current settings for systematic studies.
2070 printf("\n\nSettings for systematics for task %s:\n", GetName());
2071 printf("SplinesThr:\t%f\n", GetSystematicScalingSplinesThreshold());
2072 printf("SplinesBelowThr:\t%f\n", GetSystematicScalingSplinesBelowThreshold());
2073 printf("SplinesAboveThr:\t%f\n", GetSystematicScalingSplinesAboveThreshold());
2074 printf("EtaCorrMomThr:\t%f\n", GetSystematicScalingEtaCorrectionMomentumThr());
2075 printf("EtaCorrLowP:\t%f\n", GetSystematicScalingEtaCorrectionLowMomenta());
2076 printf("EtaCorrHighP:\t%f\n", GetSystematicScalingEtaCorrectionHighMomenta());
2077 printf("SigmaParaThr:\t%f\n", GetSystematicScalingEtaSigmaParaThreshold());
2078 printf("SigmaParaBelowThr:\t%f\n", GetSystematicScalingEtaSigmaParaBelowThreshold());
2079 printf("SigmaParaAboveThr:\t%f\n", GetSystematicScalingEtaSigmaParaAboveThreshold());
2080 printf("MultCorr:\t%f\n", GetSystematicScalingMultCorrection());
2081 printf("TOF mode: %d\n", GetTOFmode());
2087 //_____________________________________________________________________________
2088 Bool_t AliAnalysisTaskPID::ProcessTrack(const AliVTrack* track, Int_t particlePDGcode, Double_t centralityPercentile,
2091 // Process the track (generate expected response, fill histos, etc.).
2092 // particlePDGcode == 0 means data. Otherwise, the corresponding MC ID will be assumed.
2094 //Printf("Debug: Task %s is starting to process track: dEdx %f, pTPC %f, eta %f, ncl %d\n", GetName(), track->GetTPCsignal(), track->GetTPCmomentum(),
2095 // track->Eta(), track->GetTPCsignalN());
2098 printf("File: %s, Line: %d: ProcessTrack\n", (char*)__FILE__, __LINE__);
2104 printf("File: %s, Line: %d: ProcessTrack -> Processing started\n", (char*)__FILE__, __LINE__);
2106 const Bool_t isMC = (particlePDGcode == 0) ? kFALSE : kTRUE;
2111 if (TMath::Abs(particlePDGcode) == 211) {//Pion
2114 else if (TMath::Abs(particlePDGcode) == 321) {//Kaon
2117 else if (TMath::Abs(particlePDGcode) == 2212) {//Proton
2120 else if (TMath::Abs(particlePDGcode) == 11) {//Electron
2123 else if (TMath::Abs(particlePDGcode) == 13) {//Muon
2126 else // In MC-ID case, set to underflow bin such that the response from this track is only used for unidentified signal generation
2127 // or signal generation with PID response and the track is still there (as in data) - e.g. relevant w.r.t. deuterons.
2128 // This is important to be as much as possible consistent with data. And the tracks can still be removed by disabling the
2129 // underflow bin for the projections
2134 //Double_t p = track->GetP();
2135 const Double_t pTPC = track->GetTPCmomentum();
2136 Double_t pT = track->Pt();
2138 Double_t z = -1, xi = -1;
2139 GetJetTrackObservables(pT, jetPt, z, xi);
2142 Double_t trackCharge = track->Charge();
2145 Double_t dEdxTPC = fPIDResponse->IsTunedOnData() ? fPIDResponse->GetTPCsignalTunedOnData(track) : track->GetTPCsignal();
2148 Printf("Skipping track with strange dEdx value: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), pTPC,
2149 track->Eta(), track->GetTPCsignalN());
2153 // Completely remove tracks from light nuclei defined by el and pr <dEdx> as function of pT.
2154 // A very loose cut to be sure not to throw away electrons and/or protons
2155 Double_t nSigmaPr = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton);
2156 Double_t nSigmaEl = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
2158 if ((pTPC >= 0.3 && (nSigmaPr > 10 && nSigmaEl > 10)) ||
2159 (pTPC < 0.3 && (nSigmaPr > 15 && nSigmaEl > 15))) {
2160 Printf("Skipping track which seems to be a light nuclei: dEdx %f, pTPC %f, pT %f, eta %f, ncl %d, nSigmaPr %f, nSigmaEl %f\n",
2161 track->GetTPCsignal(), pTPC, pT, track->Eta(), track->GetTPCsignalN(), nSigmaPr, nSigmaEl);
2166 Double_t dEdxEl, dEdxKa, dEdxPi, dEdxMu, dEdxPr;
2167 Double_t sigmaEl, sigmaKa, sigmaPi, sigmaMu, sigmaPr;
2169 if (fDoAnySystematicStudiesOnTheExpectedSignal) {
2170 // Get the uncorrected signal first and the corresponding correction factors.
2171 // Then modify the correction factors and properly recalculate the corrected dEdx
2173 // Get pure spline values for dEdx_expected, without any correction
2174 dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2175 dEdxKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2176 dEdxPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2177 dEdxMu = !fTakeIntoAccountMuons ? -1 :
2178 fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2179 dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault, kFALSE, kFALSE);
2181 // Scale splines, if desired
2182 if ((TMath::Abs(fSystematicScalingSplinesBelowThreshold - 1.0) > fgkEpsilon) ||
2183 (TMath::Abs(fSystematicScalingSplinesAboveThreshold - 1.0) > fgkEpsilon)) {
2185 // Tune turn-on of correction for pions (not so relevant for the other species, since at very large momenta)
2187 Double_t scaleFactor = 1.;
2188 Double_t usedSystematicScalingSplines = 1.;
2190 bg = pTPC / AliPID::ParticleMass(AliPID::kElectron);
2191 scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2192 usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2193 + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2194 dEdxEl *= usedSystematicScalingSplines;
2196 bg = pTPC / AliPID::ParticleMass(AliPID::kKaon);
2197 scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2198 usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2199 + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2200 dEdxKa *= usedSystematicScalingSplines;
2202 bg = pTPC / AliPID::ParticleMass(AliPID::kPion);
2203 scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2204 usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2205 + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2206 dEdxPi *= usedSystematicScalingSplines;
2208 if (fTakeIntoAccountMuons) {
2209 bg = pTPC / AliPID::ParticleMass(AliPID::kMuon);
2210 scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2211 usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2212 + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2213 dEdxMu *= usedSystematicScalingSplines;
2217 bg = pTPC / AliPID::ParticleMass(AliPID::kProton);
2218 scaleFactor = 0.5 * (1. + TMath::Erf((bg - fSystematicScalingSplinesThreshold) / 4.));
2219 usedSystematicScalingSplines = fSystematicScalingSplinesBelowThreshold * (1 - scaleFactor)
2220 + fSystematicScalingSplinesAboveThreshold * scaleFactor;
2221 dEdxPr *= usedSystematicScalingSplines;
2224 // Get the eta correction factors for the (modified) expected dEdx
2225 Double_t etaCorrEl = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxEl) : 1.;
2226 Double_t etaCorrKa = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxKa) : 1.;
2227 Double_t etaCorrPi = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPi) : 1.;
2228 Double_t etaCorrMu = fTakeIntoAccountMuons && !fPIDResponse->UseTPCEtaCorrection() ?
2229 fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxMu) : 1.;
2230 Double_t etaCorrPr = fPIDResponse->UseTPCEtaCorrection() ? fPIDResponse->GetTPCResponse().GetEtaCorrectionFast(track, dEdxPr) : 1.;
2232 // Scale eta correction factors, if desired (and eta correction maps are to be used, otherwise it is not possible!)
2233 if (fPIDResponse->UseTPCEtaCorrection() &&
2234 (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - 1.0) > fgkEpsilon ||
2235 TMath::Abs(fSystematicScalingEtaCorrectionLowMomenta - 1.0) > fgkEpsilon)) {
2236 // Since we do not want to scale the splines with this, but only the eta variation, only scale the deviation of the correction factor!
2237 // 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!
2240 // Due to additional azimuthal effects, there is an additional eta dependence for low momenta which is not corrected successfully so far.
2241 // One can assign a different (higher) systematic scale factor for this low-p region and a threshold which separates low- and high-p.
2242 // 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
2243 Double_t usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionHighMomenta;
2245 if (TMath::Abs(fSystematicScalingEtaCorrectionHighMomenta - fSystematicScalingEtaCorrectionLowMomenta) > fgkEpsilon) {
2246 const Double_t fractionHighMomentumScaleFactor = 0.5 * (1. + TMath::Erf((pTPC - fSystematicScalingEtaCorrectionMomentumThr) / 0.1));
2247 usedSystematicScalingEtaCorrection = fSystematicScalingEtaCorrectionLowMomenta * (1 - fractionHighMomentumScaleFactor)
2248 + fSystematicScalingEtaCorrectionHighMomenta * fractionHighMomentumScaleFactor;
2251 etaCorrEl = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrEl - 1.0);
2252 etaCorrKa = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrKa - 1.0);
2253 etaCorrPi = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPi - 1.0);
2254 etaCorrMu = fTakeIntoAccountMuons ? (1.0 + usedSystematicScalingEtaCorrection * (etaCorrMu - 1.0)) : 1.0;
2255 etaCorrPr = 1.0 + usedSystematicScalingEtaCorrection * (etaCorrPr - 1.0);
2258 // Get the multiplicity correction factors for the (modified) expected dEdx
2259 const Int_t currEvtMultiplicity = fPIDResponse->GetTPCResponse().GetCurrentEventMultiplicity();
2261 Double_t multiplicityCorrEl = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2262 dEdxEl * etaCorrEl, currEvtMultiplicity) : 1.;
2263 Double_t multiplicityCorrKa = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2264 dEdxKa * etaCorrKa, currEvtMultiplicity) : 1.;
2265 Double_t multiplicityCorrPi = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2266 dEdxPi * etaCorrPi, currEvtMultiplicity) : 1.;
2267 Double_t multiplicityCorrMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCMultiplicityCorrection() ?
2268 fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track, dEdxMu * etaCorrMu, currEvtMultiplicity) : 1.;
2269 Double_t multiplicityCorrPr = fPIDResponse->UseTPCMultiplicityCorrection() ? fPIDResponse->GetTPCResponse().GetMultiplicityCorrectionFast(track,
2270 dEdxPr * etaCorrPr, currEvtMultiplicity) : 1.;
2272 Double_t multiplicityCorrSigmaEl = fPIDResponse->UseTPCMultiplicityCorrection() ?
2273 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxEl * etaCorrEl, currEvtMultiplicity) : 1.;
2274 Double_t multiplicityCorrSigmaKa = fPIDResponse->UseTPCMultiplicityCorrection() ?
2275 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxKa * etaCorrKa, currEvtMultiplicity) : 1.;
2276 Double_t multiplicityCorrSigmaPi = fPIDResponse->UseTPCMultiplicityCorrection() ?
2277 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxPi * etaCorrPi, currEvtMultiplicity) : 1.;
2278 Double_t multiplicityCorrSigmaMu = fTakeIntoAccountMuons && fPIDResponse->UseTPCMultiplicityCorrection() ?
2279 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxMu * etaCorrMu, currEvtMultiplicity) : 1.;
2280 Double_t multiplicityCorrSigmaPr = fPIDResponse->UseTPCMultiplicityCorrection() ?
2281 fPIDResponse->GetTPCResponse().GetMultiplicitySigmaCorrectionFast(dEdxPr * etaCorrPr, currEvtMultiplicity) : 1.;
2283 // Scale multiplicity correction factors, if desired (and multiplicity correction functions are to be used, otherwise it is not possible!)
2284 if (fPIDResponse->UseTPCMultiplicityCorrection() && TMath::Abs(fSystematicScalingMultCorrection - 1.0) > fgkEpsilon) {
2285 // Since we do not want to scale the splines with this, but only the multiplicity variation, only scale the deviation of the correction factor!
2286 // 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!
2288 multiplicityCorrEl = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrEl - 1.0);
2289 multiplicityCorrKa = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrKa - 1.0);
2290 multiplicityCorrPi = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrPi - 1.0);
2291 multiplicityCorrMu = fTakeIntoAccountMuons ? (1.0 + fSystematicScalingMultCorrection * (multiplicityCorrMu - 1.0)) : 1.0;
2292 multiplicityCorrPr = 1.0 + fSystematicScalingMultCorrection * (multiplicityCorrPr - 1.0);
2295 // eta correction must be enabled in order to use the new sigma parametrisation maps. Since this is the absolute sigma
2296 // for a track calculated with the unscaled paramaters, we have to devide out dEdxExpectedEtaCorrected and then need
2297 // to scale with the multiplicitySigmaCorrFactor * fSystematicScalingEtaSigmaPara. In the end, one has to scale with the
2298 // (modified) dEdx to get the absolute sigma
2299 // This means there is no extra parameter for the multiplicitySigmaCorrFactor, but only for the sigma map itself.
2300 // This is valid, since it appears only as a product. One has to assume a larger systematic shift in case of additional
2301 // multiplicity dependence....
2303 Bool_t doSigmaSystematics = (TMath::Abs(fSystematicScalingEtaSigmaParaBelowThreshold - 1.0) > fgkEpsilon) ||
2304 (TMath::Abs(fSystematicScalingEtaSigmaParaAboveThreshold - 1.0) > fgkEpsilon);
2307 Double_t dEdxElExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2309 Double_t systematicScalingEtaSigmaParaEl = 1.;
2310 if (doSigmaSystematics) {
2311 Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxElExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2312 systematicScalingEtaSigmaParaEl = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2313 + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2315 Double_t sigmaRelEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2317 / dEdxElExpected * systematicScalingEtaSigmaParaEl * multiplicityCorrSigmaEl;
2320 Double_t dEdxKaExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2322 Double_t systematicScalingEtaSigmaParaKa = 1.;
2323 if (doSigmaSystematics) {
2324 Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxKaExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2325 systematicScalingEtaSigmaParaKa = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2326 + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2328 Double_t sigmaRelKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2330 / dEdxKaExpected * systematicScalingEtaSigmaParaKa * multiplicityCorrSigmaKa;
2333 Double_t dEdxPiExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2335 Double_t systematicScalingEtaSigmaParaPi = 1.;
2336 if (doSigmaSystematics) {
2337 Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxPiExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2338 systematicScalingEtaSigmaParaPi = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2339 + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2341 Double_t sigmaRelPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2343 / dEdxPiExpected * systematicScalingEtaSigmaParaPi * multiplicityCorrSigmaPi;
2346 Double_t sigmaRelMu = 999.;
2347 if (fTakeIntoAccountMuons) {
2348 Double_t dEdxMuExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2350 Double_t systematicScalingEtaSigmaParaMu = 1.;
2351 if (doSigmaSystematics) {
2352 Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxMuExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2353 systematicScalingEtaSigmaParaMu = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2354 + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2356 sigmaRelMu = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault, kTRUE, kFALSE)
2357 / dEdxMuExpected * systematicScalingEtaSigmaParaMu * multiplicityCorrSigmaMu;
2361 Double_t dEdxPrExpected = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2363 Double_t systematicScalingEtaSigmaParaPr = 1.;
2364 if (doSigmaSystematics) {
2365 Double_t scaleFactor = 0.5 * (1. + TMath::Erf((dEdxPrExpected - fSystematicScalingEtaSigmaParaThreshold) / 25.));
2366 systematicScalingEtaSigmaParaPr = fSystematicScalingEtaSigmaParaBelowThreshold * (1 - scaleFactor)
2367 + fSystematicScalingEtaSigmaParaAboveThreshold * scaleFactor;
2369 Double_t sigmaRelPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2371 / dEdxPrExpected * systematicScalingEtaSigmaParaPr * multiplicityCorrSigmaPr;
2373 // Now scale the (possibly modified) spline values with the (possibly modified) correction factors
2374 dEdxEl *= etaCorrEl * multiplicityCorrEl;
2375 dEdxKa *= etaCorrKa * multiplicityCorrKa;
2376 dEdxPi *= etaCorrPi * multiplicityCorrPi;
2377 dEdxMu *= etaCorrMu * multiplicityCorrMu;
2378 dEdxPr *= etaCorrPr * multiplicityCorrPr;
2380 // Finally, get the absolute sigma
2381 sigmaEl = sigmaRelEl * dEdxEl;
2382 sigmaKa = sigmaRelKa * dEdxKa;
2383 sigmaPi = sigmaRelPi * dEdxPi;
2384 sigmaMu = sigmaRelMu * dEdxMu;
2385 sigmaPr = sigmaRelPr * dEdxPr;
2389 // No systematic studies on expected signal - just take it as it comve from the TPCPIDResponse
2390 dEdxEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2391 fPIDResponse->UseTPCEtaCorrection(),
2392 fPIDResponse->UseTPCMultiplicityCorrection());
2393 dEdxKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2394 fPIDResponse->UseTPCEtaCorrection(),
2395 fPIDResponse->UseTPCMultiplicityCorrection());
2396 dEdxPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2397 fPIDResponse->UseTPCEtaCorrection(),
2398 fPIDResponse->UseTPCMultiplicityCorrection());
2399 dEdxMu = !fTakeIntoAccountMuons ? -1 :
2400 fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2401 fPIDResponse->UseTPCEtaCorrection(),
2402 fPIDResponse->UseTPCMultiplicityCorrection());
2403 dEdxPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2404 fPIDResponse->UseTPCEtaCorrection(),
2405 fPIDResponse->UseTPCMultiplicityCorrection());
2407 sigmaEl = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
2408 fPIDResponse->UseTPCEtaCorrection(),
2409 fPIDResponse->UseTPCMultiplicityCorrection());
2410 sigmaKa = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
2411 fPIDResponse->UseTPCEtaCorrection(),
2412 fPIDResponse->UseTPCMultiplicityCorrection());
2413 sigmaPi = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
2414 fPIDResponse->UseTPCEtaCorrection(),
2415 fPIDResponse->UseTPCMultiplicityCorrection());
2416 sigmaMu = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kMuon, AliTPCPIDResponse::kdEdxDefault,
2417 fPIDResponse->UseTPCEtaCorrection(),
2418 fPIDResponse->UseTPCMultiplicityCorrection());
2419 sigmaPr = fPIDResponse->GetTPCResponse().GetExpectedSigma(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
2420 fPIDResponse->UseTPCEtaCorrection(),
2421 fPIDResponse->UseTPCMultiplicityCorrection());
2424 Double_t deltaPrimeElectron = (dEdxEl > 0) ? dEdxTPC / dEdxEl : -1;
2426 Printf("Error: Expected TPC signal <= 0 for electron hypothesis");
2430 Double_t deltaPrimeKaon = (dEdxKa > 0) ? dEdxTPC / dEdxKa : -1;
2432 Printf("Error: Expected TPC signal <= 0 for kaon hypothesis");
2436 Double_t deltaPrimePion = (dEdxPi > 0) ? dEdxTPC / dEdxPi : -1;
2438 Printf("Error: Expected TPC signal <= 0 for pion hypothesis");
2442 Double_t deltaPrimeProton = (dEdxPr > 0) ? dEdxTPC / dEdxPr : -1;
2444 Printf("Error: Expected TPC signal <= 0 for proton hypothesis");
2449 printf("File: %s, Line: %d: ProcessTrack -> Compute probabilities\n", (char*)__FILE__, __LINE__);
2452 // Use probabilities to weigh the response generation for the different species.
2453 // Also determine the most probable particle type.
2454 Double_t prob[AliPID::kSPECIESC];
2455 for (Int_t i = 0; i < AliPID::kSPECIESC; i++)
2458 fPIDcombined->ComputeProbabilities(track, fPIDResponse, prob);
2460 // If muons are not to be taken into account, just set their probability to zero and normalise the remaining probabilities later
2461 if (!fTakeIntoAccountMuons)
2462 prob[AliPID::kMuon] = 0;
2464 // Priors might cause problems in case of mismatches, i.e. mismatch particles will be identified as pions.
2465 // Can be easily caught for low p_TPC by just setting pion (and muon) probs to zero, if dEdx is larger than
2466 // expected dEdx of electrons and kaons
2467 if (pTPC < 1. && dEdxTPC > dEdxEl && dEdxTPC > dEdxKa) {
2468 prob[AliPID::kMuon] = 0;
2469 prob[AliPID::kPion] = 0;
2473 // Bug: One can set the number of species for PIDcombined, but PIDcombined will call PIDresponse, which writes without testing
2474 // the probs for kSPECIESC (including light nuclei) into the array.
2475 // In this case, when only kSPECIES are considered, the probabilities have to be rescaled!
2478 // 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
2479 // high-pT region (also contribution there completely negligable!)
2480 // 2. If dEdx is larger than the expected dEdx + 5 sigma of el and pr, it is most likely a light nucleus
2481 // (if dEdx is more than 5 sigma away from expectation, flat TPC probs are set... and the light nuclei splines are not
2484 // The highest expected dEdx is either for pr or for el. Assign 100% probability to the species with highest expected dEdx then.
2485 // In all other cases: Throw away light nuclei probs and rescale other probs to 100%.
2487 // Find most probable species for the ID
2488 Int_t maxIndex = TMath::LocMax(AliPID::kSPECIESC, prob);
2490 if ((prob[maxIndex] > 0 && maxIndex >= AliPID::kSPECIES && dEdxTPC > dEdxEl + 5. * sigmaEl) ||
2491 (dEdxTPC > dEdxEl + 5. * sigmaEl && dEdxTPC > dEdxPr + 5. * sigmaPr)) {
2492 for (Int_t i = 0; i < AliPID::kSPECIESC; i++)
2495 if (dEdxEl > dEdxPr) {
2496 prob[AliPID::kElectron] = 1.;
2497 maxIndex = AliPID::kElectron;
2500 prob[AliPID::kProton] = 1.;
2501 maxIndex = AliPID::kProton;
2505 for (Int_t i = AliPID::kSPECIES; i < AliPID::kSPECIESC; i++)
2508 Double_t probSum = 0.;
2509 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2513 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2517 // If all probabilities are equal (should only happen if no priors are used), set most probable PID to pion
2518 // because LocMax returns just the first maximal value (i.e. AliPID::kElectron)
2519 if (TMath::Abs(prob[maxIndex] - 1. / AliPID::kSPECIES) < 1e-5)
2520 maxIndex = AliPID::kPion;
2522 // In all other cases, the maximum remains untouched from the scaling (and is < AliPID::kSPECIES by construction)
2526 // Clean up the most probable PID for low momenta (not an issue for the probabilities, since fractions small,
2527 // but to get clean (and nice looking) templates (if most probable PID is used instead of generated responses), it should be done).
2528 // Problem: If more than 5 sigma away (and since priors also included), most probable PID for dEdx around protons could be pions.
2529 // Idea: Only clean selection for K and p possible. So, just take for the most probable PID the probs from TPC only calculated
2530 // by hand without restriction to 5 sigma and without priors. This will not affect the template generation, but afterwards one
2531 // can replace the K and p templates with the most probable PID distributions, so all other species templates remain the same.
2532 // I.e. it doesn't hurt if the most probable PID is distributed incorrectly between el, pi, mu.
2533 Bool_t maxIndexSetManually = kFALSE;
2536 Double_t probManualTPC[AliPID::kSPECIES];
2537 for (Int_t i = 0; i < AliPID::kSPECIES; i++)
2538 probManualTPC[i] = 0;
2540 probManualTPC[AliPID::kElectron] = TMath::Exp(-0.5*(dEdxTPC-dEdxEl)*(dEdxTPC-dEdxEl)/(sigmaEl*sigmaEl))/sigmaEl;
2541 // Muons are not important here, just ignore and save processing time
2542 probManualTPC[AliPID::kMuon] = 0.;//TMath::Exp(-0.5*(dEdxTPC-dEdxMu)*(dEdxTPC-dEdxMu)/(sigmaMu*sigmaMu))/sigmaMu;
2543 probManualTPC[AliPID::kPion] = TMath::Exp(-0.5*(dEdxTPC-dEdxPi)*(dEdxTPC-dEdxPi)/(sigmaPi*sigmaPi))/sigmaPi;
2544 probManualTPC[AliPID::kKaon] = TMath::Exp(-0.5*(dEdxTPC-dEdxKa)*(dEdxTPC-dEdxKa)/(sigmaKa*sigmaKa))/sigmaKa;
2545 probManualTPC[AliPID::kProton] = TMath::Exp(-0.5*(dEdxTPC-dEdxPr)*(dEdxTPC-dEdxPr)/(sigmaPr*sigmaPr))/sigmaPr;
2547 const Int_t maxIndexManualTPC = TMath::LocMax(AliPID::kSPECIES, probManualTPC);
2548 // Should always be > 0, but in case the dEdx is far off such that the machine accuracy sets every prob to zero,
2549 // better take the "old" result
2550 if (probManualTPC[maxIndexManualTPC] > 0.)
2551 maxIndex = maxIndexManualTPC;
2553 maxIndexSetManually = kTRUE;
2557 // Translate from AliPID numbering to numbering of this class
2558 if (prob[maxIndex] > 0 || maxIndexSetManually) {
2559 if (maxIndex == AliPID::kElectron)
2561 else if (maxIndex == AliPID::kKaon)
2563 else if (maxIndex == AliPID::kMuon)
2565 else if (maxIndex == AliPID::kPion)
2567 else if (maxIndex == AliPID::kProton)
2573 // Only take track into account for expectation values, if valid pid response is available.. Otherwise: Set to underflow bin.
2579 //For testing: Swap p<->pT to analyse pure p-dependence => Needs to be removed later
2585 TOFpidInfo tofPIDinfo = GetTOFType(track, fTOFmode);
2587 Double_t entry[fStoreAdditionalJetInformation ? kDataNumAxes : kDataNumAxes - fgkNumJetAxes];
2588 entry[kDataMCID] = binMC;
2589 entry[kDataSelectSpecies] = 0;
2590 entry[kDataPt] = pT;
2591 entry[kDataDeltaPrimeSpecies] = deltaPrimeElectron;
2592 entry[kDataCentrality] = centralityPercentile;
2594 if (fStoreAdditionalJetInformation) {
2595 entry[kDataJetPt] = jetPt;
2597 entry[kDataXi] = xi;
2600 entry[GetIndexOfChargeAxisData()] = trackCharge;
2601 entry[GetIndexOfTOFpidInfoAxisData()] = tofPIDinfo;
2603 fhPIDdataAll->Fill(entry);
2605 entry[kDataSelectSpecies] = 1;
2606 entry[kDataDeltaPrimeSpecies] = deltaPrimeKaon;
2607 fhPIDdataAll->Fill(entry);
2609 entry[kDataSelectSpecies] = 2;
2610 entry[kDataDeltaPrimeSpecies] = deltaPrimePion;
2611 fhPIDdataAll->Fill(entry);
2613 entry[kDataSelectSpecies] = 3;
2614 entry[kDataDeltaPrimeSpecies] = deltaPrimeProton;
2615 fhPIDdataAll->Fill(entry);
2618 // Construct the expected shape for the transition p -> pT
2620 Double_t genEntry[fStoreAdditionalJetInformation ? kGenNumAxes : kGenNumAxes - fgkNumJetAxes];
2621 genEntry[kGenMCID] = binMC;
2622 genEntry[kGenSelectSpecies] = 0;
2623 genEntry[kGenPt] = pT;
2624 genEntry[kGenDeltaPrimeSpecies] = -999;
2625 genEntry[kGenCentrality] = centralityPercentile;
2627 if (fStoreAdditionalJetInformation) {
2628 genEntry[kGenJetPt] = jetPt;
2629 genEntry[kGenZ] = z;
2630 genEntry[kGenXi] = xi;
2633 genEntry[GetIndexOfChargeAxisGen()] = trackCharge;
2634 genEntry[GetIndexOfTOFpidInfoAxisGen()] = tofPIDinfo;
2636 // Generate numGenEntries "responses" with fluctuations around the expected values.
2637 // fgkMaxNumGenEntries = 500 turned out to give reasonable templates even for highest track pT in 15-20 GeV/c jet pT bin.
2638 Int_t numGenEntries = fgkMaxNumGenEntries; // fgkMaxNumGenEntries = 500
2640 /*OLD: Different number of responses depending on pT and jet pT for fgkMaxNumGenEntries = 1000
2641 * => Problem: If threshold to higher number of responses inside a bin (or after rebinning), then the template
2642 * is biased to the higher pT.
2643 // Generate numGenEntries "responses" with fluctuations around the expected values.
2644 // The higher the (transverse) momentum, the more "responses" will be generated in order not to run out of statistics too fast.
2645 Int_t numGenEntries = 10;
2647 // Jets have even worse statistics, therefore, scale numGenEntries further
2649 numGenEntries *= 20;
2650 else if (jetPt >= 20)
2651 numGenEntries *= 10;
2652 else if (jetPt >= 10)
2657 // Do not generate more entries than available in memory!
2658 if (numGenEntries > fgkMaxNumGenEntries)// fgkMaxNumGenEntries = 1000
2659 numGenEntries = fgkMaxNumGenEntries;
2663 ErrorCode errCode = kNoErrors;
2666 printf("File: %s, Line: %d: ProcessTrack -> Generate Responses\n", (char*)__FILE__, __LINE__);
2669 errCode = GenerateDetectorResponse(errCode, 1., sigmaEl / dEdxEl, fGenRespElDeltaPrimeEl, numGenEntries);
2670 errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxKa, sigmaEl / dEdxKa, fGenRespElDeltaPrimeKa, numGenEntries);
2671 errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxPi, sigmaEl / dEdxPi, fGenRespElDeltaPrimePi, numGenEntries);
2672 errCode = GenerateDetectorResponse(errCode, dEdxEl / dEdxPr, sigmaEl / dEdxPr, fGenRespElDeltaPrimePr, numGenEntries);
2675 errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxEl, sigmaKa / dEdxEl, fGenRespKaDeltaPrimeEl, numGenEntries);
2676 errCode = GenerateDetectorResponse(errCode, 1., sigmaKa / dEdxKa, fGenRespKaDeltaPrimeKa, numGenEntries);
2677 errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxPi, sigmaKa / dEdxPi, fGenRespKaDeltaPrimePi, numGenEntries);
2678 errCode = GenerateDetectorResponse(errCode, dEdxKa / dEdxPr, sigmaKa / dEdxPr, fGenRespKaDeltaPrimePr, numGenEntries);
2681 errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxEl, sigmaPi / dEdxEl, fGenRespPiDeltaPrimeEl, numGenEntries);
2682 errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxKa, sigmaPi / dEdxKa, fGenRespPiDeltaPrimeKa, numGenEntries);
2683 errCode = GenerateDetectorResponse(errCode, 1., sigmaPi / dEdxPi, fGenRespPiDeltaPrimePi, numGenEntries);
2684 errCode = GenerateDetectorResponse(errCode, dEdxPi / dEdxPr, sigmaPi / dEdxPr, fGenRespPiDeltaPrimePr, numGenEntries);
2686 // Muons, if desired
2687 if (fTakeIntoAccountMuons) {
2688 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxEl, sigmaMu / dEdxEl, fGenRespMuDeltaPrimeEl, numGenEntries);
2689 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxKa, sigmaMu / dEdxKa, fGenRespMuDeltaPrimeKa, numGenEntries);
2690 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxPi, sigmaMu / dEdxPi, fGenRespMuDeltaPrimePi, numGenEntries);
2691 errCode = GenerateDetectorResponse(errCode, dEdxMu / dEdxPr, sigmaMu / dEdxPr, fGenRespMuDeltaPrimePr, numGenEntries);
2695 errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxEl, sigmaPr / dEdxEl, fGenRespPrDeltaPrimeEl, numGenEntries);
2696 errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxKa, sigmaPr / dEdxKa, fGenRespPrDeltaPrimeKa, numGenEntries);
2697 errCode = GenerateDetectorResponse(errCode, dEdxPr / dEdxPi, sigmaPr / dEdxPi, fGenRespPrDeltaPrimePi, numGenEntries);
2698 errCode = GenerateDetectorResponse(errCode, 1., sigmaPr / dEdxPr, fGenRespPrDeltaPrimePr, numGenEntries);
2700 if (errCode != kNoErrors) {
2701 if (errCode == kWarning && fDebug > 1) {
2702 Printf("Warning: Questionable detector response for track, most likely due to very low number of PID clusters! Debug output (dEdx_expected, sigma_expected):");
2705 Printf("Error: Failed to generate detector response for track - skipped! Debug output (dEdx_expected, sigma_expected):");
2708 Printf("Pr: %e, %e", dEdxPr, sigmaPr);
2709 Printf("Pi: %e, %e", dEdxPi, sigmaPi);
2710 Printf("El: %e, %e", dEdxEl, sigmaEl);
2711 Printf("Mu: %e, %e", dEdxMu, sigmaMu);
2712 Printf("Ka: %e, %e", dEdxKa, sigmaKa);
2713 Printf("track: dEdx %f, pTPC %f, eta %f, ncl %d\n", track->GetTPCsignal(), track->GetTPCmomentum(), track->Eta(),
2714 track->GetTPCsignalN());
2717 if (errCode != kWarning) {
2718 fhSkippedTracksForSignalGeneration->Fill(track->GetTPCmomentum(), track->GetTPCsignalN());
2719 return kFALSE;// Skip generated response in case of error
2723 for (Int_t n = 0; n < numGenEntries; n++) {
2724 if (!isMC || !fUseMCidForGeneration) {
2725 // If no MC info is available or shall not be used, use weighting with priors to generate entries for the different species
2726 Double_t rnd = fRandom->Rndm(); // Produce uniformly distributed floating point in ]0, 1]
2728 // Consider generated response as originating from...
2729 if (rnd <= prob[AliPID::kElectron])
2730 genEntry[kGenMCID] = 0; // ... an electron
2731 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon])
2732 genEntry[kGenMCID] = 1; // ... a kaon
2733 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon])
2734 genEntry[kGenMCID] = 2; // ... a muon -> NOTE: prob[AliPID::kMuon] = 0 in case of fTakeIntoAccountMuons = kFALSE
2735 else if (rnd <= prob[AliPID::kElectron] + prob[AliPID::kKaon] + prob[AliPID::kMuon] + prob[AliPID::kPion])
2736 genEntry[kGenMCID] = 3; // ... a pion
2738 genEntry[kGenMCID] = 4; // ... a proton
2742 genEntry[kGenSelectSpecies] = 0;
2743 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeEl[n];
2744 fhGenEl->Fill(genEntry);
2746 genEntry[kGenSelectSpecies] = 1;
2747 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimeKa[n];
2748 fhGenEl->Fill(genEntry);
2750 genEntry[kGenSelectSpecies] = 2;
2751 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePi[n];
2752 fhGenEl->Fill(genEntry);
2754 genEntry[kGenSelectSpecies] = 3;
2755 genEntry[kGenDeltaPrimeSpecies] = fGenRespElDeltaPrimePr[n];
2756 fhGenEl->Fill(genEntry);
2759 genEntry[kGenSelectSpecies] = 0;
2760 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeEl[n];
2761 fhGenKa->Fill(genEntry);
2763 genEntry[kGenSelectSpecies] = 1;
2764 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimeKa[n];
2765 fhGenKa->Fill(genEntry);
2767 genEntry[kGenSelectSpecies] = 2;
2768 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePi[n];
2769 fhGenKa->Fill(genEntry);
2771 genEntry[kGenSelectSpecies] = 3;
2772 genEntry[kGenDeltaPrimeSpecies] = fGenRespKaDeltaPrimePr[n];
2773 fhGenKa->Fill(genEntry);
2776 genEntry[kGenSelectSpecies] = 0;
2777 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeEl[n];
2778 fhGenPi->Fill(genEntry);
2780 genEntry[kGenSelectSpecies] = 1;
2781 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimeKa[n];
2782 fhGenPi->Fill(genEntry);
2784 genEntry[kGenSelectSpecies] = 2;
2785 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePi[n];
2786 fhGenPi->Fill(genEntry);
2788 genEntry[kGenSelectSpecies] = 3;
2789 genEntry[kGenDeltaPrimeSpecies] = fGenRespPiDeltaPrimePr[n];
2790 fhGenPi->Fill(genEntry);
2792 if (fTakeIntoAccountMuons) {
2793 // Muons, if desired
2794 genEntry[kGenSelectSpecies] = 0;
2795 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeEl[n];
2796 fhGenMu->Fill(genEntry);
2798 genEntry[kGenSelectSpecies] = 1;
2799 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimeKa[n];
2800 fhGenMu->Fill(genEntry);
2802 genEntry[kGenSelectSpecies] = 2;
2803 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePi[n];
2804 fhGenMu->Fill(genEntry);
2806 genEntry[kGenSelectSpecies] = 3;
2807 genEntry[kGenDeltaPrimeSpecies] = fGenRespMuDeltaPrimePr[n];
2808 fhGenMu->Fill(genEntry);
2812 genEntry[kGenSelectSpecies] = 0;
2813 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeEl[n];
2814 fhGenPr->Fill(genEntry);
2816 genEntry[kGenSelectSpecies] = 1;
2817 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimeKa[n];
2818 fhGenPr->Fill(genEntry);
2820 genEntry[kGenSelectSpecies] = 2;
2821 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePi[n];
2822 fhGenPr->Fill(genEntry);
2824 genEntry[kGenSelectSpecies] = 3;
2825 genEntry[kGenDeltaPrimeSpecies] = fGenRespPrDeltaPrimePr[n];
2826 fhGenPr->Fill(genEntry);
2830 printf("File: %s, Line: %d: ProcessTrack -> Done\n", (char*)__FILE__, __LINE__);
2836 //_____________________________________________________________________________
2837 Bool_t AliAnalysisTaskPID::SetConvolutedGaussLambdaParameter(Double_t lambda)
2839 // Set the lambda parameter of the convolution to the desired value. Automatically
2840 // calculates the parameters for the transition (restricted) gauss -> convoluted gauss.
2841 fConvolutedGaussTransitionPars[2] = lambda;
2843 // Save old parameters and settings of function to restore them later:
2844 Double_t* oldFuncParams = new Double_t[fkConvolutedGausNPar];
2845 for (Int_t i = 0; i < fkConvolutedGausNPar; i++)
2846 oldFuncParams[i] = fConvolutedGausDeltaPrime->GetParameter(i);
2847 Int_t oldFuncNpx = fConvolutedGausDeltaPrime->GetNpx();
2848 Double_t oldFuncRangeLow = 0, oldFuncRangeUp = 100;
2849 fConvolutedGausDeltaPrime->GetRange(oldFuncRangeLow, oldFuncRangeUp);
2851 // Choose some sufficiently large range
2852 const Double_t rangeStart = 0.5;
2853 const Double_t rangeEnd = 2.0;
2855 // To get the parameters for the transition, just choose arbitrarily input parameters for mu and sigma
2856 // (it makes sense to choose typical values). The ratio sigma_gauss / sigma_convolution is independent
2857 // of mu and as well the difference mu_gauss - mu_convolution.
2858 Double_t muInput = 1.0;
2859 Double_t sigmaInput = fgkSigmaReferenceForTransitionPars;
2862 // Step 1: Generate distribution with input parameters
2863 const Int_t nPar = 3;
2864 Double_t inputPar[nPar] = { muInput, sigmaInput, lambda };
2866 TH1D* hInput = new TH1D("hInput", "Input distribution", 2000, rangeStart, rangeEnd);
2868 fConvolutedGausDeltaPrime->SetParameters(inputPar);
2869 fConvolutedGausDeltaPrime->SetRange(rangeStart, rangeEnd);
2870 fConvolutedGausDeltaPrime->SetNpx(2000);
2873 // The resolution and mean of the AliPIDResponse are determined in finite intervals
2874 // of ncl (also second order effects due to finite dEdx and tanTheta).
2875 // Take this into account for the transition parameters via assuming an approximately flat
2876 // distritubtion in ncl in this interval.
2877 // NOTE: The ncl interval should be the same as the one used for the sigma map creation!
2878 const Int_t nclStart = 151;
2879 const Int_t nclEnd = 160;
2880 const Int_t nclSteps = (nclEnd - nclStart) + 1;
2881 for (Int_t ncl = nclStart; ncl <= nclEnd; ncl++) {
2882 // Resolution scales with 1/sqrt(ncl)
2883 fConvolutedGausDeltaPrime->SetParameter(1, inputPar[1] * sqrt(nclEnd) / sqrt(ncl));
2884 TH1* hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
2886 for (Int_t i = 0; i < 50000000 / nclSteps; i++)
2887 hInput->Fill(hProbDensity->GetRandom());
2891 TH1* hProbDensity = fConvolutedGausDeltaPrime->GetHistogram();
2893 for (Int_t i = 0; i < 50000000; i++)
2894 hInput->Fill(hProbDensity->GetRandom());
2896 // Step 2: Fit generated distribution with restricted gaussian
2897 Int_t maxBin = hInput->GetMaximumBin();
2898 Double_t max = hInput->GetBinContent(maxBin);
2900 UChar_t usedBins = 1;
2902 max += hInput->GetBinContent(maxBin - 1);
2905 if (maxBin < hInput->GetNbinsX()) {
2906 max += hInput->GetBinContent(maxBin + 1);
2911 // NOTE: The range (<-> fraction of maximum) should be the same
2912 // as for the spline and eta maps creation
2913 const Double_t lowThreshold = hInput->GetXaxis()->GetBinLowEdge(hInput->FindFirstBinAbove(0.1 * max));
2914 const Double_t highThreshold = hInput->GetXaxis()->GetBinUpEdge(hInput->FindLastBinAbove(0.1 * max));
2916 TFitResultPtr fitResGaussFirstStep = hInput->Fit("gaus", "QNRS", "", lowThreshold, highThreshold);
2918 TFitResultPtr fitResGauss;
2920 if ((Int_t)fitResGaussFirstStep == 0) {
2921 TF1 fGauss("fGauss", "[0]*TMath::Gaus(x, [1], [2], 1)", rangeStart, rangeEnd);
2922 fGauss.SetParameter(0, fitResGaussFirstStep->GetParams()[0]);
2923 fGauss.SetParError(0, fitResGaussFirstStep->GetErrors()[0]);
2924 fGauss.SetParameter(2, fitResGaussFirstStep->GetParams()[2]);
2925 fGauss.SetParError(2, fitResGaussFirstStep->GetErrors()[2]);
2927 fGauss.FixParameter(1, fitResGaussFirstStep->GetParams()[1]);
2928 fitResGauss = hInput->Fit(&fGauss, "QNS", "s", rangeStart, rangeEnd);
2931 fitResGauss = hInput->Fit("gaus", "QNRS", "same", rangeStart, rangeEnd);
2933 //OLD TFitResultPtr fitResGauss = hInput->Fit("gaus", "QNRS", "", hInput->GetXaxis()->GetBinLowEdge(hInput->FindFirstBinAbove(0.1 * max)),
2934 // hInput->GetXaxis()->GetBinUpEdge(hInput->FindLastBinAbove(0.1 * max)));
2937 // Step 3: Use parameters from gaussian fit to obtain parameters for the transition "restricted gauss" -> "convoluted gauss"
2939 // 3.1 The ratio sigmaInput / sigma_gaussFit ONLY depends on lambda (which is fixed per period) -> Calculate this first
2940 // for an arbitrary (here: typical) sigma. The ratio is then ~the same for ALL sigma for given lambda!
2941 if ((Int_t)fitResGauss != 0) {
2942 AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Gauss Fit failed!\n");
2943 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2944 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2945 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2948 delete [] oldFuncParams;
2953 if (fitResGauss->GetParams()[2] <= 0) {
2954 AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Sigma of gauss fit <= 0!\n");
2955 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
2956 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
2957 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
2960 delete [] oldFuncParams;
2965 // sigma correction factor
2966 fConvolutedGaussTransitionPars[1] = sigmaInput / fitResGauss->GetParams()[2];
2968 // 3.2 Now that sigma und lambda are determined, one can calculate mu by shifting the maximum to the desired position,
2969 // i.e. the maximum of the input distribution should coincide with that of the re-calculated distribution,
2970 // which is achieved by getting the same mu for the same sigma.
2971 // NOTE: For fixed lambda, the shift is proportional to sigma and independent of mu!
2972 // So, one can calculate the shift for an arbitrary fixed (here: typical)
2973 // sigma and then simply use this shift for any other sigma by scaling it correspondingly!!!
2975 // Mu shift correction:
2976 // Shift in mu (difference between mean of gauss and mean of convolution) is proportional to sigma!
2977 // Thus, choose a reference sigma (typical -> 0.05), such that for arbitrary sigma one can simple scale
2978 // this shift correction with sigma / referenceSigma.
2979 fConvolutedGaussTransitionPars[0] = (fitResGauss->GetParams()[1] - muInput);
2982 /*Changed on 03.07.2013
2984 // Maximum of fConvolutedGausDeltaPrime should agree with maximum of input
2985 Double_t par[nPar] = { fitResGauss->GetParams()[1], // just as a guess of the maximum position
2989 fConvolutedGausDeltaPrime->SetParameters(par);
2991 Double_t maxXInput = fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001, muInput - 3. * sigmaInput),
2992 muInput + 10. * sigmaInput,
2995 // Maximum shifts proportional to sigma and is linear in mu (~mean of gauss)
2996 // Maximum should be typically located within [gaussMean, gaussMean + 3 gaussSigma].
2997 // -> Larger search range for safety reasons (also: sigma and/or mean might be not 100% accurate).
2998 Double_t maxXconvoluted = fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001,
2999 fitResGauss->GetParams()[1] - 3. * fitResGauss->GetParams()[2]),
3000 fitResGauss->GetParams()[1] + 10. * fitResGauss->GetParams()[2],
3002 if (maxXconvoluted <= 0) {
3003 AliError("Not able to calculate parameters for the transition \"restricted gauss\" -> \"convoluted gauss\": Maximum of fConvolutedGausDeltaPrime <= 0!\n");
3005 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
3006 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
3007 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
3010 delete [] oldFuncParams;
3015 // maxX perfectly shifts as par[0] (scaled by sigma) -> Can shift maxX to input value.
3016 // Mu shift correction:
3017 fConvolutedGaussTransitionPars[0] = maxXconvoluted - maxXInput;
3022 fConvolutedGausDeltaPrime->SetParameters(oldFuncParams);
3023 fConvolutedGausDeltaPrime->SetNpx(oldFuncNpx);
3024 fConvolutedGausDeltaPrime->SetRange(oldFuncRangeLow, oldFuncRangeUp);
3027 delete [] oldFuncParams;
3033 //_____________________________________________________________________________
3034 AliAnalysisTaskPID::ErrorCode AliAnalysisTaskPID::SetParamsForConvolutedGaus(Double_t gausMean, Double_t gausSigma)
3036 // Set parameters for convoluted gauss using parameters for a pure gaussian.
3037 // If SetConvolutedGaussLambdaParameter has not been called before to initialise the translation parameters,
3038 // some default parameters will be used and an error will show up.
3041 printf("File: %s, Line: %d: SetParamsForConvolutedGaus: mean %e, sigma %e\n", (char*)__FILE__, __LINE__, gausMean, gausSigma);
3043 if (fConvolutedGaussTransitionPars[1] < -998) {
3044 AliError("Transition parameters not initialised! Default parameters will be used. Please call SetConvolutedGaussLambdaParameter(...) before any calculations!");
3045 SetConvolutedGaussLambdaParameter(2.0);
3046 AliError(Form("Parameters set to:\n[0]: %f\n[1]: %f\n[2]: %f\n", fConvolutedGaussTransitionPars[0],
3047 fConvolutedGaussTransitionPars[1], fConvolutedGaussTransitionPars[2]));
3050 Double_t par[fkConvolutedGausNPar];
3051 par[2] = fConvolutedGaussTransitionPars[2];
3052 par[1] = fConvolutedGaussTransitionPars[1] * gausSigma;
3053 // maxX perfectly shifts as par[0] (scaled by sigma) -> Can shift maxX so that it sits at the right place.
3054 par[0] = gausMean - fConvolutedGaussTransitionPars[0] * par[1] / fgkSigmaReferenceForTransitionPars;
3056 ErrorCode errCode = kNoErrors;
3057 fConvolutedGausDeltaPrime->SetParameters(par);
3060 printf("File: %s, Line: %d: SetParamsForConvolutedGaus -> Parameters set to: %e, %e, %e (transition pars: %e, %e, %e, %e)\n",
3061 (char*)__FILE__, __LINE__, par[0], par[1], par[2], fConvolutedGaussTransitionPars[0], fConvolutedGaussTransitionPars[1],
3062 fConvolutedGaussTransitionPars[2], fgkSigmaReferenceForTransitionPars);
3064 fConvolutedGausDeltaPrime->SetNpx(20); // Small value speeds up following algorithm (valid, since extrema far apart)
3066 // Accuracy of 10^-5 is enough to get 0.1% precise peak for MIPS w.r.t. to dEdx = 2000 of protons
3067 // (should boost up the algorithm, because 10^-10 is the default value!)
3068 Double_t maxX= fConvolutedGausDeltaPrime->GetMaximumX(TMath::Max(0.001, gausMean - 2. * gausSigma),
3069 gausMean + 6. * gausSigma, 1.0E-5);
3071 const Double_t maximum = fConvolutedGausDeltaPrime->Eval(maxX);
3072 const Double_t maximumFraction = maximum * fAccuracyNonGaussianTail;
3074 // Estimate lower boundary for subsequent search:
3075 Double_t lowBoundSearchBoundLow = TMath::Max(1e-4, maxX - 5. * gausSigma);
3076 Double_t lowBoundSearchBoundUp = maxX;
3078 Bool_t lowerBoundaryFixedAtZero = kFALSE;
3080 while (fConvolutedGausDeltaPrime->Eval(lowBoundSearchBoundLow) >= maximumFraction) {
3081 if (lowBoundSearchBoundLow <= 0) {
3082 // This should only happen to low dEdx particles with very few clusters and therefore large sigma, such that the gauss goes below zero deltaPrime
3083 if (maximum <= 0) { // Something is weired
3084 printf("Error generating signal: maximum is <= 0!\n");
3088 const Double_t valueAtZero = fConvolutedGausDeltaPrime->Eval(0);
3089 if (valueAtZero / maximum > 0.05) {
3090 // Too large fraction below zero deltaPrime. Signal generation cannot be reliable in this case
3091 printf("Error generating signal: Too large fraction below zero deltaPrime: convGauss(0) / convGauss(max) = %e / %e = %e!\n",
3092 valueAtZero, maximum, valueAtZero / maximum);
3098 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",
3099 fConvolutedGausDeltaPrime->Eval(maxX), fAccuracyNonGaussianTail, maximumFraction, maxX, par[0], par[1], par[2], gausMean, gausSigma);
3102 lowerBoundaryFixedAtZero = kTRUE;
3104 if (errCode != kError)
3110 lowBoundSearchBoundUp -= gausSigma;
3111 lowBoundSearchBoundLow -= gausSigma;
3113 if (lowBoundSearchBoundLow < 0) {
3114 lowBoundSearchBoundLow = 0;
3115 lowBoundSearchBoundUp += gausSigma;
3119 // Determine lower boundary inside estimated range. For small values of the maximum: Need more precision, since finer binning!
3120 Double_t rangeStart = lowerBoundaryFixedAtZero ? 0 :
3121 fConvolutedGausDeltaPrime->GetX(maximumFraction, lowBoundSearchBoundLow, lowBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
3123 // .. and the same for the upper boundary
3124 Double_t rangeEnd = 0;
3125 // If distribution starts beyond upper boundary, everything ends up in the overflow bin. So, just reduce range and Npx to minimum
3126 if (rangeStart > fkDeltaPrimeUpLimit) {
3127 rangeEnd = rangeStart + 0.00001;
3128 fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
3129 fConvolutedGausDeltaPrime->SetNpx(4);
3132 // Estimate upper boundary for subsequent search:
3133 Double_t upBoundSearchBoundUp = maxX + 5 * gausSigma;
3134 Double_t upBoundSearchBoundLow = maxX;
3135 while (fConvolutedGausDeltaPrime->Eval(upBoundSearchBoundUp) >= maximumFraction) {
3136 upBoundSearchBoundUp += gausSigma;
3137 upBoundSearchBoundLow += gausSigma;
3140 // For small values of the maximum: Need more precision, since finer binning!
3141 rangeEnd = fConvolutedGausDeltaPrime->GetX(maximumFraction, upBoundSearchBoundLow, upBoundSearchBoundUp, (maxX < 0.4) ? 1e-5 : 0.001);
3143 fConvolutedGausDeltaPrime->SetRange(rangeStart,rangeEnd);
3144 fConvolutedGausDeltaPrime->SetNpx(fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindBin(rangeEnd)
3145 - fhPIDdataAll->GetAxis(kDataDeltaPrimeSpecies)->FindBin(rangeStart) + 1);
3146 //fConvolutedGausDeltaPrime->SetNpx((rangeEnd - rangeStart) / fDeltaPrimeBinWidth + 1);
3150 printf("File: %s, Line: %d: SetParamsForConvolutedGaus -> range %f - %f, error code %d\n", (char*)__FILE__, __LINE__,
3151 rangeStart, rangeEnd, errCode);
3157 //________________________________________________________________________
3158 void AliAnalysisTaskPID::SetUpGenHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const
3160 // Sets bin limits for axes which are not standard binned and the axes titles.
3162 hist->SetBinEdges(kGenPt, binsPt);
3163 hist->SetBinEdges(kGenDeltaPrimeSpecies, binsDeltaPrime);
3164 hist->SetBinEdges(kGenCentrality, binsCent);
3166 if (fStoreAdditionalJetInformation)
3167 hist->SetBinEdges(kGenJetPt, binsJetPt);
3170 hist->GetAxis(kGenMCID)->SetTitle("MC PID");
3171 hist->GetAxis(kGenMCID)->SetBinLabel(1, "e");
3172 hist->GetAxis(kGenMCID)->SetBinLabel(2, "K");
3173 hist->GetAxis(kGenMCID)->SetBinLabel(3, "#mu");
3174 hist->GetAxis(kGenMCID)->SetBinLabel(4, "#pi");
3175 hist->GetAxis(kGenMCID)->SetBinLabel(5, "p");
3177 hist->GetAxis(kGenSelectSpecies)->SetTitle("Select Species");
3178 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(1, "e");
3179 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(2, "K");
3180 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(3, "#pi");
3181 hist->GetAxis(kGenSelectSpecies)->SetBinLabel(4, "p");
3183 hist->GetAxis(kGenPt)->SetTitle("p_{T} (GeV/c)");
3185 hist->GetAxis(kGenDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
3187 hist->GetAxis(kGenCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
3189 if (fStoreAdditionalJetInformation) {
3190 hist->GetAxis(kGenJetPt)->SetTitle("p_{T}^{jet} (GeV/c)");
3192 hist->GetAxis(kGenZ)->SetTitle("z = p_{T}^{track} / p_{T}^{jet}");
3194 hist->GetAxis(kGenXi)->SetTitle("#xi = ln(p_{T}^{jet} / p_{T}^{track})");
3197 hist->GetAxis(GetIndexOfChargeAxisGen())->SetTitle("Charge (e_{0})");
3199 hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetTitle("TOF PID Info");
3200 // Offset is (TMath::Abs(kNoTOFinfo) + 1), such that bin 1 (= first label) corresponds to kNoTOFinfo (< 0)
3201 hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kNoTOFinfo + (TMath::Abs(kNoTOFinfo) + 1), "No TOF Info");
3202 hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kNoTOFpid + (TMath::Abs(kNoTOFinfo) + 1), "No TOF PID");
3203 hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kTOFpion + (TMath::Abs(kNoTOFinfo) + 1), "#pi");
3204 hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kTOFkaon + (TMath::Abs(kNoTOFinfo) + 1), "K");
3205 hist->GetAxis(GetIndexOfTOFpidInfoAxisGen())->SetBinLabel(kTOFproton + (TMath::Abs(kNoTOFinfo) + 1), "p");
3209 //________________________________________________________________________
3210 void AliAnalysisTaskPID::SetUpGenYieldHist(THnSparse* hist, Double_t* binsPt, Double_t* binsCent, Double_t* binsJetPt) const
3212 // Sets bin limits for axes which are not standard binned and the axes titles.
3214 hist->SetBinEdges(kGenYieldPt, binsPt);
3215 hist->SetBinEdges(kGenYieldCentrality, binsCent);
3216 if (fStoreAdditionalJetInformation)
3217 hist->SetBinEdges(kGenYieldJetPt, binsJetPt);
3219 for (Int_t i = 0; i < 5; i++)
3220 hist->GetAxis(kGenYieldMCID)->SetBinLabel(i + 1, AliPID::ParticleLatexName(i));
3223 hist->GetAxis(kGenYieldMCID)->SetTitle("MC PID");
3224 hist->GetAxis(kGenYieldPt)->SetTitle("p_{T}^{gen} (GeV/c)");
3225 hist->GetAxis(kGenYieldCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
3227 if (fStoreAdditionalJetInformation) {
3228 hist->GetAxis(kGenYieldJetPt)->SetTitle("p_{T}^{jet, gen} (GeV/c)");
3230 hist->GetAxis(kGenYieldZ)->SetTitle("z = p_{T}^{track} / p_{T}^{jet}");
3232 hist->GetAxis(kGenYieldXi)->SetTitle("#xi = ln(p_{T}^{jet} / p_{T}^{track})");
3235 hist->GetAxis(GetIndexOfChargeAxisGenYield())->SetTitle("Charge (e_{0})");
3239 //________________________________________________________________________
3240 void AliAnalysisTaskPID::SetUpHist(THnSparse* hist, Double_t* binsPt, Double_t* binsDeltaPrime, Double_t* binsCent, Double_t* binsJetPt) const
3242 // Sets bin limits for axes which are not standard binned and the axes titles.
3244 hist->SetBinEdges(kDataPt, binsPt);
3245 hist->SetBinEdges(kDataDeltaPrimeSpecies, binsDeltaPrime);
3246 hist->SetBinEdges(kDataCentrality, binsCent);
3248 if (fStoreAdditionalJetInformation)
3249 hist->SetBinEdges(kDataJetPt, binsJetPt);
3252 hist->GetAxis(kDataMCID)->SetTitle("MC PID");
3253 hist->GetAxis(kDataMCID)->SetBinLabel(1, "e");
3254 hist->GetAxis(kDataMCID)->SetBinLabel(2, "K");
3255 hist->GetAxis(kDataMCID)->SetBinLabel(3, "#mu");
3256 hist->GetAxis(kDataMCID)->SetBinLabel(4, "#pi");
3257 hist->GetAxis(kDataMCID)->SetBinLabel(5, "p");
3259 hist->GetAxis(kDataSelectSpecies)->SetTitle("Select Species");
3260 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(1, "e");
3261 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(2, "K");
3262 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(3, "#pi");
3263 hist->GetAxis(kDataSelectSpecies)->SetBinLabel(4, "p");
3265 hist->GetAxis(kDataPt)->SetTitle("p_{T} (GeV/c)");
3267 hist->GetAxis(kDataDeltaPrimeSpecies)->SetTitle("TPC #Delta'_{species} (arb. unit)");
3269 hist->GetAxis(kDataCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));
3271 if (fStoreAdditionalJetInformation) {
3272 hist->GetAxis(kDataJetPt)->SetTitle("p_{T}^{jet} (GeV/c)");
3274 hist->GetAxis(kDataZ)->SetTitle("z = p_{T}^{track} / p_{T}^{jet}");
3276 hist->GetAxis(kDataXi)->SetTitle("#xi = ln(p_{T}^{jet} / p_{T}^{track})");
3279 hist->GetAxis(GetIndexOfChargeAxisData())->SetTitle("Charge (e_{0})");
3281 hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetTitle("TOF PID Info");
3282 // Offset is (TMath::Abs(kNoTOFinfo) + 1), such that bin 1 (= first label) corresponds to kNoTOFinfo (< 0)
3283 hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kNoTOFinfo + (TMath::Abs(kNoTOFinfo) + 1), "No TOF Info");
3284 hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kNoTOFpid + (TMath::Abs(kNoTOFinfo) + 1), "No TOF PID");
3285 hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kTOFpion + (TMath::Abs(kNoTOFinfo) + 1), "#pi");
3286 hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kTOFkaon + (TMath::Abs(kNoTOFinfo) + 1), "K");
3287 hist->GetAxis(GetIndexOfTOFpidInfoAxisData())->SetBinLabel(kTOFproton + (TMath::Abs(kNoTOFinfo) + 1), "p");
3291 //________________________________________________________________________
3292 void AliAnalysisTaskPID::SetUpPtResHist(THnSparse* hist, Double_t* binsPt, Double_t* binsJetPt, Double_t* binsCent) const
3294 // Sets bin limits for axes which are not standard binned and the axes titles.
3296 hist->SetBinEdges(kPtResJetPt, binsJetPt);
3297 hist->SetBinEdges(kPtResGenPt, binsPt);
3298 hist->SetBinEdges(kPtResRecPt, binsPt);
3299 hist->SetBinEdges(kPtResCentrality, binsCent);
3302 hist->GetAxis(kPtResJetPt)->SetTitle("p_{T}^{jet, rec} (GeV/c)");
3303 hist->GetAxis(kPtResGenPt)->SetTitle("p_{T}^{gen} (GeV/c)");
3304 hist->GetAxis(kPtResRecPt)->SetTitle("p_{T}^{rec} (GeV/c)");
3306 hist->GetAxis(kPtResCharge)->SetTitle("Charge (e_{0})");
3307 hist->GetAxis(kPtResCentrality)->SetTitle(Form("Centrality Percentile (%s)", fCentralityEstimator.Data()));