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