]>
Commit | Line | Data |
---|---|---|
10d100d4 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | //----------------------------------------------------------------- | |
17 | // Implementation of the TPC PID class | |
18 | // Very naive one... Should be made better by the detector experts... | |
19 | // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch | |
20 | // With many additions and modifications suggested by | |
21 | // Alexander Kalweit, GSI, alexander.philipp.kalweit@cern.ch | |
22 | // Dariusz Miskowiec, GSI, D.Miskowiec@gsi.de | |
644666df | 23 | // ...and some modifications by |
24 | // Mikolaj Krzewicki, GSI, mikolaj.krzewicki@cern.ch | |
f84b18dd | 25 | // ...and some modifications plus eta correction functions by |
26 | // Benjamin Hess, University of Tuebingen, bhess@cern.ch | |
10d100d4 | 27 | //----------------------------------------------------------------- |
28 | ||
d2aa6df0 | 29 | #include <TGraph.h> |
30 | #include <TObjArray.h> | |
31 | #include <TSpline.h> | |
644666df | 32 | #include <TBits.h> |
00a38d07 | 33 | #include <TMath.h> |
f84b18dd | 34 | #include <TH2D.h> |
d2aa6df0 | 35 | |
644666df | 36 | #include <AliLog.h> |
10d100d4 | 37 | #include "AliExternalTrackParam.h" |
644666df | 38 | #include "AliVTrack.h" |
d2aa6df0 | 39 | #include "AliTPCPIDResponse.h" |
644666df | 40 | #include "AliTPCdEdxInfo.h" |
912f7fdc | 41 | #include "TFile.h" |
42 | #include "TSpline.h" | |
644666df | 43 | |
44 | ClassImp(AliTPCPIDResponse); | |
d2aa6df0 | 45 | |
912f7fdc | 46 | |
47 | AliTPCPIDResponse *AliTPCPIDResponse::fgInstance =0; | |
48 | ||
644666df | 49 | const char* AliTPCPIDResponse::fgkGainScenarioName[fgkNumberOfGainScenarios+1]= |
50 | { | |
51 | "", //default - no name | |
52 | "1", //all high | |
53 | "2", //OROC only | |
54 | "unknown" | |
55 | }; | |
10d100d4 | 56 | |
57 | //_________________________________________________________________________ | |
58 | AliTPCPIDResponse::AliTPCPIDResponse(): | |
644666df | 59 | TNamed(), |
d2aa6df0 | 60 | fMIP(50.), |
644666df | 61 | fRes0(), |
62 | fResN2(), | |
d2aa6df0 | 63 | fKp1(0.0283086), |
64 | fKp2(2.63394e+01), | |
65 | fKp3(5.04114e-11), | |
66 | fKp4(2.12543), | |
67 | fKp5(4.88663), | |
68 | fUseDatabase(kFALSE), | |
644666df | 69 | fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios), |
70 | fVoltageMap(72), | |
71 | fLowGainIROCthreshold(-40), | |
72 | fBadIROCthreshhold(-70), | |
73 | fLowGainOROCthreshold(-40), | |
74 | fBadOROCthreshhold(-40), | |
75 | fMaxBadLengthFraction(0.5), | |
f84b18dd | 76 | fMagField(0.), |
77 | fhEtaCorr(0x0), | |
78 | fhEtaSigmaPar1(0x0), | |
87da0205 | 79 | fSigmaPar0(0.0), |
80 | fCurrentEventMultiplicity(0), | |
81 | fCorrFuncMultiplicity(0x0), | |
82 | fCorrFuncMultiplicityTanTheta(0x0), | |
03696add | 83 | fCorrFuncSigmaMultiplicity(0x0), |
84 | fSplineArray() | |
10d100d4 | 85 | { |
86 | // | |
87 | // The default constructor | |
88 | // | |
5a9dc560 | 89 | |
90 | AliLog::SetClassDebugLevel("AliTPCPIDResponse", AliLog::kInfo); | |
91 | ||
644666df | 92 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=0.07;fResN2[i]=0.0;} |
87da0205 | 93 | |
94 | fCorrFuncMultiplicity = new TF1("fCorrFuncMultiplicity", | |
95 | "[0] + [1]*TMath::Max([4], TMath::Min(x, [3])) + [2] * TMath::Power(TMath::Max([4], TMath::Min(x, [3])), 2)", | |
96 | 0., 0.2); | |
97 | fCorrFuncMultiplicityTanTheta = new TF1("fCorrFuncMultiplicityTanTheta", "[0] * (x - [2]) + [1] * (x * x - [2] * [2])", -1.5, 1.5); | |
98 | fCorrFuncSigmaMultiplicity = new TF1("fCorrFuncSigmaMultiplicity", | |
99 | "TMath::Max(0, [0] + [1]*TMath::Min(x, [3]) + [2] * TMath::Power(TMath::Min(x, [3]), 2))", 0., 0.2); | |
100 | ||
101 | ||
102 | ResetMultiplicityCorrectionFunctions(); | |
912f7fdc | 103 | fgInstance=this; |
10d100d4 | 104 | } |
1b45e564 | 105 | /*TODO remove? |
10d100d4 | 106 | //_________________________________________________________________________ |
d2aa6df0 | 107 | AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param): |
644666df | 108 | TNamed(), |
d2aa6df0 | 109 | fMIP(param[0]), |
644666df | 110 | fRes0(), |
111 | fResN2(), | |
d2aa6df0 | 112 | fKp1(0.0283086), |
113 | fKp2(2.63394e+01), | |
114 | fKp3(5.04114e-11), | |
115 | fKp4(2.12543), | |
116 | fKp5(4.88663), | |
117 | fUseDatabase(kFALSE), | |
644666df | 118 | fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios), |
119 | fVoltageMap(72), | |
120 | fLowGainIROCthreshold(-40), | |
121 | fBadIROCthreshhold(-70), | |
122 | fLowGainOROCthreshold(-40), | |
123 | fBadOROCthreshhold(-40), | |
124 | fMaxBadLengthFraction(0.5), | |
f84b18dd | 125 | fMagField(0.), |
126 | fhEtaCorr(0x0), | |
127 | fhEtaSigmaPar1(0x0), | |
128 | fSigmaPar0(0.0) | |
10d100d4 | 129 | { |
130 | // | |
131 | // The main constructor | |
132 | // | |
644666df | 133 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=param[1];fResN2[i]=param[2];} |
134 | } | |
1b45e564 | 135 | */ |
f84b18dd | 136 | |
137 | //_________________________________________________________________________ | |
138 | AliTPCPIDResponse::~AliTPCPIDResponse() | |
139 | { | |
140 | // | |
141 | // Destructor | |
142 | // | |
143 | ||
144 | delete fhEtaCorr; | |
145 | fhEtaCorr = 0x0; | |
146 | ||
147 | delete fhEtaSigmaPar1; | |
148 | fhEtaSigmaPar1 = 0x0; | |
87da0205 | 149 | |
150 | delete fCorrFuncMultiplicity; | |
151 | fCorrFuncMultiplicity = 0x0; | |
152 | ||
153 | delete fCorrFuncMultiplicityTanTheta; | |
154 | fCorrFuncMultiplicityTanTheta = 0x0; | |
155 | ||
156 | delete fCorrFuncSigmaMultiplicity; | |
157 | fCorrFuncSigmaMultiplicity = 0x0; | |
912f7fdc | 158 | if (fgInstance==this) fgInstance=0; |
f84b18dd | 159 | } |
160 | ||
161 | ||
644666df | 162 | //_________________________________________________________________________ |
163 | AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& that): | |
164 | TNamed(that), | |
165 | fMIP(that.fMIP), | |
166 | fRes0(), | |
167 | fResN2(), | |
168 | fKp1(that.fKp1), | |
169 | fKp2(that.fKp2), | |
170 | fKp3(that.fKp3), | |
171 | fKp4(that.fKp4), | |
172 | fKp5(that.fKp5), | |
173 | fUseDatabase(that.fUseDatabase), | |
174 | fResponseFunctions(that.fResponseFunctions), | |
175 | fVoltageMap(that.fVoltageMap), | |
176 | fLowGainIROCthreshold(that.fLowGainIROCthreshold), | |
177 | fBadIROCthreshhold(that.fBadIROCthreshhold), | |
178 | fLowGainOROCthreshold(that.fLowGainOROCthreshold), | |
179 | fBadOROCthreshhold(that.fBadOROCthreshhold), | |
180 | fMaxBadLengthFraction(that.fMaxBadLengthFraction), | |
f84b18dd | 181 | fMagField(that.fMagField), |
182 | fhEtaCorr(0x0), | |
183 | fhEtaSigmaPar1(0x0), | |
87da0205 | 184 | fSigmaPar0(that.fSigmaPar0), |
185 | fCurrentEventMultiplicity(that.fCurrentEventMultiplicity), | |
186 | fCorrFuncMultiplicity(0x0), | |
187 | fCorrFuncMultiplicityTanTheta(0x0), | |
03696add | 188 | fCorrFuncSigmaMultiplicity(0x0), |
189 | fSplineArray() | |
644666df | 190 | { |
191 | //copy ctor | |
192 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];} | |
f84b18dd | 193 | |
194 | // Copy eta maps | |
1b45e564 | 195 | if (that.fhEtaCorr) { |
3800f96f | 196 | fhEtaCorr = new TH2D(*(that.fhEtaCorr)); |
197 | fhEtaCorr->SetDirectory(0); | |
198 | } | |
f84b18dd | 199 | |
1b45e564 | 200 | if (that.fhEtaSigmaPar1) { |
3800f96f | 201 | fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1)); |
202 | fhEtaSigmaPar1->SetDirectory(0); | |
203 | } | |
87da0205 | 204 | |
205 | // Copy multiplicity correction functions | |
206 | if (that.fCorrFuncMultiplicity) { | |
207 | fCorrFuncMultiplicity = new TF1(*(that.fCorrFuncMultiplicity)); | |
208 | } | |
209 | ||
210 | if (that.fCorrFuncMultiplicityTanTheta) { | |
211 | fCorrFuncMultiplicityTanTheta = new TF1(*(that.fCorrFuncMultiplicityTanTheta)); | |
212 | } | |
213 | ||
214 | if (that.fCorrFuncSigmaMultiplicity) { | |
215 | fCorrFuncSigmaMultiplicity = new TF1(*(that.fCorrFuncSigmaMultiplicity)); | |
216 | } | |
644666df | 217 | } |
218 | ||
219 | //_________________________________________________________________________ | |
220 | AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that) | |
221 | { | |
222 | //assignment | |
223 | if (&that==this) return *this; | |
224 | TNamed::operator=(that); | |
225 | fMIP=that.fMIP; | |
226 | fKp1=that.fKp1; | |
227 | fKp2=that.fKp2; | |
228 | fKp3=that.fKp3; | |
229 | fKp4=that.fKp4; | |
230 | fKp5=that.fKp5; | |
231 | fUseDatabase=that.fUseDatabase; | |
232 | fResponseFunctions=that.fResponseFunctions; | |
233 | fVoltageMap=that.fVoltageMap; | |
234 | fLowGainIROCthreshold=that.fLowGainIROCthreshold; | |
235 | fBadIROCthreshhold=that.fBadIROCthreshhold; | |
236 | fLowGainOROCthreshold=that.fLowGainOROCthreshold; | |
237 | fBadOROCthreshhold=that.fBadOROCthreshhold; | |
238 | fMaxBadLengthFraction=that.fMaxBadLengthFraction; | |
239 | fMagField=that.fMagField; | |
87da0205 | 240 | fCurrentEventMultiplicity=that.fCurrentEventMultiplicity; |
644666df | 241 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];} |
f84b18dd | 242 | |
243 | delete fhEtaCorr; | |
f85a3764 | 244 | fhEtaCorr=0x0; |
1b45e564 | 245 | if (that.fhEtaCorr) { |
3800f96f | 246 | fhEtaCorr = new TH2D(*(that.fhEtaCorr)); |
247 | fhEtaCorr->SetDirectory(0); | |
248 | } | |
f84b18dd | 249 | |
250 | delete fhEtaSigmaPar1; | |
f85a3764 | 251 | fhEtaSigmaPar1=0x0; |
1b45e564 | 252 | if (that.fhEtaSigmaPar1) { |
3800f96f | 253 | fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1)); |
254 | fhEtaSigmaPar1->SetDirectory(0); | |
255 | } | |
f84b18dd | 256 | |
257 | fSigmaPar0 = that.fSigmaPar0; | |
87da0205 | 258 | |
259 | delete fCorrFuncMultiplicity; | |
260 | fCorrFuncMultiplicity = 0x0; | |
261 | if (that.fCorrFuncMultiplicity) { | |
262 | fCorrFuncMultiplicity = new TF1(*(that.fCorrFuncMultiplicity)); | |
263 | } | |
264 | ||
265 | delete fCorrFuncMultiplicityTanTheta; | |
266 | fCorrFuncMultiplicityTanTheta = 0x0; | |
267 | if (that.fCorrFuncMultiplicityTanTheta) { | |
268 | fCorrFuncMultiplicityTanTheta = new TF1(*(that.fCorrFuncMultiplicityTanTheta)); | |
269 | } | |
270 | ||
271 | delete fCorrFuncSigmaMultiplicity; | |
272 | fCorrFuncSigmaMultiplicity = 0x0; | |
273 | if (that.fCorrFuncSigmaMultiplicity) { | |
274 | fCorrFuncSigmaMultiplicity = new TF1(*(that.fCorrFuncSigmaMultiplicity)); | |
275 | } | |
f84b18dd | 276 | |
644666df | 277 | return *this; |
10d100d4 | 278 | } |
279 | ||
d2aa6df0 | 280 | //_________________________________________________________________________ |
10d100d4 | 281 | Double_t AliTPCPIDResponse::Bethe(Double_t betaGamma) const { |
282 | // | |
283 | // This is the Bethe-Bloch function normalised to 1 at the minimum | |
284 | // WARNING | |
285 | // Simulated and reconstructed Bethe-Bloch differs | |
286 | // Simulated curve is the dNprim/dx | |
287 | // Reconstructed is proportianal dNtot/dx | |
288 | // Temporary fix for production - Simple linear correction function | |
289 | // Future 2 Bethe Bloch formulas needed | |
290 | // 1. for simulation | |
291 | // 2. for reconstructed PID | |
292 | // | |
d2aa6df0 | 293 | |
294 | // const Float_t kmeanCorrection =0.1; | |
10d100d4 | 295 | Double_t bb= |
296 | AliExternalTrackParam::BetheBlochAleph(betaGamma,fKp1,fKp2,fKp3,fKp4,fKp5); | |
297 | return bb*fMIP; | |
298 | } | |
299 | ||
300 | //_________________________________________________________________________ | |
301 | void AliTPCPIDResponse::SetBetheBlochParameters(Double_t kp1, | |
302 | Double_t kp2, | |
303 | Double_t kp3, | |
304 | Double_t kp4, | |
305 | Double_t kp5) { | |
306 | // | |
307 | // Set the parameters of the ALEPH Bethe-Bloch formula | |
308 | // | |
309 | fKp1=kp1; | |
310 | fKp2=kp2; | |
311 | fKp3=kp3; | |
312 | fKp4=kp4; | |
313 | fKp5=kp5; | |
314 | } | |
644666df | 315 | |
10d100d4 | 316 | //_________________________________________________________________________ |
317 | void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2) { | |
318 | // | |
319 | // Set the relative resolution sigma_rel = res0 * sqrt(1+resN2/npoint) | |
320 | // | |
644666df | 321 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=res0;fResN2[i]=resN2;} |
10d100d4 | 322 | } |
323 | ||
324 | //_________________________________________________________________________ | |
205d83d8 | 325 | Double_t AliTPCPIDResponse::GetExpectedSignal(Float_t mom, |
d2aa6df0 | 326 | AliPID::EParticleType n) const { |
10d100d4 | 327 | // |
f84b18dd | 328 | // Deprecated function (for backward compatibility). Please use |
329 | // GetExpectedSignal(const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource, | |
87da0205 | 330 | // Bool_t correctEta, Bool_t correctMultiplicity); |
f84b18dd | 331 | // instead! |
332 | // | |
333 | // | |
10d100d4 | 334 | // Calculates the expected PID signal as the function of |
335 | // the information stored in the track, for the specified particle type | |
336 | // | |
337 | // At the moment, these signals are just the results of calling the | |
338 | // Bethe-Bloch formula. | |
339 | // This can be improved. By taking into account the number of | |
340 | // assigned clusters and/or the track dip angle, for example. | |
341 | // | |
342 | ||
f84b18dd | 343 | //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...) |
344 | // !!! Splines for light nuclei need to be normalised to this factor !!! | |
345 | const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(n),2.3); | |
346 | ||
539a5a59 | 347 | Double_t mass=AliPID::ParticleMassZ(n); |
f84b18dd | 348 | if (!fUseDatabase) return Bethe(mom/mass) * chargeFactor; |
d2aa6df0 | 349 | // |
644666df | 350 | const TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n); |
00a38d07 | 351 | |
f84b18dd | 352 | if (!responseFunction) return Bethe(mom/mass) * chargeFactor; |
353 | ||
00a38d07 | 354 | return fMIP*responseFunction->Eval(mom/mass)*chargeFactor; |
d2aa6df0 | 355 | |
10d100d4 | 356 | } |
357 | ||
358 | //_________________________________________________________________________ | |
205d83d8 | 359 | Double_t AliTPCPIDResponse::GetExpectedSigma(Float_t mom, |
360 | Int_t nPoints, | |
644666df | 361 | AliPID::EParticleType n) const { |
10d100d4 | 362 | // |
f84b18dd | 363 | // Deprecated function (for backward compatibility). Please use |
364 | // GetExpectedSigma(onst AliVTrack* track, AliPID::EParticleType species, | |
365 | // ETPCdEdxSource dedxSource, Bool_t correctEta) instead! | |
366 | // | |
367 | // | |
10d100d4 | 368 | // Calculates the expected sigma of the PID signal as the function of |
369 | // the information stored in the track, for the specified particle type | |
370 | // | |
644666df | 371 | |
372 | if (nPoints != 0) | |
373 | return GetExpectedSignal(mom,n)*fRes0[0]*sqrt(1. + fResN2[0]/nPoints); | |
374 | else | |
375 | return GetExpectedSignal(mom,n)*fRes0[0]; | |
376 | } | |
377 | ||
378 | ////////////////////////////////////////////////////NEW////////////////////////////// | |
379 | ||
380 | //_________________________________________________________________________ | |
381 | void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario) { | |
382 | // | |
383 | // Set the relative resolution sigma_rel = res0 * sqrt(1+resN2/npoint) | |
384 | // | |
bd46387d | 385 | if ((Int_t)gainScenario>=(Int_t)fgkNumberOfGainScenarios) return; //TODO: better handling! |
644666df | 386 | fRes0[gainScenario]=res0; |
387 | fResN2[gainScenario]=resN2; | |
388 | } | |
389 | ||
f84b18dd | 390 | |
644666df | 391 | //_________________________________________________________________________ |
f84b18dd | 392 | Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track, |
644666df | 393 | AliPID::EParticleType species, |
f84b18dd | 394 | Double_t /*dEdx*/, |
395 | const TSpline3* responseFunction, | |
87da0205 | 396 | Bool_t correctEta, |
397 | Bool_t correctMultiplicity) const | |
644666df | 398 | { |
399 | // Calculates the expected PID signal as the function of | |
f84b18dd | 400 | // the information stored in the track and the given parameters, |
401 | // for the specified particle type | |
644666df | 402 | // |
403 | // At the moment, these signals are just the results of calling the | |
87da0205 | 404 | // Bethe-Bloch formula plus, if desired, taking into account the eta dependence |
405 | // and the multiplicity dependence (for PbPb). | |
644666df | 406 | // This can be improved. By taking into account the number of |
407 | // assigned clusters and/or the track dip angle, for example. | |
10d100d4 | 408 | // |
409 | ||
f84b18dd | 410 | Double_t mom=track->GetTPCmomentum(); |
97d17530 | 411 | Double_t mass=AliPID::ParticleMassZ(species); |
644666df | 412 | |
f84b18dd | 413 | //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...) |
414 | // !!! Splines for light nuclei need to be normalised to this factor !!! | |
415 | const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(species),2.3); | |
416 | ||
417 | if (!responseFunction) | |
418 | return Bethe(mom/mass) * chargeFactor; | |
644666df | 419 | |
f84b18dd | 420 | Double_t dEdxSplines = fMIP*responseFunction->Eval(mom/mass) * chargeFactor; |
87da0205 | 421 | |
422 | if (!correctEta && !correctMultiplicity) | |
f84b18dd | 423 | return dEdxSplines; |
424 | ||
87da0205 | 425 | Double_t corrFactorEta = 1.0; |
426 | Double_t corrFactorMultiplicity = 1.0; | |
427 | ||
428 | if (correctEta) { | |
ef7661fd | 429 | corrFactorEta = GetEtaCorrectionFast(track, dEdxSplines); |
87da0205 | 430 | //TODO Alternatively take current track dEdx |
ef7661fd | 431 | //corrFactorEta = GetEtaCorrectionFast(track, dEdx); |
87da0205 | 432 | } |
433 | ||
434 | if (correctMultiplicity) | |
ef7661fd | 435 | corrFactorMultiplicity = GetMultiplicityCorrectionFast(track, dEdxSplines * corrFactorEta, fCurrentEventMultiplicity); |
87da0205 | 436 | |
437 | return dEdxSplines * corrFactorEta * corrFactorMultiplicity; | |
644666df | 438 | } |
439 | ||
f84b18dd | 440 | |
644666df | 441 | //_________________________________________________________________________ |
442 | Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track, | |
443 | AliPID::EParticleType species, | |
f84b18dd | 444 | ETPCdEdxSource dedxSource, |
87da0205 | 445 | Bool_t correctEta, |
446 | Bool_t correctMultiplicity) const | |
644666df | 447 | { |
448 | // Calculates the expected PID signal as the function of | |
449 | // the information stored in the track, for the specified particle type | |
450 | // | |
451 | // At the moment, these signals are just the results of calling the | |
87da0205 | 452 | // Bethe-Bloch formula plus, if desired, taking into account the eta dependence |
453 | // and the multiplicity dependence (for PbPb). | |
644666df | 454 | // This can be improved. By taking into account the number of |
455 | // assigned clusters and/or the track dip angle, for example. | |
456 | // | |
457 | ||
f84b18dd | 458 | if (!fUseDatabase) { |
459 | //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...) | |
460 | // !!! Splines for light nuclei need to be normalised to this factor !!! | |
461 | const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(species),2.3); | |
644666df | 462 | |
97d17530 | 463 | return Bethe(track->GetTPCmomentum() / AliPID::ParticleMassZ(species)) * chargeFactor; |
f84b18dd | 464 | } |
644666df | 465 | |
f84b18dd | 466 | Double_t dEdx = -1; |
467 | Int_t nPoints = -1; | |
468 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
469 | TSpline3* responseFunction = 0x0; | |
470 | ||
471 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) { | |
472 | // Something is wrong with the track -> Return obviously invalid value | |
473 | return -999; | |
474 | } | |
475 | ||
476 | // Charge factor already taken into account inside the following function call | |
87da0205 | 477 | return GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity); |
644666df | 478 | } |
479 | ||
480 | //_________________________________________________________________________ | |
481 | TSpline3* AliTPCPIDResponse::GetResponseFunction( AliPID::EParticleType type, | |
482 | AliTPCPIDResponse::ETPCgainScenario gainScenario ) const | |
483 | { | |
484 | //get response function | |
485 | return dynamic_cast<TSpline3*>(fResponseFunctions.At(ResponseFunctionIndex(type,gainScenario))); | |
486 | } | |
487 | ||
488 | //_________________________________________________________________________ | |
489 | TSpline3* AliTPCPIDResponse::GetResponseFunction( const AliVTrack* track, | |
490 | AliPID::EParticleType species, | |
f84b18dd | 491 | ETPCdEdxSource dedxSource) const |
644666df | 492 | { |
493 | //the splines are stored in an array, different scenarios | |
494 | ||
f84b18dd | 495 | Double_t dEdx = -1; |
496 | Int_t nPoints = -1; | |
497 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
498 | TSpline3* responseFunction = 0x0; | |
499 | ||
500 | if (ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
501 | return responseFunction; | |
644666df | 502 | |
503 | return NULL; | |
504 | } | |
505 | ||
506 | //_________________________________________________________________________ | |
507 | void AliTPCPIDResponse::ResetSplines() | |
508 | { | |
509 | //reset the array with splines | |
510 | for (Int_t i=0;i<fResponseFunctions.GetEntriesFast();i++) | |
511 | { | |
512 | fResponseFunctions.AddAt(NULL,i); | |
513 | } | |
514 | } | |
515 | //_________________________________________________________________________ | |
516 | Int_t AliTPCPIDResponse::ResponseFunctionIndex( AliPID::EParticleType species, | |
517 | ETPCgainScenario gainScenario ) const | |
518 | { | |
519 | //get the index in fResponseFunctions given type and scenario | |
520 | return Int_t(species)+Int_t(gainScenario)*fgkNumberOfParticleSpecies; | |
521 | } | |
522 | ||
523 | //_________________________________________________________________________ | |
524 | void AliTPCPIDResponse::SetResponseFunction( TObject* o, | |
525 | AliPID::EParticleType species, | |
526 | ETPCgainScenario gainScenario ) | |
527 | { | |
528 | fResponseFunctions.AddAtAndExpand(o,ResponseFunctionIndex(species,gainScenario)); | |
529 | } | |
530 | ||
f84b18dd | 531 | |
644666df | 532 | //_________________________________________________________________________ |
f84b18dd | 533 | Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, |
534 | AliPID::EParticleType species, | |
535 | ETPCgainScenario gainScenario, | |
536 | Double_t dEdx, | |
537 | Int_t nPoints, | |
538 | const TSpline3* responseFunction, | |
87da0205 | 539 | Bool_t correctEta, |
540 | Bool_t correctMultiplicity) const | |
644666df | 541 | { |
542 | // Calculates the expected sigma of the PID signal as the function of | |
f84b18dd | 543 | // the information stored in the track and the given parameters, |
544 | // for the specified particle type | |
644666df | 545 | // |
546 | ||
f84b18dd | 547 | if (!responseFunction) |
548 | return 999; | |
87da0205 | 549 | |
550 | //TODO Check whether it makes sense to set correctMultiplicity to kTRUE while correctEta might be kFALSE | |
5a9dc560 | 551 | |
552 | // If eta correction (=> new sigma parametrisation) is requested, but no sigma map is available, print error message | |
553 | if (correctEta && !fhEtaSigmaPar1) { | |
554 | AliError("New sigma parametrisation requested, but sigma map not initialised (usually via AliPIDResponse). Old sigma parametrisation will be used!"); | |
555 | } | |
556 | ||
f84b18dd | 557 | // If no sigma map is available or if no eta correction is requested (sigma maps only for corrected eta!), use the old parametrisation |
558 | if (!fhEtaSigmaPar1 || !correctEta) { | |
559 | if (nPoints != 0) | |
87da0205 | 560 | return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity) * |
f84b18dd | 561 | fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints); |
562 | else | |
87da0205 | 563 | return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, correctMultiplicity)*fRes0[gainScenario]; |
f84b18dd | 564 | } |
565 | ||
566 | if (nPoints > 0) { | |
87da0205 | 567 | // Use eta correction (+ eta-dependent sigma) |
ef7661fd | 568 | Double_t sigmaPar1 = GetSigmaPar1Fast(track, species, dEdx, responseFunction); |
f84b18dd | 569 | |
87da0205 | 570 | if (correctMultiplicity) { |
571 | // In addition, take into account multiplicity dependence of mean and sigma of dEdx | |
572 | Double_t dEdxExpectedEtaCorrected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE); | |
573 | ||
574 | // GetMultiplicityCorrection and GetMultiplicitySigmaCorrection both need the eta corrected dEdxExpected | |
ef7661fd | 575 | Double_t multiplicityCorrFactor = GetMultiplicityCorrectionFast(track, dEdxExpectedEtaCorrected, fCurrentEventMultiplicity); |
576 | Double_t multiplicitySigmaCorrFactor = GetMultiplicitySigmaCorrectionFast(dEdxExpectedEtaCorrected, fCurrentEventMultiplicity); | |
87da0205 | 577 | |
578 | // multiplicityCorrFactor to correct dEdxExpected for multiplicity. In addition: Correction factor for sigma | |
579 | return (dEdxExpectedEtaCorrected * multiplicityCorrFactor) | |
580 | * (sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints) * multiplicitySigmaCorrFactor); | |
581 | } | |
582 | else { | |
583 | return GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE)* | |
584 | sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints); | |
585 | } | |
f84b18dd | 586 | } |
587 | else { | |
588 | // One should never have/take tracks with 0 dEdx clusters! | |
589 | return 999; | |
590 | } | |
644666df | 591 | } |
592 | ||
f84b18dd | 593 | |
644666df | 594 | //_________________________________________________________________________ |
595 | Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, | |
596 | AliPID::EParticleType species, | |
f84b18dd | 597 | ETPCdEdxSource dedxSource, |
87da0205 | 598 | Bool_t correctEta, |
599 | Bool_t correctMultiplicity) const | |
644666df | 600 | { |
601 | // Calculates the expected sigma of the PID signal as the function of | |
602 | // the information stored in the track, for the specified particle type | |
603 | // and dedx scenario | |
604 | // | |
605 | ||
f84b18dd | 606 | Double_t dEdx = -1; |
607 | Int_t nPoints = -1; | |
608 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
609 | TSpline3* responseFunction = 0x0; | |
610 | ||
611 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
612 | return 999; //TODO: better handling! | |
613 | ||
87da0205 | 614 | return GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity); |
644666df | 615 | } |
616 | ||
f84b18dd | 617 | |
644666df | 618 | //_________________________________________________________________________ |
619 | Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track, | |
620 | AliPID::EParticleType species, | |
f84b18dd | 621 | ETPCdEdxSource dedxSource, |
87da0205 | 622 | Bool_t correctEta, |
623 | Bool_t correctMultiplicity) const | |
644666df | 624 | { |
f84b18dd | 625 | //Calculates the number of sigmas of the PID signal from the expected value |
626 | //for a given particle species in the presence of multiple gain scenarios | |
644666df | 627 | //inside the TPC |
567624b5 | 628 | |
f84b18dd | 629 | Double_t dEdx = -1; |
630 | Int_t nPoints = -1; | |
631 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
632 | TSpline3* responseFunction = 0x0; | |
633 | ||
634 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
635 | return -999; //TODO: Better handling! | |
567624b5 | 636 | |
87da0205 | 637 | Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity); |
638 | Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta, correctMultiplicity); | |
f84b18dd | 639 | // 999 will be returned by GetExpectedSigma e.g. in case of 0 dEdx clusters |
640 | if (sigma >= 998) | |
641 | return -999; | |
642 | else | |
643 | return (dEdx-bethe)/sigma; | |
644666df | 644 | } |
645 | ||
567624b5 | 646 | //_________________________________________________________________________ |
647 | Float_t AliTPCPIDResponse::GetSignalDelta(const AliVTrack* track, | |
648 | AliPID::EParticleType species, | |
649 | ETPCdEdxSource dedxSource, | |
1d59271b | 650 | Bool_t correctEta, |
87da0205 | 651 | Bool_t correctMultiplicity, |
652 | Bool_t ratio/*=kFALSE*/)const | |
567624b5 | 653 | { |
654 | //Calculates the number of sigmas of the PID signal from the expected value | |
655 | //for a given particle species in the presence of multiple gain scenarios | |
656 | //inside the TPC | |
657 | ||
658 | Double_t dEdx = -1; | |
659 | Int_t nPoints = -1; | |
660 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
661 | TSpline3* responseFunction = 0x0; | |
662 | ||
663 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
664 | return -9999.; //TODO: Better handling! | |
665 | ||
87da0205 | 666 | const Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta, correctMultiplicity); |
1d59271b | 667 | |
668 | Double_t delta=-9999.; | |
669 | if (!ratio) delta=dEdx-bethe; | |
670 | else if (bethe>1.e-20) delta=dEdx/bethe; | |
87da0205 | 671 | |
1d59271b | 672 | return delta; |
567624b5 | 673 | } |
674 | ||
644666df | 675 | //_________________________________________________________________________ |
676 | Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track, | |
677 | AliPID::EParticleType species, | |
f84b18dd | 678 | ETPCdEdxSource dedxSource, |
1b45e564 | 679 | Double_t& dEdx, |
680 | Int_t& nPoints, | |
681 | ETPCgainScenario& gainScenario, | |
f85a3764 | 682 | TSpline3** responseFunction) const |
644666df | 683 | { |
684 | // Calculates the right parameters for PID | |
685 | // dEdx parametrization for the proper gain scenario, dEdx | |
686 | // and NPoints used for dEdx | |
687 | // based on the track geometry (which chambers it crosses) for the specified particle type | |
688 | // and preferred source of dedx. | |
689 | // returns true on success | |
690 | ||
f84b18dd | 691 | |
692 | if (dedxSource == kdEdxDefault) { | |
693 | // Fast handling for default case. In addition: Keep it simple (don't call additional functions) to | |
694 | // avoid possible bugs | |
f85a3764 | 695 | |
696 | // GetTPCsignalTunedOnData will be non-positive, if it has not been set (i.e. in case of MC NOT tuned to data). | |
697 | // If this is the case, just take the normal signal | |
698 | dEdx = track->GetTPCsignalTunedOnData(); | |
699 | if (dEdx <= 0) { | |
700 | dEdx = track->GetTPCsignal(); | |
701 | } | |
702 | ||
f84b18dd | 703 | nPoints = track->GetTPCsignalN(); |
704 | gainScenario = kDefault; | |
705 | ||
706 | TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario)); | |
707 | *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast? | |
708 | ||
709 | return kTRUE; | |
710 | } | |
711 | ||
f85a3764 | 712 | //TODO Proper handle of tuneMConData for other dEdx sources |
f84b18dd | 713 | |
644666df | 714 | Double32_t signal[4]; //0: IROC, 1: OROC medium, 2:OROC long, 3: OROC all (def. truncation used) |
715 | Char_t ncl[3]; //same | |
716 | Char_t nrows[3]; //same | |
f84b18dd | 717 | const AliTPCdEdxInfo* dEdxInfo = track->GetTPCdEdxInfo(); |
644666df | 718 | |
719 | if (!dEdxInfo && dedxSource!=kdEdxDefault) //in one case its ok if we dont have the info | |
720 | { | |
721 | AliError("AliTPCdEdxInfo not available"); | |
644666df | 722 | return kFALSE; |
723 | } | |
724 | ||
725 | if (dEdxInfo) dEdxInfo->GetTPCSignalRegionInfo(signal,ncl,nrows); | |
726 | ||
727 | //check if we cross a bad OROC in which case we reject | |
f85a3764 | 728 | EChamberStatus trackOROCStatus = TrackStatus(track,2); |
729 | if (trackOROCStatus==kChamberOff || trackOROCStatus==kChamberLowGain) | |
644666df | 730 | { |
644666df | 731 | return kFALSE; |
732 | } | |
733 | ||
734 | switch (dedxSource) | |
735 | { | |
644666df | 736 | case kdEdxOROC: |
737 | { | |
f85a3764 | 738 | if (trackOROCStatus==kChamberInvalid) return kFALSE; //never reached OROC |
f84b18dd | 739 | dEdx = signal[3]; |
740 | nPoints = ncl[2]+ncl[1]; | |
741 | gainScenario = kOROChigh; | |
644666df | 742 | break; |
743 | } | |
744 | case kdEdxHybrid: | |
745 | { | |
746 | //if we cross a bad IROC we use OROC dedx, if we dont we use combined | |
747 | EChamberStatus status = TrackStatus(track,1); | |
748 | if (status!=kChamberHighGain) | |
749 | { | |
f84b18dd | 750 | dEdx = signal[3]; |
751 | nPoints = ncl[2]+ncl[1]; | |
752 | gainScenario = kOROChigh; | |
644666df | 753 | } |
754 | else | |
755 | { | |
f84b18dd | 756 | dEdx = track->GetTPCsignal(); |
757 | nPoints = track->GetTPCsignalN(); | |
758 | gainScenario = kALLhigh; | |
644666df | 759 | } |
760 | break; | |
761 | } | |
762 | default: | |
763 | { | |
f84b18dd | 764 | dEdx = 0.; |
765 | nPoints = 0; | |
766 | gainScenario = kGainScenarioInvalid; | |
644666df | 767 | return kFALSE; |
768 | } | |
769 | } | |
f84b18dd | 770 | TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario)); |
771 | *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast? | |
772 | ||
644666df | 773 | return kTRUE; |
774 | } | |
775 | ||
f84b18dd | 776 | |
777 | //_________________________________________________________________________ | |
ef7661fd | 778 | Double_t AliTPCPIDResponse::GetEtaCorrectionFast(const AliVTrack *track, Double_t dEdxSplines) const |
f84b18dd | 779 | { |
ef7661fd | 780 | // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly. |
f84b18dd | 781 | // |
782 | // Get eta correction for the given parameters. | |
783 | // | |
784 | ||
5a9dc560 | 785 | if (!fhEtaCorr) { |
786 | // Calling this function means to request eta correction in some way. Print error message, if no map is available! | |
787 | AliError("Eta correction requested, but map not initialised (usually via AliPIDResponse). Returning eta correction factor 1!"); | |
f84b18dd | 788 | return 1.; |
5a9dc560 | 789 | } |
f84b18dd | 790 | |
791 | Double_t tpcSignal = dEdxSplines; | |
792 | ||
793 | if (tpcSignal < 1.) | |
794 | return 1.; | |
795 | ||
87da0205 | 796 | Double_t tanTheta = GetTrackTanTheta(track); |
7294ba92 | 797 | Int_t binX = fhEtaCorr->GetXaxis()->FindFixBin(tanTheta); |
798 | Int_t binY = fhEtaCorr->GetYaxis()->FindFixBin(1. / tpcSignal); | |
f84b18dd | 799 | |
800 | if (binX == 0) | |
801 | binX = 1; | |
802 | if (binX > fhEtaCorr->GetXaxis()->GetNbins()) | |
803 | binX = fhEtaCorr->GetXaxis()->GetNbins(); | |
804 | ||
805 | if (binY == 0) | |
806 | binY = 1; | |
807 | if (binY > fhEtaCorr->GetYaxis()->GetNbins()) | |
808 | binY = fhEtaCorr->GetYaxis()->GetNbins(); | |
809 | ||
810 | return fhEtaCorr->GetBinContent(binX, binY); | |
811 | } | |
812 | ||
813 | ||
814 | //_________________________________________________________________________ | |
815 | Double_t AliTPCPIDResponse::GetEtaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const | |
816 | { | |
817 | // | |
818 | // Get eta correction for the given track. | |
819 | // | |
820 | ||
821 | if (!fhEtaCorr) | |
822 | return 1.; | |
823 | ||
824 | Double_t dEdx = -1; | |
825 | Int_t nPoints = -1; | |
826 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
827 | TSpline3* responseFunction = 0x0; | |
828 | ||
829 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
830 | return 1.; | |
831 | ||
87da0205 | 832 | // For the eta correction, do NOT take the multiplicity corrected value of dEdx |
833 | Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE); | |
f84b18dd | 834 | |
835 | //TODO Alternatively take current track dEdx | |
ef7661fd | 836 | //return GetEtaCorrectionFast(track, dEdx); |
f84b18dd | 837 | |
ef7661fd | 838 | return GetEtaCorrectionFast(track, dEdxSplines); |
f84b18dd | 839 | } |
840 | ||
841 | ||
842 | //_________________________________________________________________________ | |
843 | Double_t AliTPCPIDResponse::GetEtaCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const | |
844 | { | |
845 | // | |
846 | // Get eta corrected dEdx for the given track. For the correction, the expected dEdx of | |
847 | // the specified species will be used. If the species is set to AliPID::kUnknown, the | |
848 | // dEdx of the track is used instead. | |
849 | // WARNING: In the latter case, the eta correction might not be as good as if the | |
850 | // expected dEdx is used, which is the way the correction factor is designed | |
851 | // for. | |
852 | // In any case, one has to decide carefully to which expected signal one wants to | |
853 | // compare the corrected value - to the corrected or uncorrected. | |
854 | // Anyhow, a safer way of looking e.g. at the n-sigma is to call | |
855 | // the corresponding function GetNumberOfSigmas! | |
856 | // | |
857 | ||
858 | Double_t dEdx = -1; | |
859 | Int_t nPoints = -1; | |
860 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
861 | TSpline3* responseFunction = 0x0; | |
862 | ||
863 | // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case | |
864 | // it is not used anyway, so this causes no trouble. | |
865 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
866 | return -1.; | |
867 | ||
868 | Double_t etaCorr = 0.; | |
869 | ||
870 | if (species < AliPID::kUnknown) { | |
87da0205 | 871 | // For the eta correction, do NOT take the multiplicity corrected value of dEdx |
872 | Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE); | |
ef7661fd | 873 | etaCorr = GetEtaCorrectionFast(track, dEdxSplines); |
f84b18dd | 874 | } |
875 | else { | |
ef7661fd | 876 | etaCorr = GetEtaCorrectionFast(track, dEdx); |
f84b18dd | 877 | } |
878 | ||
879 | if (etaCorr <= 0) | |
880 | return -1.; | |
881 | ||
882 | return dEdx / etaCorr; | |
883 | } | |
884 | ||
885 | ||
f84b18dd | 886 | //_________________________________________________________________________ |
ef7661fd | 887 | Double_t AliTPCPIDResponse::GetSigmaPar1Fast(const AliVTrack *track, AliPID::EParticleType species, Double_t dEdx, |
888 | const TSpline3* responseFunction) const | |
f84b18dd | 889 | { |
ef7661fd | 890 | // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly. |
f84b18dd | 891 | // |
892 | // Get parameter 1 of sigma parametrisation of TPC dEdx from the histogram for the given track. | |
893 | // | |
894 | ||
5a9dc560 | 895 | if (!fhEtaSigmaPar1) { |
896 | // Calling this function means to request new sigma parametrisation in some way. Print error message, if no map is available! | |
897 | AliError("New sigma parametrisation requested, but sigma map not initialised (usually via AliPIDResponse). Returning error value for sigma parameter1 = 999!"); | |
f84b18dd | 898 | return 999; |
5a9dc560 | 899 | } |
f84b18dd | 900 | |
901 | // The sigma maps are created with data sets that are already eta corrected and for which the | |
902 | // splines have been re-created. Therefore, the value for the lookup needs to be the value of | |
903 | // the splines without any additional eta correction. | |
904 | // NOTE: This is due to the method the maps are created. The track dEdx (not the expected one!) | |
905 | // is corrected to uniquely related a momemtum bin with an expected dEdx, where the expected dEdx | |
906 | // equals the track dEdx for all eta (thanks to the correction and the re-fit of the splines). | |
907 | // Consequently, looking up the uncorrected expected dEdx at a given tanTheta yields the correct | |
908 | // sigma parameter! | |
909 | // Also: It has to be the spline dEdx, since one wants to get the sigma for the assumption(!) | |
910 | // of such a particle, which by assumption then has this dEdx value | |
87da0205 | 911 | |
912 | // For the eta correction, do NOT take the multiplicity corrected value of dEdx | |
913 | Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE); | |
f84b18dd | 914 | |
915 | if (dEdxExpected < 1.) | |
916 | return 999; | |
917 | ||
ef7661fd | 918 | Double_t tanTheta = GetTrackTanTheta(track); |
7294ba92 | 919 | Int_t binX = fhEtaSigmaPar1->GetXaxis()->FindFixBin(tanTheta); |
920 | Int_t binY = fhEtaSigmaPar1->GetYaxis()->FindFixBin(1. / dEdxExpected); | |
f84b18dd | 921 | |
922 | if (binX == 0) | |
923 | binX = 1; | |
924 | if (binX > fhEtaSigmaPar1->GetXaxis()->GetNbins()) | |
925 | binX = fhEtaSigmaPar1->GetXaxis()->GetNbins(); | |
926 | ||
927 | if (binY == 0) | |
928 | binY = 1; | |
929 | if (binY > fhEtaSigmaPar1->GetYaxis()->GetNbins()) | |
930 | binY = fhEtaSigmaPar1->GetYaxis()->GetNbins(); | |
931 | ||
932 | return fhEtaSigmaPar1->GetBinContent(binX, binY); | |
933 | } | |
934 | ||
935 | ||
936 | //_________________________________________________________________________ | |
937 | Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const | |
938 | { | |
939 | // | |
87da0205 | 940 | // Get parameter 1 of sigma parametrisation of TPC dEdx from the histogram for the given track. |
f84b18dd | 941 | // |
942 | ||
943 | if (!fhEtaSigmaPar1) | |
944 | return 999; | |
945 | ||
946 | Double_t dEdx = -1; | |
947 | Int_t nPoints = -1; | |
948 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
949 | TSpline3* responseFunction = 0x0; | |
950 | ||
951 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
952 | return 999; | |
953 | ||
ef7661fd | 954 | return GetSigmaPar1Fast(track, species, dEdx, responseFunction); |
f84b18dd | 955 | } |
956 | ||
957 | ||
958 | //_________________________________________________________________________ | |
959 | Bool_t AliTPCPIDResponse::SetEtaCorrMap(TH2D* hMap) | |
960 | { | |
961 | // | |
962 | // Load map for TPC eta correction (a copy is stored and will be deleted automatically). | |
963 | // If hMap is 0x0,the eta correction will be disabled and kFALSE is returned. | |
964 | // If the map can be set, kTRUE is returned. | |
965 | // | |
966 | ||
967 | delete fhEtaCorr; | |
968 | ||
969 | if (!hMap) { | |
970 | fhEtaCorr = 0x0; | |
971 | ||
972 | return kFALSE; | |
973 | } | |
974 | ||
975 | fhEtaCorr = (TH2D*)(hMap->Clone()); | |
976 | fhEtaCorr->SetDirectory(0); | |
977 | ||
978 | return kTRUE; | |
979 | } | |
980 | ||
87da0205 | 981 | |
f84b18dd | 982 | //_________________________________________________________________________ |
983 | Bool_t AliTPCPIDResponse::SetSigmaParams(TH2D* hSigmaPar1Map, Double_t sigmaPar0) | |
984 | { | |
985 | // | |
986 | // Load map for TPC sigma map (a copy is stored and will be deleted automatically): | |
987 | // Parameter 1 is stored as a 2D map (1/dEdx vs. tanTheta_local) and parameter 0 is | |
988 | // a constant. If hSigmaPar1Map is 0x0, the old sigma parametrisation will be used | |
989 | // (and sigmaPar0 is ignored!) and kFALSE is returned. | |
990 | // If the map can be set, sigmaPar0 is also set and kTRUE will be returned. | |
991 | // | |
992 | ||
993 | delete fhEtaSigmaPar1; | |
994 | ||
995 | if (!hSigmaPar1Map) { | |
996 | fhEtaSigmaPar1 = 0x0; | |
997 | fSigmaPar0 = 0.0; | |
998 | ||
999 | return kFALSE; | |
1000 | } | |
1001 | ||
1002 | fhEtaSigmaPar1 = (TH2D*)(hSigmaPar1Map->Clone()); | |
1003 | fhEtaSigmaPar1->SetDirectory(0); | |
1004 | fSigmaPar0 = sigmaPar0; | |
1005 | ||
1006 | return kTRUE; | |
1007 | } | |
1008 | ||
1009 | ||
87da0205 | 1010 | //_________________________________________________________________________ |
1011 | Double_t AliTPCPIDResponse::GetTrackTanTheta(const AliVTrack *track) const | |
1012 | { | |
1013 | // Extract the tanTheta from the information available in the AliVTrack | |
1014 | ||
1015 | // For ESD tracks, the local tanTheta could be used (esdTrack->GetInnerParam()->GetTgl()). | |
1016 | // However, this value is not available for AODs and, thus, not for AliVTrack. | |
1017 | // Fortunately, the following formula allows to approximate the local tanTheta with the | |
1018 | // global theta angle -> This is for by far most of the tracks the same, but gives at | |
1019 | // maybe the percent level differences within +- 0.2 in tanTheta -> Which is still ok. | |
1020 | ||
ef7661fd | 1021 | /* |
1022 | const AliExternalTrackParam* innerParam = track->GetInnerParam(); | |
1023 | Double_t tanTheta = 0; | |
1024 | if (innerParam) | |
1025 | tanTheta = innerParam->GetTgl(); | |
1026 | else | |
1027 | tanTheta = TMath::Tan(-track->Theta() + TMath::Pi() / 2.0); | |
1028 | ||
1029 | // Constant in formula for B in kGauss (factor 0.1 to convert B from Tesla to kGauss), | |
1030 | // pT in GeV/c (factor c*1e-9 to arrive at GeV/c) and curvature in 1/cm (factor 0.01 to get from m to cm) | |
1031 | const Double_t constant = TMath::C()* 1e-9 * 0.1 * 0.01; | |
1032 | const Double_t curvature = fMagField * constant / track->Pt(); // in 1./cm | |
1033 | ||
1034 | Double_t averageddzdr = 0.; | |
1035 | Int_t nParts = 0; | |
1036 | ||
1037 | for (Double_t r = 85; r < 245; r++) { | |
1038 | Double_t sinPhiLocal = TMath::Abs(r*curvature*0.5); | |
1039 | ||
1040 | // Cut on |sin(phi)| as during reco | |
1041 | if (TMath::Abs(sinPhiLocal) <= 0.95) { | |
1042 | const Double_t phiLocal = TMath::ASin(sinPhiLocal); | |
1043 | const Double_t tanPhiLocal = TMath::Tan(phiLocal); | |
1044 | ||
1045 | averageddzdr += tanTheta * TMath::Sqrt(1. + tanPhiLocal * tanPhiLocal); | |
1046 | nParts++; | |
1047 | } | |
1048 | } | |
1049 | ||
1050 | if (nParts > 0) | |
1051 | averageddzdr /= nParts; | |
1052 | else { | |
1053 | AliError("Problems during determination of dz/dr. Returning pure tanTheta as best estimate!"); | |
1054 | return tanTheta; | |
1055 | } | |
1056 | ||
1057 | //printf("pT: %f\nFactor/magField(kGs)/curvature^-1: %f / %f /%f\ntanThetaGlobalFromTheta/tanTheta/Averageddzdr: %f / %f / %f\n\n", | |
1058 | // track->Pt(), constant, fMagField, 1./curvature, TMath::Tan(-track->Theta() + TMath::Pi() / 2.0), tanTheta, averageddzdr); | |
1059 | ||
1060 | return averageddzdr; | |
1061 | */ | |
1062 | ||
87da0205 | 1063 | |
1064 | // Alternatively (in average, the effect is found to be negligable!): | |
1065 | // Take local tanTheta from TPC inner wall, if available (currently only for ESDs available) | |
ef7661fd | 1066 | //const AliExternalTrackParam* innerParam = track->GetInnerParam(); |
1067 | //if (innerParam) { | |
1068 | // return innerParam->GetTgl(); | |
1069 | //} | |
87da0205 | 1070 | |
1071 | return TMath::Tan(-track->Theta() + TMath::Pi() / 2.0); | |
1072 | } | |
1073 | ||
1074 | ||
1075 | //_________________________________________________________________________ | |
205d83d8 | 1076 | Double_t AliTPCPIDResponse::GetMultiplicityCorrectionFast(const AliVTrack *track, Double_t dEdxExpected, Int_t multiplicity) const |
87da0205 | 1077 | { |
ef7661fd | 1078 | // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly. |
1079 | // | |
87da0205 | 1080 | // Calculate the multiplicity correction factor for this track for the given multiplicity. |
1081 | // The parameter dEdxExpected should take into account the eta correction already! | |
1082 | ||
1083 | // Multiplicity depends on pure dEdx. Therefore, correction factor depends indirectly on eta | |
1084 | // => Use eta corrected value for dEdexExpected. | |
1085 | ||
1086 | if (dEdxExpected <= 0 || multiplicity <= 0) | |
1087 | return 1.0; | |
1088 | ||
1089 | const Double_t dEdxExpectedInv = 1. / dEdxExpected; | |
1090 | Double_t relSlope = fCorrFuncMultiplicity->Eval(dEdxExpectedInv); | |
1091 | ||
1092 | const Double_t tanTheta = GetTrackTanTheta(track); | |
1093 | relSlope += fCorrFuncMultiplicityTanTheta->Eval(tanTheta); | |
1094 | ||
1095 | return (1. + relSlope * multiplicity); | |
1096 | } | |
1097 | ||
1098 | ||
1099 | //_________________________________________________________________________ | |
1100 | Double_t AliTPCPIDResponse::GetMultiplicityCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const | |
1101 | { | |
1102 | // | |
1103 | // Get multiplicity correction for the given track (w.r.t. the mulitplicity of the current event) | |
1104 | // | |
1105 | ||
1106 | //TODO Should return error value, if no eta correction is enabled (needs eta correction). OR: Reset correction in case of no eta correction in PIDresponse | |
1107 | ||
1108 | // No correction in case of multiplicity <= 0 | |
1109 | if (fCurrentEventMultiplicity <= 0) | |
1110 | return 1.; | |
1111 | ||
1112 | Double_t dEdx = -1; | |
1113 | Int_t nPoints = -1; | |
1114 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
1115 | TSpline3* responseFunction = 0x0; | |
1116 | ||
1117 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
1118 | return 1.; | |
1119 | ||
1120 | //TODO Does it make sense to use the multiplicity correction WITHOUT eta correction?! Are the fit parameters still valid? | |
1121 | // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course) | |
1122 | Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE); | |
1123 | ||
ef7661fd | 1124 | return GetMultiplicityCorrectionFast(track, dEdxExpected, fCurrentEventMultiplicity); |
87da0205 | 1125 | } |
1126 | ||
1127 | ||
1128 | //_________________________________________________________________________ | |
1129 | Double_t AliTPCPIDResponse::GetMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const | |
1130 | { | |
1131 | // | |
1132 | // Get multiplicity corrected dEdx for the given track. For the correction, the expected dEdx of | |
1133 | // the specified species will be used. If the species is set to AliPID::kUnknown, the | |
1134 | // dEdx of the track is used instead. | |
1135 | // WARNING: In the latter case, the correction might not be as good as if the | |
1136 | // expected dEdx is used, which is the way the correction factor is designed | |
1137 | // for. | |
1138 | // In any case, one has to decide carefully to which expected signal one wants to | |
1139 | // compare the corrected value - to the corrected or uncorrected. | |
1140 | // Anyhow, a safer way of looking e.g. at the n-sigma is to call | |
1141 | // the corresponding function GetNumberOfSigmas! | |
1142 | // | |
1143 | ||
1144 | //TODO Should return error value, if no eta correction is enabled (needs eta correction). OR: Reset correction in case of no eta correction in PIDresponse | |
1145 | ||
1146 | Double_t dEdx = -1; | |
1147 | Int_t nPoints = -1; | |
1148 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
1149 | TSpline3* responseFunction = 0x0; | |
1150 | ||
1151 | // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case | |
1152 | // it is not used anyway, so this causes no trouble. | |
1153 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
1154 | return -1.; | |
1155 | ||
1156 | ||
1157 | // No correction in case of multiplicity <= 0 | |
1158 | if (fCurrentEventMultiplicity <= 0) | |
1159 | return dEdx; | |
1160 | ||
1161 | Double_t multiplicityCorr = 0.; | |
1162 | ||
1163 | // TODO Normally, one should use the eta corrected values in BOTH of the following cases. Does it make sense to use the uncorrected dEdx values? | |
1164 | if (species < AliPID::kUnknown) { | |
1165 | // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course). | |
1166 | // However, one needs the eta corrected value! | |
1167 | Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE); | |
ef7661fd | 1168 | multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdxSplines, fCurrentEventMultiplicity); |
87da0205 | 1169 | } |
1170 | else { | |
1171 | // One needs the eta corrected value to determine the multiplicity correction factor! | |
ef7661fd | 1172 | Double_t etaCorr = GetEtaCorrectionFast(track, dEdx); |
1173 | multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdx * etaCorr, fCurrentEventMultiplicity); | |
87da0205 | 1174 | } |
1175 | ||
1176 | if (multiplicityCorr <= 0) | |
1177 | return -1.; | |
1178 | ||
1179 | return dEdx / multiplicityCorr; | |
1180 | } | |
1181 | ||
1182 | ||
1183 | //_________________________________________________________________________ | |
1184 | Double_t AliTPCPIDResponse::GetEtaAndMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, | |
1185 | ETPCdEdxSource dedxSource) const | |
1186 | { | |
1187 | // | |
1188 | // Get multiplicity and eta corrected dEdx for the given track. For the correction, | |
1189 | // the expected dEdx of the specified species will be used. If the species is set | |
1190 | // to AliPID::kUnknown, the dEdx of the track is used instead. | |
1191 | // WARNING: In the latter case, the correction might not be as good as if the | |
1192 | // expected dEdx is used, which is the way the correction factor is designed | |
1193 | // for. | |
1194 | // In any case, one has to decide carefully to which expected signal one wants to | |
1195 | // compare the corrected value - to the corrected or uncorrected. | |
1196 | // Anyhow, a safer way of looking e.g. at the n-sigma is to call | |
1197 | // the corresponding function GetNumberOfSigmas! | |
1198 | // | |
1199 | ||
1200 | Double_t dEdx = -1; | |
1201 | Int_t nPoints = -1; | |
1202 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
1203 | TSpline3* responseFunction = 0x0; | |
1204 | ||
1205 | // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case | |
1206 | // it is not used anyway, so this causes no trouble. | |
1207 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
1208 | return -1.; | |
1209 | ||
1210 | Double_t multiplicityCorr = 0.; | |
1211 | Double_t etaCorr = 0.; | |
1212 | ||
1213 | if (species < AliPID::kUnknown) { | |
1214 | // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course) | |
1215 | Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE, kFALSE); | |
ef7661fd | 1216 | etaCorr = GetEtaCorrectionFast(track, dEdxSplines); |
1217 | multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdxSplines * etaCorr, fCurrentEventMultiplicity); | |
87da0205 | 1218 | } |
1219 | else { | |
ef7661fd | 1220 | etaCorr = GetEtaCorrectionFast(track, dEdx); |
1221 | multiplicityCorr = GetMultiplicityCorrectionFast(track, dEdx * etaCorr, fCurrentEventMultiplicity); | |
87da0205 | 1222 | } |
1223 | ||
1224 | if (multiplicityCorr <= 0 || etaCorr <= 0) | |
1225 | return -1.; | |
1226 | ||
1227 | return dEdx / multiplicityCorr / etaCorr; | |
1228 | } | |
1229 | ||
1230 | ||
1231 | //_________________________________________________________________________ | |
205d83d8 | 1232 | Double_t AliTPCPIDResponse::GetMultiplicitySigmaCorrectionFast(Double_t dEdxExpected, Int_t multiplicity) const |
87da0205 | 1233 | { |
ef7661fd | 1234 | // NOTE: For expert use only -> Non-experts are advised to use the function without the "Fast" suffix or stick to AliPIDResponse directly. |
1235 | // | |
87da0205 | 1236 | // Calculate the multiplicity sigma correction factor for the corresponding expected dEdx and for the given multiplicity. |
1237 | // The parameter dEdxExpected should take into account the eta correction already! | |
1238 | ||
1239 | // Multiplicity dependence of sigma depends on the real dEdx at zero multiplicity, | |
1240 | // i.e. the eta (only) corrected dEdxExpected value has to be used | |
1241 | // since all maps etc. have been created for ~zero multiplicity | |
1242 | ||
1243 | if (dEdxExpected <= 0 || multiplicity <= 0) | |
1244 | return 1.0; | |
1245 | ||
1246 | Double_t relSigmaSlope = fCorrFuncSigmaMultiplicity->Eval(1. / dEdxExpected); | |
1247 | ||
1248 | return (1. + relSigmaSlope * multiplicity); | |
1249 | } | |
1250 | ||
1251 | ||
1252 | //_________________________________________________________________________ | |
1253 | Double_t AliTPCPIDResponse::GetMultiplicitySigmaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const | |
1254 | { | |
1255 | // | |
1256 | // Get multiplicity sigma correction for the given track (w.r.t. the mulitplicity of the current event) | |
1257 | // | |
1258 | ||
1259 | //TODO Should return error value, if no eta correction is enabled (needs eta correction). OR: Reset correction in case of no eta correction in PIDresponse | |
1260 | ||
1261 | // No correction in case of multiplicity <= 0 | |
1262 | if (fCurrentEventMultiplicity <= 0) | |
1263 | return 1.; | |
1264 | ||
1265 | Double_t dEdx = -1; | |
1266 | Int_t nPoints = -1; | |
1267 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
1268 | TSpline3* responseFunction = 0x0; | |
1269 | ||
1270 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
1271 | return 1.; | |
1272 | ||
1273 | //TODO Does it make sense to use the multiplicity correction WITHOUT eta correction?! Are the fit parameters still valid? | |
1274 | // To get the expected signal to determine the multiplicity correction, do NOT ask for the multiplicity corrected value (of course) | |
1275 | Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE, kFALSE); | |
1276 | ||
ef7661fd | 1277 | return GetMultiplicitySigmaCorrectionFast(dEdxExpected, fCurrentEventMultiplicity); |
87da0205 | 1278 | } |
1279 | ||
1280 | ||
1281 | //_________________________________________________________________________ | |
1282 | void AliTPCPIDResponse::ResetMultiplicityCorrectionFunctions() | |
1283 | { | |
1284 | // Default values: No correction, i.e. overall correction factor should be one | |
1285 | ||
1286 | // This function should always return 0.0, if multcorr disabled | |
1287 | fCorrFuncMultiplicity->SetParameter(0, 0.); | |
1288 | fCorrFuncMultiplicity->SetParameter(1, 0.); | |
1289 | fCorrFuncMultiplicity->SetParameter(2, 0.); | |
1290 | fCorrFuncMultiplicity->SetParameter(3, 0.); | |
1291 | fCorrFuncMultiplicity->SetParameter(4, 0.); | |
1292 | ||
1293 | // This function should always return 0., if multCorr disabled | |
1294 | fCorrFuncMultiplicityTanTheta->SetParameter(0, 0.); | |
1295 | fCorrFuncMultiplicityTanTheta->SetParameter(1, 0.); | |
1296 | fCorrFuncMultiplicityTanTheta->SetParameter(2, 0.); | |
1297 | ||
1298 | // This function should always return 0.0, if mutlcorr disabled | |
1299 | fCorrFuncSigmaMultiplicity->SetParameter(0, 0.); | |
1300 | fCorrFuncSigmaMultiplicity->SetParameter(1, 0.); | |
1301 | fCorrFuncSigmaMultiplicity->SetParameter(2, 0.); | |
1302 | fCorrFuncSigmaMultiplicity->SetParameter(3, 0.); | |
1303 | } | |
1304 | ||
1305 | ||
644666df | 1306 | //_________________________________________________________________________ |
f85a3764 | 1307 | Bool_t AliTPCPIDResponse::sectorNumbersInOut(Double_t* trackPositionInner, |
1308 | Double_t* trackPositionOuter, | |
1309 | Float_t& inphi, | |
644666df | 1310 | Float_t& outphi, |
f85a3764 | 1311 | Int_t& in, |
1b45e564 | 1312 | Int_t& out ) const |
644666df | 1313 | { |
1314 | //calculate the sector numbers (equivalent to IROC chamber numbers) a track crosses | |
1315 | //for OROC chamber numbers add 36 | |
1316 | //returned angles are between (0,2pi) | |
644666df | 1317 | |
1318 | inphi = TMath::ATan2(trackPositionInner[1],trackPositionInner[0]); | |
1319 | outphi = TMath::ATan2(trackPositionOuter[1], trackPositionOuter[0]); | |
1320 | ||
1321 | if (inphi<0) {inphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi | |
1322 | if (outphi<0) {outphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi | |
1323 | ||
1324 | in = sectorNumber(inphi); | |
1325 | out = sectorNumber(outphi); | |
1326 | ||
1327 | //for the C side (positive z) offset by 18 | |
1328 | if (trackPositionInner[2]>0.0) | |
1329 | { | |
1330 | in+=18; | |
1331 | out+=18; | |
1332 | } | |
1333 | return kTRUE; | |
1334 | } | |
1b45e564 | 1335 | |
1336 | ||
644666df | 1337 | //_____________________________________________________________________________ |
1338 | Int_t AliTPCPIDResponse::sectorNumber(Double_t phi) const | |
1339 | { | |
1340 | //calculate sector number | |
1341 | const Float_t width=TMath::TwoPi()/18.; | |
1342 | return TMath::Floor(phi/width); | |
1343 | } | |
1344 | ||
1b45e564 | 1345 | |
644666df | 1346 | //_____________________________________________________________________________ |
1347 | void AliTPCPIDResponse::Print(Option_t* /*option*/) const | |
1348 | { | |
1349 | //Print info | |
1350 | fResponseFunctions.Print(); | |
1351 | } | |
1352 | ||
1b45e564 | 1353 | |
644666df | 1354 | //_____________________________________________________________________________ |
1355 | AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack* track, Int_t layer) const | |
1356 | { | |
1357 | //status of the track: if it crosses any bad chambers return kChamberOff; | |
1358 | //IROC:layer=1, OROC:layer=2 | |
1359 | if (layer<1 || layer>2) layer=1; | |
1360 | Int_t in=0; | |
1361 | Int_t out=0; | |
1362 | Float_t inphi=0.; | |
1363 | Float_t outphi=0.; | |
1364 | Float_t innerRadius = (layer==1)?83.0:133.7; | |
1365 | Float_t outerRadius = (layer==1)?133.5:247.7; | |
f85a3764 | 1366 | |
1367 | ///////////////////////////////////////////////////////////////////////////// | |
1368 | //find out where track enters and leaves the layer. | |
1369 | // | |
1b45e564 | 1370 | Double_t trackPositionInner[3]; |
1371 | Double_t trackPositionOuter[3]; | |
1372 | ||
f85a3764 | 1373 | //if there is no inner param this could mean we're using the AOD track, |
1374 | //we still can extrapolate from the vertex - so use those params. | |
1375 | const AliExternalTrackParam* ip = track->GetInnerParam(); | |
50ff8c36 | 1376 | if (ip) track=dynamic_cast<const AliVTrack*>(ip); |
f85a3764 | 1377 | |
1378 | Bool_t trackAtInner = track->GetXYZAt(innerRadius, fMagField, trackPositionInner); | |
1379 | Bool_t trackAtOuter = track->GetXYZAt(outerRadius, fMagField, trackPositionOuter); | |
1380 | ||
1381 | if (!trackAtInner) | |
1382 | { | |
1383 | //if we dont even enter inner radius we do nothing and return invalid | |
1384 | inphi=0.0; | |
1385 | outphi=0.0; | |
1386 | in=0; | |
1387 | out=0; | |
1388 | return kChamberInvalid; | |
1389 | } | |
1390 | ||
1391 | if (!trackAtOuter) | |
1392 | { | |
1393 | //if we don't reach the outer radius check that the apex is indeed within the outer radius and use apex position | |
1394 | Bool_t haveApex = TrackApex(track, fMagField, trackPositionOuter); | |
1395 | Float_t apexRadius = TMath::Sqrt(trackPositionOuter[0]*trackPositionOuter[0]+trackPositionOuter[1]*trackPositionOuter[1]); | |
1396 | if ( haveApex && apexRadius<=outerRadius && apexRadius>innerRadius) | |
1397 | { | |
1398 | //printf("pt: %.2f, apexRadius: %.2f(%s), x: %.2f, y: %.2f\n",track->Pt(),apexRadius,(haveApex)?"OK":"BAD",trackPositionOuter[0],trackPositionOuter[1]); | |
1399 | } | |
1400 | else | |
1401 | { | |
1402 | inphi=0.0; | |
1403 | outphi=0.0; | |
1404 | in=0; | |
1405 | out=0; | |
1406 | return kChamberInvalid; | |
1407 | } | |
1408 | } | |
1409 | ||
1410 | ||
1b45e564 | 1411 | if (!sectorNumbersInOut(trackPositionInner, |
1412 | trackPositionOuter, | |
1413 | inphi, | |
1414 | outphi, | |
1415 | in, | |
f85a3764 | 1416 | out)) return kChamberInvalid; |
644666df | 1417 | |
f85a3764 | 1418 | ///////////////////////////////////////////////////////////////////////////// |
1b45e564 | 1419 | //now we have the location of the track we can check |
f85a3764 | 1420 | //if it is in a good/bad chamber |
1421 | // | |
644666df | 1422 | Bool_t sideA = kTRUE; |
1b45e564 | 1423 | |
644666df | 1424 | if (((in/18)%2==1) && ((out/18)%2==1)) sideA=kFALSE; |
1425 | ||
1426 | in=in%18; | |
1427 | out=out%18; | |
1428 | ||
1429 | if (TMath::Abs(in-out)>9) | |
1430 | { | |
1431 | if (TMath::Max(in,out)==out) | |
1432 | { | |
1433 | Int_t tmp=out; | |
1434 | out = in; | |
1435 | in = tmp; | |
1436 | Float_t tmpphi=outphi; | |
1437 | outphi=inphi; | |
1438 | inphi=tmpphi; | |
1439 | } | |
1440 | in-=18; | |
1441 | inphi-=TMath::TwoPi(); | |
1442 | } | |
10d100d4 | 1443 | else |
644666df | 1444 | { |
1445 | if (TMath::Max(in,out)==in) | |
1446 | { | |
1447 | Int_t tmp=out; | |
1448 | out=in; | |
1449 | in=tmp; | |
1450 | Float_t tmpphi=outphi; | |
1451 | outphi=inphi; | |
1452 | inphi=tmpphi; | |
1453 | } | |
1454 | } | |
1455 | ||
1456 | Float_t trackLengthInBad = 0.; | |
1457 | Float_t trackLengthInLowGain = 0.; | |
1458 | Float_t trackLengthTotal = TMath::Abs(outphi-inphi); | |
f85a3764 | 1459 | Float_t lengthFractionInBadSectors = 0.; |
644666df | 1460 | |
1461 | const Float_t sectorWidth = TMath::TwoPi()/18.; | |
1b45e564 | 1462 | |
644666df | 1463 | for (Int_t i=in; i<=out; i++) |
1464 | { | |
1465 | int j=i; | |
1466 | if (i<0) j+=18; //correct for the negative values | |
1467 | if (!sideA) j+=18; //move to the correct side | |
1b45e564 | 1468 | |
644666df | 1469 | Float_t deltaPhi = 0.; |
1470 | Float_t phiEdge=sectorWidth*i; | |
1471 | if (inphi>phiEdge) {deltaPhi=phiEdge+sectorWidth-inphi;} | |
1472 | else if ((outphi>=phiEdge) && (outphi<(phiEdge+sectorWidth))) {deltaPhi=outphi-phiEdge;} | |
1473 | else {deltaPhi=sectorWidth;} | |
1b45e564 | 1474 | |
644666df | 1475 | Float_t v = fVoltageMap[(layer==1)?(j):(j+36)]; |
1b45e564 | 1476 | if (v<=fBadIROCthreshhold) |
1477 | { | |
1478 | trackLengthInBad+=deltaPhi; | |
f85a3764 | 1479 | lengthFractionInBadSectors=1.; |
1480 | } | |
1b45e564 | 1481 | if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold) |
f85a3764 | 1482 | { |
1483 | trackLengthInLowGain+=deltaPhi; | |
1484 | lengthFractionInBadSectors=1.; | |
1485 | } | |
644666df | 1486 | } |
1487 | ||
1488 | //for now low gain and bad (off) chambers are treated equally | |
f85a3764 | 1489 | if (trackLengthTotal>0) |
1490 | lengthFractionInBadSectors = (trackLengthInLowGain+trackLengthInBad)/trackLengthTotal; | |
644666df | 1491 | |
1492 | //printf("### side: %s, pt: %.2f, pz: %.2f, in: %i, out: %i, phiIN: %.2f, phiOUT: %.2f, rIN: %.2f, rOUT: %.2f\n",(sideA)?"A":"C",track->Pt(),track->Pz(),in,out,inphi,outphi,innerRadius,outerRadius); | |
1b45e564 | 1493 | |
1494 | if (lengthFractionInBadSectors>fMaxBadLengthFraction) | |
644666df | 1495 | { |
1496 | //printf("%%%%%%%% %s kChamberLowGain\n",(layer==1)?"IROC":"OROC"); | |
1497 | return kChamberLowGain; | |
1498 | } | |
1b45e564 | 1499 | |
644666df | 1500 | return kChamberHighGain; |
1501 | } | |
1502 | ||
1b45e564 | 1503 | |
644666df | 1504 | //_____________________________________________________________________________ |
1505 | Float_t AliTPCPIDResponse::MaxClusterRadius(const AliVTrack* track) const | |
1506 | { | |
1507 | //return the radius of the outermost padrow containing a cluster in TPC | |
1508 | //for the track | |
1509 | const TBits* clusterMap=track->GetTPCClusterMapPtr(); | |
1510 | if (!clusterMap) return 0.; | |
1511 | ||
1512 | //from AliTPCParam, radius of first IROC row | |
1b45e564 | 1513 | const Float_t rfirstIROC = 8.52249984741210938e+01; |
644666df | 1514 | const Float_t padrowHeightIROC = 0.75; |
1515 | const Float_t rfirstOROC0 = 1.35100006103515625e+02; | |
1516 | const Float_t padrowHeightOROC0 = 1.0; | |
1517 | const Float_t rfirstOROC1 = 1.99350006103515625e+02; | |
1518 | const Float_t padrowHeightOROC1 = 1.5; | |
1519 | ||
1520 | Int_t maxPadRow=160; | |
1521 | while ((--maxPadRow)>0 && !clusterMap->TestBitNumber(maxPadRow)){} | |
1522 | if (maxPadRow>126) return rfirstOROC1+(maxPadRow-126-1)*padrowHeightOROC1; | |
1523 | if (maxPadRow>62) return rfirstOROC0+(maxPadRow-62-1)*padrowHeightOROC0; | |
1524 | if (maxPadRow>0) return rfirstIROC+(maxPadRow-1)*padrowHeightIROC; | |
1525 | return 0.0; | |
10d100d4 | 1526 | } |
644666df | 1527 | |
1b45e564 | 1528 | |
644666df | 1529 | //_____________________________________________________________________________ |
1530 | Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const | |
1531 | { | |
1532 | //calculate the coordinates of the apex of the track | |
1533 | Double_t x[3]; | |
1534 | track->GetXYZ(x); | |
1535 | Double_t p[3]; | |
1536 | track->GetPxPyPz(p); | |
1537 | Double_t r = 1./track->OneOverPt()/0.0299792458/magField; //signed - will determine the direction of b | |
1b45e564 | 1538 | //printf("b: %.2f, x:%.2f, y:%.2f, pt: %.2f, px:%.2f, py%.2f, r: %.2f\n",magField, x[0],x[1],track->Pt(), p[0],p[1],r); |
644666df | 1539 | //find orthogonal vector (Gram-Schmidt) |
1540 | Double_t alpha = (p[0]*x[0] + p[1]*x[1])/(p[0]*p[0] + p[1]*p[1]); | |
1541 | Double_t b[2]; | |
1b45e564 | 1542 | b[0]=x[0]-alpha*p[0]; |
1543 | b[1]=x[1]-alpha*p[1]; | |
1544 | ||
644666df | 1545 | Double_t norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]); |
1546 | if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE; | |
1b45e564 | 1547 | b[0]/=norm; |
1548 | b[1]/=norm; | |
1549 | b[0]*=r; | |
1550 | b[1]*=r; | |
1551 | b[0]+=x[0]; | |
644666df | 1552 | b[1]+=x[1]; |
1553 | //printf("center: x:%.2f, y:%.2f\n",b[0],b[1]); | |
1b45e564 | 1554 | |
644666df | 1555 | norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]); |
1556 | if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE; | |
1557 | ||
1558 | position[0]=b[0]+b[0]*TMath::Abs(r)/norm; | |
1559 | position[1]=b[1]+b[1]*TMath::Abs(r)/norm; | |
1560 | position[2]=0.; | |
1561 | return kTRUE; | |
f85a3764 | 1562 | } |
912f7fdc | 1563 | |
1564 | Double_t AliTPCPIDResponse::EvaldEdxSpline(Double_t bg,Int_t entry){ | |
1565 | // | |
1566 | // Evaluate the dEdx response for given entry | |
1567 | // | |
1568 | TSpline * spline = (TSpline*)fSplineArray.At(entry); | |
1569 | if (spline) return spline->Eval(bg); | |
1570 | return 0; | |
1571 | } | |
1572 | ||
1573 | ||
1574 | Bool_t AliTPCPIDResponse::RegisterSpline(const char * name, Int_t index){ | |
1575 | // | |
1576 | // register spline to be used for drawing comparisons | |
1577 | // | |
1578 | TFile * fTPCBB = TFile::Open("$ALICE_ROOT/OADB/COMMON/PID/data/TPCPIDResponse.root"); | |
1579 | TObjArray *arrayTPCPID= (TObjArray*) fTPCBB->Get("TPCPIDResponse"); | |
1580 | if (fSplineArray.GetEntriesFast()<index) fSplineArray.Expand(index*2); | |
1581 | TSpline3 *spline=0; | |
1582 | if (arrayTPCPID){ | |
1583 | spline = (TSpline3*)arrayTPCPID->FindObject(name); | |
1584 | if (spline) fSplineArray.AddAt(spline->Clone(),index); | |
1585 | } | |
1586 | delete arrayTPCPID; | |
1587 | delete fTPCBB; | |
1588 | return (spline!=0); | |
1589 | } |