]>
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), | |
74 | fSigmaPar0(0.0) | |
10d100d4 | 75 | { |
76 | // | |
77 | // The default constructor | |
78 | // | |
644666df | 79 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=0.07;fResN2[i]=0.0;} |
10d100d4 | 80 | } |
1b45e564 | 81 | /*TODO remove? |
10d100d4 | 82 | //_________________________________________________________________________ |
d2aa6df0 | 83 | AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param): |
644666df | 84 | TNamed(), |
d2aa6df0 | 85 | fMIP(param[0]), |
644666df | 86 | fRes0(), |
87 | fResN2(), | |
d2aa6df0 | 88 | fKp1(0.0283086), |
89 | fKp2(2.63394e+01), | |
90 | fKp3(5.04114e-11), | |
91 | fKp4(2.12543), | |
92 | fKp5(4.88663), | |
93 | fUseDatabase(kFALSE), | |
644666df | 94 | fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios), |
95 | fVoltageMap(72), | |
96 | fLowGainIROCthreshold(-40), | |
97 | fBadIROCthreshhold(-70), | |
98 | fLowGainOROCthreshold(-40), | |
99 | fBadOROCthreshhold(-40), | |
100 | fMaxBadLengthFraction(0.5), | |
f84b18dd | 101 | fMagField(0.), |
102 | fhEtaCorr(0x0), | |
103 | fhEtaSigmaPar1(0x0), | |
104 | fSigmaPar0(0.0) | |
10d100d4 | 105 | { |
106 | // | |
107 | // The main constructor | |
108 | // | |
644666df | 109 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=param[1];fResN2[i]=param[2];} |
110 | } | |
1b45e564 | 111 | */ |
f84b18dd | 112 | |
113 | //_________________________________________________________________________ | |
114 | AliTPCPIDResponse::~AliTPCPIDResponse() | |
115 | { | |
116 | // | |
117 | // Destructor | |
118 | // | |
119 | ||
120 | delete fhEtaCorr; | |
121 | fhEtaCorr = 0x0; | |
122 | ||
123 | delete fhEtaSigmaPar1; | |
124 | fhEtaSigmaPar1 = 0x0; | |
125 | } | |
126 | ||
127 | ||
644666df | 128 | //_________________________________________________________________________ |
129 | AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& that): | |
130 | TNamed(that), | |
131 | fMIP(that.fMIP), | |
132 | fRes0(), | |
133 | fResN2(), | |
134 | fKp1(that.fKp1), | |
135 | fKp2(that.fKp2), | |
136 | fKp3(that.fKp3), | |
137 | fKp4(that.fKp4), | |
138 | fKp5(that.fKp5), | |
139 | fUseDatabase(that.fUseDatabase), | |
140 | fResponseFunctions(that.fResponseFunctions), | |
141 | fVoltageMap(that.fVoltageMap), | |
142 | fLowGainIROCthreshold(that.fLowGainIROCthreshold), | |
143 | fBadIROCthreshhold(that.fBadIROCthreshhold), | |
144 | fLowGainOROCthreshold(that.fLowGainOROCthreshold), | |
145 | fBadOROCthreshhold(that.fBadOROCthreshhold), | |
146 | fMaxBadLengthFraction(that.fMaxBadLengthFraction), | |
f84b18dd | 147 | fMagField(that.fMagField), |
148 | fhEtaCorr(0x0), | |
149 | fhEtaSigmaPar1(0x0), | |
150 | fSigmaPar0(that.fSigmaPar0) | |
644666df | 151 | { |
152 | //copy ctor | |
153 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];} | |
f84b18dd | 154 | |
155 | // Copy eta maps | |
1b45e564 | 156 | if (that.fhEtaCorr) { |
3800f96f | 157 | fhEtaCorr = new TH2D(*(that.fhEtaCorr)); |
158 | fhEtaCorr->SetDirectory(0); | |
159 | } | |
f84b18dd | 160 | |
1b45e564 | 161 | if (that.fhEtaSigmaPar1) { |
3800f96f | 162 | fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1)); |
163 | fhEtaSigmaPar1->SetDirectory(0); | |
164 | } | |
644666df | 165 | } |
166 | ||
167 | //_________________________________________________________________________ | |
168 | AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that) | |
169 | { | |
170 | //assignment | |
171 | if (&that==this) return *this; | |
172 | TNamed::operator=(that); | |
173 | fMIP=that.fMIP; | |
174 | fKp1=that.fKp1; | |
175 | fKp2=that.fKp2; | |
176 | fKp3=that.fKp3; | |
177 | fKp4=that.fKp4; | |
178 | fKp5=that.fKp5; | |
179 | fUseDatabase=that.fUseDatabase; | |
180 | fResponseFunctions=that.fResponseFunctions; | |
181 | fVoltageMap=that.fVoltageMap; | |
182 | fLowGainIROCthreshold=that.fLowGainIROCthreshold; | |
183 | fBadIROCthreshhold=that.fBadIROCthreshhold; | |
184 | fLowGainOROCthreshold=that.fLowGainOROCthreshold; | |
185 | fBadOROCthreshhold=that.fBadOROCthreshhold; | |
186 | fMaxBadLengthFraction=that.fMaxBadLengthFraction; | |
187 | fMagField=that.fMagField; | |
188 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];} | |
f84b18dd | 189 | |
190 | delete fhEtaCorr; | |
f85a3764 | 191 | fhEtaCorr=0x0; |
1b45e564 | 192 | if (that.fhEtaCorr) { |
3800f96f | 193 | fhEtaCorr = new TH2D(*(that.fhEtaCorr)); |
194 | fhEtaCorr->SetDirectory(0); | |
195 | } | |
f84b18dd | 196 | |
197 | delete fhEtaSigmaPar1; | |
f85a3764 | 198 | fhEtaSigmaPar1=0x0; |
1b45e564 | 199 | if (that.fhEtaSigmaPar1) { |
3800f96f | 200 | fhEtaSigmaPar1 = new TH2D(*(that.fhEtaSigmaPar1)); |
201 | fhEtaSigmaPar1->SetDirectory(0); | |
202 | } | |
f84b18dd | 203 | |
204 | fSigmaPar0 = that.fSigmaPar0; | |
205 | ||
644666df | 206 | return *this; |
10d100d4 | 207 | } |
208 | ||
d2aa6df0 | 209 | //_________________________________________________________________________ |
10d100d4 | 210 | Double_t AliTPCPIDResponse::Bethe(Double_t betaGamma) const { |
211 | // | |
212 | // This is the Bethe-Bloch function normalised to 1 at the minimum | |
213 | // WARNING | |
214 | // Simulated and reconstructed Bethe-Bloch differs | |
215 | // Simulated curve is the dNprim/dx | |
216 | // Reconstructed is proportianal dNtot/dx | |
217 | // Temporary fix for production - Simple linear correction function | |
218 | // Future 2 Bethe Bloch formulas needed | |
219 | // 1. for simulation | |
220 | // 2. for reconstructed PID | |
221 | // | |
d2aa6df0 | 222 | |
223 | // const Float_t kmeanCorrection =0.1; | |
10d100d4 | 224 | Double_t bb= |
225 | AliExternalTrackParam::BetheBlochAleph(betaGamma,fKp1,fKp2,fKp3,fKp4,fKp5); | |
226 | return bb*fMIP; | |
227 | } | |
228 | ||
229 | //_________________________________________________________________________ | |
230 | void AliTPCPIDResponse::SetBetheBlochParameters(Double_t kp1, | |
231 | Double_t kp2, | |
232 | Double_t kp3, | |
233 | Double_t kp4, | |
234 | Double_t kp5) { | |
235 | // | |
236 | // Set the parameters of the ALEPH Bethe-Bloch formula | |
237 | // | |
238 | fKp1=kp1; | |
239 | fKp2=kp2; | |
240 | fKp3=kp3; | |
241 | fKp4=kp4; | |
242 | fKp5=kp5; | |
243 | } | |
644666df | 244 | |
10d100d4 | 245 | //_________________________________________________________________________ |
246 | void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2) { | |
247 | // | |
248 | // Set the relative resolution sigma_rel = res0 * sqrt(1+resN2/npoint) | |
249 | // | |
644666df | 250 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=res0;fResN2[i]=resN2;} |
10d100d4 | 251 | } |
252 | ||
253 | //_________________________________________________________________________ | |
254 | Double_t AliTPCPIDResponse::GetExpectedSignal(const Float_t mom, | |
d2aa6df0 | 255 | AliPID::EParticleType n) const { |
10d100d4 | 256 | // |
f84b18dd | 257 | // Deprecated function (for backward compatibility). Please use |
258 | // GetExpectedSignal(const AliVTrack* track, AliPID::EParticleType species, ETPCdEdxSource dedxSource, | |
259 | // Bool_t correctEta = kTRUE); | |
260 | // instead! | |
261 | // | |
262 | // | |
10d100d4 | 263 | // Calculates the expected PID signal as the function of |
264 | // the information stored in the track, for the specified particle type | |
265 | // | |
266 | // At the moment, these signals are just the results of calling the | |
267 | // Bethe-Bloch formula. | |
268 | // This can be improved. By taking into account the number of | |
269 | // assigned clusters and/or the track dip angle, for example. | |
270 | // | |
271 | ||
f84b18dd | 272 | //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...) |
273 | // !!! Splines for light nuclei need to be normalised to this factor !!! | |
274 | const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(n),2.3); | |
275 | ||
539a5a59 | 276 | Double_t mass=AliPID::ParticleMassZ(n); |
f84b18dd | 277 | if (!fUseDatabase) return Bethe(mom/mass) * chargeFactor; |
d2aa6df0 | 278 | // |
644666df | 279 | const TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n); |
00a38d07 | 280 | |
f84b18dd | 281 | if (!responseFunction) return Bethe(mom/mass) * chargeFactor; |
282 | ||
00a38d07 | 283 | return fMIP*responseFunction->Eval(mom/mass)*chargeFactor; |
d2aa6df0 | 284 | |
10d100d4 | 285 | } |
286 | ||
287 | //_________________________________________________________________________ | |
288 | Double_t AliTPCPIDResponse::GetExpectedSigma(const Float_t mom, | |
f84b18dd | 289 | const Int_t nPoints, |
644666df | 290 | AliPID::EParticleType n) const { |
10d100d4 | 291 | // |
f84b18dd | 292 | // Deprecated function (for backward compatibility). Please use |
293 | // GetExpectedSigma(onst AliVTrack* track, AliPID::EParticleType species, | |
294 | // ETPCdEdxSource dedxSource, Bool_t correctEta) instead! | |
295 | // | |
296 | // | |
10d100d4 | 297 | // Calculates the expected sigma of the PID signal as the function of |
298 | // the information stored in the track, for the specified particle type | |
299 | // | |
644666df | 300 | |
301 | if (nPoints != 0) | |
302 | return GetExpectedSignal(mom,n)*fRes0[0]*sqrt(1. + fResN2[0]/nPoints); | |
303 | else | |
304 | return GetExpectedSignal(mom,n)*fRes0[0]; | |
305 | } | |
306 | ||
307 | ////////////////////////////////////////////////////NEW////////////////////////////// | |
308 | ||
309 | //_________________________________________________________________________ | |
310 | void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario) { | |
311 | // | |
312 | // Set the relative resolution sigma_rel = res0 * sqrt(1+resN2/npoint) | |
313 | // | |
314 | if ((Int_t)gainScenario>(Int_t)fgkNumberOfGainScenarios) return; //TODO: better handling! | |
315 | fRes0[gainScenario]=res0; | |
316 | fResN2[gainScenario]=resN2; | |
317 | } | |
318 | ||
f84b18dd | 319 | |
644666df | 320 | //_________________________________________________________________________ |
f84b18dd | 321 | Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track, |
644666df | 322 | AliPID::EParticleType species, |
f84b18dd | 323 | Double_t /*dEdx*/, |
324 | const TSpline3* responseFunction, | |
325 | Bool_t correctEta) const | |
644666df | 326 | { |
327 | // Calculates the expected PID signal as the function of | |
f84b18dd | 328 | // the information stored in the track and the given parameters, |
329 | // for the specified particle type | |
644666df | 330 | // |
331 | // At the moment, these signals are just the results of calling the | |
f84b18dd | 332 | // Bethe-Bloch formula plus, if desired, taking into account the eta dependence. |
644666df | 333 | // This can be improved. By taking into account the number of |
334 | // assigned clusters and/or the track dip angle, for example. | |
10d100d4 | 335 | // |
336 | ||
f84b18dd | 337 | Double_t mom=track->GetTPCmomentum(); |
97d17530 | 338 | Double_t mass=AliPID::ParticleMassZ(species); |
644666df | 339 | |
f84b18dd | 340 | //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...) |
341 | // !!! Splines for light nuclei need to be normalised to this factor !!! | |
342 | const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(species),2.3); | |
343 | ||
344 | if (!responseFunction) | |
345 | return Bethe(mom/mass) * chargeFactor; | |
644666df | 346 | |
f84b18dd | 347 | Double_t dEdxSplines = fMIP*responseFunction->Eval(mom/mass) * chargeFactor; |
348 | ||
349 | if (!correctEta) | |
350 | return dEdxSplines; | |
351 | ||
352 | //TODO Alternatively take current track dEdx | |
353 | //return dEdxSplines * GetEtaCorrection(track, dEdx); | |
354 | return dEdxSplines * GetEtaCorrection(track, dEdxSplines); | |
644666df | 355 | } |
356 | ||
f84b18dd | 357 | |
644666df | 358 | //_________________________________________________________________________ |
359 | Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track, | |
360 | AliPID::EParticleType species, | |
f84b18dd | 361 | ETPCdEdxSource dedxSource, |
362 | Bool_t correctEta) const | |
644666df | 363 | { |
364 | // Calculates the expected PID signal as the function of | |
365 | // the information stored in the track, for the specified particle type | |
366 | // | |
367 | // At the moment, these signals are just the results of calling the | |
f84b18dd | 368 | // Bethe-Bloch formula plus, if desired, taking into account the eta dependence. |
644666df | 369 | // This can be improved. By taking into account the number of |
370 | // assigned clusters and/or the track dip angle, for example. | |
371 | // | |
372 | ||
f84b18dd | 373 | if (!fUseDatabase) { |
374 | //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...) | |
375 | // !!! Splines for light nuclei need to be normalised to this factor !!! | |
376 | const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(species),2.3); | |
644666df | 377 | |
97d17530 | 378 | return Bethe(track->GetTPCmomentum() / AliPID::ParticleMassZ(species)) * chargeFactor; |
f84b18dd | 379 | } |
644666df | 380 | |
f84b18dd | 381 | Double_t dEdx = -1; |
382 | Int_t nPoints = -1; | |
383 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
384 | TSpline3* responseFunction = 0x0; | |
385 | ||
386 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) { | |
387 | // Something is wrong with the track -> Return obviously invalid value | |
388 | return -999; | |
389 | } | |
390 | ||
391 | // Charge factor already taken into account inside the following function call | |
392 | return GetExpectedSignal(track, species, dEdx, responseFunction, correctEta); | |
644666df | 393 | } |
394 | ||
395 | //_________________________________________________________________________ | |
396 | TSpline3* AliTPCPIDResponse::GetResponseFunction( AliPID::EParticleType type, | |
397 | AliTPCPIDResponse::ETPCgainScenario gainScenario ) const | |
398 | { | |
399 | //get response function | |
400 | return dynamic_cast<TSpline3*>(fResponseFunctions.At(ResponseFunctionIndex(type,gainScenario))); | |
401 | } | |
402 | ||
403 | //_________________________________________________________________________ | |
404 | TSpline3* AliTPCPIDResponse::GetResponseFunction( const AliVTrack* track, | |
405 | AliPID::EParticleType species, | |
f84b18dd | 406 | ETPCdEdxSource dedxSource) const |
644666df | 407 | { |
408 | //the splines are stored in an array, different scenarios | |
409 | ||
f84b18dd | 410 | Double_t dEdx = -1; |
411 | Int_t nPoints = -1; | |
412 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
413 | TSpline3* responseFunction = 0x0; | |
414 | ||
415 | if (ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
416 | return responseFunction; | |
644666df | 417 | |
418 | return NULL; | |
419 | } | |
420 | ||
421 | //_________________________________________________________________________ | |
422 | void AliTPCPIDResponse::ResetSplines() | |
423 | { | |
424 | //reset the array with splines | |
425 | for (Int_t i=0;i<fResponseFunctions.GetEntriesFast();i++) | |
426 | { | |
427 | fResponseFunctions.AddAt(NULL,i); | |
428 | } | |
429 | } | |
430 | //_________________________________________________________________________ | |
431 | Int_t AliTPCPIDResponse::ResponseFunctionIndex( AliPID::EParticleType species, | |
432 | ETPCgainScenario gainScenario ) const | |
433 | { | |
434 | //get the index in fResponseFunctions given type and scenario | |
435 | return Int_t(species)+Int_t(gainScenario)*fgkNumberOfParticleSpecies; | |
436 | } | |
437 | ||
438 | //_________________________________________________________________________ | |
439 | void AliTPCPIDResponse::SetResponseFunction( TObject* o, | |
440 | AliPID::EParticleType species, | |
441 | ETPCgainScenario gainScenario ) | |
442 | { | |
443 | fResponseFunctions.AddAtAndExpand(o,ResponseFunctionIndex(species,gainScenario)); | |
444 | } | |
445 | ||
f84b18dd | 446 | |
644666df | 447 | //_________________________________________________________________________ |
f84b18dd | 448 | Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, |
449 | AliPID::EParticleType species, | |
450 | ETPCgainScenario gainScenario, | |
451 | Double_t dEdx, | |
452 | Int_t nPoints, | |
453 | const TSpline3* responseFunction, | |
454 | Bool_t correctEta) const | |
644666df | 455 | { |
456 | // Calculates the expected sigma of the PID signal as the function of | |
f84b18dd | 457 | // the information stored in the track and the given parameters, |
458 | // for the specified particle type | |
644666df | 459 | // |
460 | ||
f84b18dd | 461 | if (!responseFunction) |
462 | return 999; | |
463 | ||
464 | // If no sigma map is available or if no eta correction is requested (sigma maps only for corrected eta!), use the old parametrisation | |
465 | if (!fhEtaSigmaPar1 || !correctEta) { | |
466 | if (nPoints != 0) | |
467 | return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE) * | |
468 | fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints); | |
469 | else | |
470 | return GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE)*fRes0[gainScenario]; | |
471 | } | |
472 | ||
473 | if (nPoints > 0) { | |
474 | Double_t sigmaPar1 = GetSigmaPar1(track, species, dEdx, responseFunction); | |
475 | ||
476 | return GetExpectedSignal(track, species, dEdx, responseFunction, kTRUE)*sqrt(fSigmaPar0 * fSigmaPar0 + sigmaPar1 * sigmaPar1 / nPoints); | |
477 | } | |
478 | else { | |
479 | // One should never have/take tracks with 0 dEdx clusters! | |
480 | return 999; | |
481 | } | |
644666df | 482 | } |
483 | ||
f84b18dd | 484 | |
644666df | 485 | //_________________________________________________________________________ |
486 | Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, | |
487 | AliPID::EParticleType species, | |
f84b18dd | 488 | ETPCdEdxSource dedxSource, |
489 | Bool_t correctEta) const | |
644666df | 490 | { |
491 | // Calculates the expected sigma of the PID signal as the function of | |
492 | // the information stored in the track, for the specified particle type | |
493 | // and dedx scenario | |
494 | // | |
495 | ||
f84b18dd | 496 | Double_t dEdx = -1; |
497 | Int_t nPoints = -1; | |
498 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
499 | TSpline3* responseFunction = 0x0; | |
500 | ||
501 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
502 | return 999; //TODO: better handling! | |
503 | ||
504 | return GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta); | |
644666df | 505 | } |
506 | ||
f84b18dd | 507 | |
644666df | 508 | //_________________________________________________________________________ |
509 | Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track, | |
510 | AliPID::EParticleType species, | |
f84b18dd | 511 | ETPCdEdxSource dedxSource, |
512 | Bool_t correctEta) const | |
644666df | 513 | { |
f84b18dd | 514 | //Calculates the number of sigmas of the PID signal from the expected value |
515 | //for a given particle species in the presence of multiple gain scenarios | |
644666df | 516 | //inside the TPC |
f84b18dd | 517 | |
518 | Double_t dEdx = -1; | |
519 | Int_t nPoints = -1; | |
520 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
521 | TSpline3* responseFunction = 0x0; | |
522 | ||
523 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
524 | return -999; //TODO: Better handling! | |
525 | ||
526 | Double_t bethe = GetExpectedSignal(track, species, dEdx, responseFunction, correctEta); | |
527 | Double_t sigma = GetExpectedSigma(track, species, gainScenario, dEdx, nPoints, responseFunction, correctEta); | |
528 | // 999 will be returned by GetExpectedSigma e.g. in case of 0 dEdx clusters | |
529 | if (sigma >= 998) | |
530 | return -999; | |
531 | else | |
532 | return (dEdx-bethe)/sigma; | |
644666df | 533 | } |
534 | ||
535 | //_________________________________________________________________________ | |
536 | Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track, | |
537 | AliPID::EParticleType species, | |
f84b18dd | 538 | ETPCdEdxSource dedxSource, |
1b45e564 | 539 | Double_t& dEdx, |
540 | Int_t& nPoints, | |
541 | ETPCgainScenario& gainScenario, | |
f85a3764 | 542 | TSpline3** responseFunction) const |
644666df | 543 | { |
544 | // Calculates the right parameters for PID | |
545 | // dEdx parametrization for the proper gain scenario, dEdx | |
546 | // and NPoints used for dEdx | |
547 | // based on the track geometry (which chambers it crosses) for the specified particle type | |
548 | // and preferred source of dedx. | |
549 | // returns true on success | |
550 | ||
f84b18dd | 551 | |
552 | if (dedxSource == kdEdxDefault) { | |
553 | // Fast handling for default case. In addition: Keep it simple (don't call additional functions) to | |
554 | // avoid possible bugs | |
f85a3764 | 555 | |
556 | // GetTPCsignalTunedOnData will be non-positive, if it has not been set (i.e. in case of MC NOT tuned to data). | |
557 | // If this is the case, just take the normal signal | |
558 | dEdx = track->GetTPCsignalTunedOnData(); | |
559 | if (dEdx <= 0) { | |
560 | dEdx = track->GetTPCsignal(); | |
561 | } | |
562 | ||
f84b18dd | 563 | nPoints = track->GetTPCsignalN(); |
564 | gainScenario = kDefault; | |
565 | ||
566 | TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario)); | |
567 | *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast? | |
568 | ||
569 | return kTRUE; | |
570 | } | |
571 | ||
f85a3764 | 572 | //TODO Proper handle of tuneMConData for other dEdx sources |
f84b18dd | 573 | |
644666df | 574 | Double32_t signal[4]; //0: IROC, 1: OROC medium, 2:OROC long, 3: OROC all (def. truncation used) |
575 | Char_t ncl[3]; //same | |
576 | Char_t nrows[3]; //same | |
f84b18dd | 577 | const AliTPCdEdxInfo* dEdxInfo = track->GetTPCdEdxInfo(); |
644666df | 578 | |
579 | if (!dEdxInfo && dedxSource!=kdEdxDefault) //in one case its ok if we dont have the info | |
580 | { | |
581 | AliError("AliTPCdEdxInfo not available"); | |
644666df | 582 | return kFALSE; |
583 | } | |
584 | ||
585 | if (dEdxInfo) dEdxInfo->GetTPCSignalRegionInfo(signal,ncl,nrows); | |
586 | ||
587 | //check if we cross a bad OROC in which case we reject | |
f85a3764 | 588 | EChamberStatus trackOROCStatus = TrackStatus(track,2); |
589 | if (trackOROCStatus==kChamberOff || trackOROCStatus==kChamberLowGain) | |
644666df | 590 | { |
644666df | 591 | return kFALSE; |
592 | } | |
593 | ||
594 | switch (dedxSource) | |
595 | { | |
644666df | 596 | case kdEdxOROC: |
597 | { | |
f85a3764 | 598 | if (trackOROCStatus==kChamberInvalid) return kFALSE; //never reached OROC |
f84b18dd | 599 | dEdx = signal[3]; |
600 | nPoints = ncl[2]+ncl[1]; | |
601 | gainScenario = kOROChigh; | |
644666df | 602 | break; |
603 | } | |
604 | case kdEdxHybrid: | |
605 | { | |
606 | //if we cross a bad IROC we use OROC dedx, if we dont we use combined | |
607 | EChamberStatus status = TrackStatus(track,1); | |
608 | if (status!=kChamberHighGain) | |
609 | { | |
f84b18dd | 610 | dEdx = signal[3]; |
611 | nPoints = ncl[2]+ncl[1]; | |
612 | gainScenario = kOROChigh; | |
644666df | 613 | } |
614 | else | |
615 | { | |
f84b18dd | 616 | dEdx = track->GetTPCsignal(); |
617 | nPoints = track->GetTPCsignalN(); | |
618 | gainScenario = kALLhigh; | |
644666df | 619 | } |
620 | break; | |
621 | } | |
622 | default: | |
623 | { | |
f84b18dd | 624 | dEdx = 0.; |
625 | nPoints = 0; | |
626 | gainScenario = kGainScenarioInvalid; | |
644666df | 627 | return kFALSE; |
628 | } | |
629 | } | |
f84b18dd | 630 | TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,gainScenario)); |
631 | *responseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast? | |
632 | ||
644666df | 633 | return kTRUE; |
634 | } | |
635 | ||
f84b18dd | 636 | |
637 | //_________________________________________________________________________ | |
638 | Double_t AliTPCPIDResponse::GetEtaCorrection(const AliVTrack *track, Double_t dEdxSplines) const | |
639 | { | |
640 | // | |
641 | // Get eta correction for the given parameters. | |
642 | // | |
643 | ||
644 | if (!fhEtaCorr) | |
645 | return 1.; | |
646 | ||
647 | Double_t tpcSignal = dEdxSplines; | |
648 | ||
649 | if (tpcSignal < 1.) | |
650 | return 1.; | |
651 | ||
652 | // For ESD tracks, the local tanTheta could be used (esdTrack->GetInnerParam()->GetTgl()). | |
653 | // However, this value is not available for AODs and, thus, not for AliVTrack. | |
654 | // Fortunately, the following formula allows to approximate the local tanTheta with the | |
655 | // global theta angle -> This is for by far most of the tracks the same, but gives at | |
656 | // maybe the percent level differences within +- 0.2 in tanTheta -> Which is still ok. | |
657 | Double_t tanTheta = TMath::Tan(-track->Theta() + TMath::Pi() / 2.0); | |
658 | Int_t binX = fhEtaCorr->GetXaxis()->FindBin(tanTheta); | |
659 | Int_t binY = fhEtaCorr->GetYaxis()->FindBin(1. / tpcSignal); | |
660 | ||
661 | if (binX == 0) | |
662 | binX = 1; | |
663 | if (binX > fhEtaCorr->GetXaxis()->GetNbins()) | |
664 | binX = fhEtaCorr->GetXaxis()->GetNbins(); | |
665 | ||
666 | if (binY == 0) | |
667 | binY = 1; | |
668 | if (binY > fhEtaCorr->GetYaxis()->GetNbins()) | |
669 | binY = fhEtaCorr->GetYaxis()->GetNbins(); | |
670 | ||
671 | return fhEtaCorr->GetBinContent(binX, binY); | |
672 | } | |
673 | ||
674 | ||
675 | //_________________________________________________________________________ | |
676 | Double_t AliTPCPIDResponse::GetEtaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const | |
677 | { | |
678 | // | |
679 | // Get eta correction for the given track. | |
680 | // | |
681 | ||
682 | if (!fhEtaCorr) | |
683 | return 1.; | |
684 | ||
685 | Double_t dEdx = -1; | |
686 | Int_t nPoints = -1; | |
687 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
688 | TSpline3* responseFunction = 0x0; | |
689 | ||
690 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
691 | return 1.; | |
692 | ||
693 | Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE); | |
694 | ||
695 | //TODO Alternatively take current track dEdx | |
696 | //return GetEtaCorrection(track, dEdx); | |
697 | ||
698 | return GetEtaCorrection(track, dEdxSplines); | |
699 | } | |
700 | ||
701 | ||
702 | //_________________________________________________________________________ | |
703 | Double_t AliTPCPIDResponse::GetEtaCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const | |
704 | { | |
705 | // | |
706 | // Get eta corrected dEdx for the given track. For the correction, the expected dEdx of | |
707 | // the specified species will be used. If the species is set to AliPID::kUnknown, the | |
708 | // dEdx of the track is used instead. | |
709 | // WARNING: In the latter case, the eta correction might not be as good as if the | |
710 | // expected dEdx is used, which is the way the correction factor is designed | |
711 | // for. | |
712 | // In any case, one has to decide carefully to which expected signal one wants to | |
713 | // compare the corrected value - to the corrected or uncorrected. | |
714 | // Anyhow, a safer way of looking e.g. at the n-sigma is to call | |
715 | // the corresponding function GetNumberOfSigmas! | |
716 | // | |
717 | ||
718 | Double_t dEdx = -1; | |
719 | Int_t nPoints = -1; | |
720 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
721 | TSpline3* responseFunction = 0x0; | |
722 | ||
723 | // Note: In case of species == AliPID::kUnknown, the responseFunction might not be set. However, in this case | |
724 | // it is not used anyway, so this causes no trouble. | |
725 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
726 | return -1.; | |
727 | ||
728 | Double_t etaCorr = 0.; | |
729 | ||
730 | if (species < AliPID::kUnknown) { | |
731 | Double_t dEdxSplines = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE); | |
732 | etaCorr = GetEtaCorrection(track, dEdxSplines); | |
733 | } | |
734 | else { | |
735 | etaCorr = GetEtaCorrection(track, dEdx); | |
736 | } | |
737 | ||
738 | if (etaCorr <= 0) | |
739 | return -1.; | |
740 | ||
741 | return dEdx / etaCorr; | |
742 | } | |
743 | ||
744 | ||
745 | ||
746 | //_________________________________________________________________________ | |
747 | Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, Double_t dEdx, const TSpline3* responseFunction) const | |
748 | { | |
749 | // | |
750 | // Get parameter 1 of sigma parametrisation of TPC dEdx from the histogram for the given track. | |
751 | // | |
752 | ||
753 | if (!fhEtaSigmaPar1) | |
754 | return 999; | |
755 | ||
756 | // The sigma maps are created with data sets that are already eta corrected and for which the | |
757 | // splines have been re-created. Therefore, the value for the lookup needs to be the value of | |
758 | // the splines without any additional eta correction. | |
759 | // NOTE: This is due to the method the maps are created. The track dEdx (not the expected one!) | |
760 | // is corrected to uniquely related a momemtum bin with an expected dEdx, where the expected dEdx | |
761 | // equals the track dEdx for all eta (thanks to the correction and the re-fit of the splines). | |
762 | // Consequently, looking up the uncorrected expected dEdx at a given tanTheta yields the correct | |
763 | // sigma parameter! | |
764 | // Also: It has to be the spline dEdx, since one wants to get the sigma for the assumption(!) | |
765 | // of such a particle, which by assumption then has this dEdx value | |
766 | ||
767 | Double_t dEdxExpected = GetExpectedSignal(track, species, dEdx, responseFunction, kFALSE); | |
768 | ||
769 | if (dEdxExpected < 1.) | |
770 | return 999; | |
771 | ||
772 | // For ESD tracks, the local tanTheta could be used (esdTrack->GetInnerParam()->GetTgl()). | |
773 | // However, this value is not available for AODs and, thus, not or AliVTrack. | |
774 | // Fortunately, the following formula allows to approximate the local tanTheta with the | |
775 | // global theta angle -> This is for by far most of the tracks the same, but gives at | |
776 | // maybe the percent level differences within +- 0.2 in tanTheta -> Which is still ok. | |
777 | Double_t tanTheta = TMath::Tan(-track->Theta() + TMath::Pi() / 2.0); | |
778 | Int_t binX = fhEtaSigmaPar1->GetXaxis()->FindBin(tanTheta); | |
779 | Int_t binY = fhEtaSigmaPar1->GetYaxis()->FindBin(1. / dEdxExpected); | |
780 | ||
781 | if (binX == 0) | |
782 | binX = 1; | |
783 | if (binX > fhEtaSigmaPar1->GetXaxis()->GetNbins()) | |
784 | binX = fhEtaSigmaPar1->GetXaxis()->GetNbins(); | |
785 | ||
786 | if (binY == 0) | |
787 | binY = 1; | |
788 | if (binY > fhEtaSigmaPar1->GetYaxis()->GetNbins()) | |
789 | binY = fhEtaSigmaPar1->GetYaxis()->GetNbins(); | |
790 | ||
791 | return fhEtaSigmaPar1->GetBinContent(binX, binY); | |
792 | } | |
793 | ||
794 | ||
795 | //_________________________________________________________________________ | |
796 | Double_t AliTPCPIDResponse::GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource) const | |
797 | { | |
798 | // | |
799 | // Get eta correction for the given track. | |
800 | // | |
801 | ||
802 | if (!fhEtaSigmaPar1) | |
803 | return 999; | |
804 | ||
805 | Double_t dEdx = -1; | |
806 | Int_t nPoints = -1; | |
807 | ETPCgainScenario gainScenario = kGainScenarioInvalid; | |
808 | TSpline3* responseFunction = 0x0; | |
809 | ||
810 | if (!ResponseFunctiondEdxN(track, species, dedxSource, dEdx, nPoints, gainScenario, &responseFunction)) | |
811 | return 999; | |
812 | ||
813 | return GetSigmaPar1(track, species, dEdx, responseFunction); | |
814 | } | |
815 | ||
816 | ||
817 | //_________________________________________________________________________ | |
818 | Bool_t AliTPCPIDResponse::SetEtaCorrMap(TH2D* hMap) | |
819 | { | |
820 | // | |
821 | // Load map for TPC eta correction (a copy is stored and will be deleted automatically). | |
822 | // If hMap is 0x0,the eta correction will be disabled and kFALSE is returned. | |
823 | // If the map can be set, kTRUE is returned. | |
824 | // | |
825 | ||
826 | delete fhEtaCorr; | |
827 | ||
828 | if (!hMap) { | |
829 | fhEtaCorr = 0x0; | |
830 | ||
831 | return kFALSE; | |
832 | } | |
833 | ||
834 | fhEtaCorr = (TH2D*)(hMap->Clone()); | |
835 | fhEtaCorr->SetDirectory(0); | |
836 | ||
837 | return kTRUE; | |
838 | } | |
839 | ||
840 | //_________________________________________________________________________ | |
841 | Bool_t AliTPCPIDResponse::SetSigmaParams(TH2D* hSigmaPar1Map, Double_t sigmaPar0) | |
842 | { | |
843 | // | |
844 | // Load map for TPC sigma map (a copy is stored and will be deleted automatically): | |
845 | // Parameter 1 is stored as a 2D map (1/dEdx vs. tanTheta_local) and parameter 0 is | |
846 | // a constant. If hSigmaPar1Map is 0x0, the old sigma parametrisation will be used | |
847 | // (and sigmaPar0 is ignored!) and kFALSE is returned. | |
848 | // If the map can be set, sigmaPar0 is also set and kTRUE will be returned. | |
849 | // | |
850 | ||
851 | delete fhEtaSigmaPar1; | |
852 | ||
853 | if (!hSigmaPar1Map) { | |
854 | fhEtaSigmaPar1 = 0x0; | |
855 | fSigmaPar0 = 0.0; | |
856 | ||
857 | return kFALSE; | |
858 | } | |
859 | ||
860 | fhEtaSigmaPar1 = (TH2D*)(hSigmaPar1Map->Clone()); | |
861 | fhEtaSigmaPar1->SetDirectory(0); | |
862 | fSigmaPar0 = sigmaPar0; | |
863 | ||
864 | return kTRUE; | |
865 | } | |
866 | ||
867 | ||
644666df | 868 | //_________________________________________________________________________ |
f85a3764 | 869 | Bool_t AliTPCPIDResponse::sectorNumbersInOut(Double_t* trackPositionInner, |
870 | Double_t* trackPositionOuter, | |
871 | Float_t& inphi, | |
644666df | 872 | Float_t& outphi, |
f85a3764 | 873 | Int_t& in, |
1b45e564 | 874 | Int_t& out ) const |
644666df | 875 | { |
876 | //calculate the sector numbers (equivalent to IROC chamber numbers) a track crosses | |
877 | //for OROC chamber numbers add 36 | |
878 | //returned angles are between (0,2pi) | |
644666df | 879 | |
880 | inphi = TMath::ATan2(trackPositionInner[1],trackPositionInner[0]); | |
881 | outphi = TMath::ATan2(trackPositionOuter[1], trackPositionOuter[0]); | |
882 | ||
883 | if (inphi<0) {inphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi | |
884 | if (outphi<0) {outphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi | |
885 | ||
886 | in = sectorNumber(inphi); | |
887 | out = sectorNumber(outphi); | |
888 | ||
889 | //for the C side (positive z) offset by 18 | |
890 | if (trackPositionInner[2]>0.0) | |
891 | { | |
892 | in+=18; | |
893 | out+=18; | |
894 | } | |
895 | return kTRUE; | |
896 | } | |
1b45e564 | 897 | |
898 | ||
644666df | 899 | //_____________________________________________________________________________ |
900 | Int_t AliTPCPIDResponse::sectorNumber(Double_t phi) const | |
901 | { | |
902 | //calculate sector number | |
903 | const Float_t width=TMath::TwoPi()/18.; | |
904 | return TMath::Floor(phi/width); | |
905 | } | |
906 | ||
1b45e564 | 907 | |
644666df | 908 | //_____________________________________________________________________________ |
909 | void AliTPCPIDResponse::Print(Option_t* /*option*/) const | |
910 | { | |
911 | //Print info | |
912 | fResponseFunctions.Print(); | |
913 | } | |
914 | ||
1b45e564 | 915 | |
644666df | 916 | //_____________________________________________________________________________ |
917 | AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack* track, Int_t layer) const | |
918 | { | |
919 | //status of the track: if it crosses any bad chambers return kChamberOff; | |
920 | //IROC:layer=1, OROC:layer=2 | |
921 | if (layer<1 || layer>2) layer=1; | |
922 | Int_t in=0; | |
923 | Int_t out=0; | |
924 | Float_t inphi=0.; | |
925 | Float_t outphi=0.; | |
926 | Float_t innerRadius = (layer==1)?83.0:133.7; | |
927 | Float_t outerRadius = (layer==1)?133.5:247.7; | |
f85a3764 | 928 | |
929 | ///////////////////////////////////////////////////////////////////////////// | |
930 | //find out where track enters and leaves the layer. | |
931 | // | |
1b45e564 | 932 | Double_t trackPositionInner[3]; |
933 | Double_t trackPositionOuter[3]; | |
934 | ||
f85a3764 | 935 | //if there is no inner param this could mean we're using the AOD track, |
936 | //we still can extrapolate from the vertex - so use those params. | |
937 | const AliExternalTrackParam* ip = track->GetInnerParam(); | |
938 | if (ip) track=ip; | |
939 | ||
940 | Bool_t trackAtInner = track->GetXYZAt(innerRadius, fMagField, trackPositionInner); | |
941 | Bool_t trackAtOuter = track->GetXYZAt(outerRadius, fMagField, trackPositionOuter); | |
942 | ||
943 | if (!trackAtInner) | |
944 | { | |
945 | //if we dont even enter inner radius we do nothing and return invalid | |
946 | inphi=0.0; | |
947 | outphi=0.0; | |
948 | in=0; | |
949 | out=0; | |
950 | return kChamberInvalid; | |
951 | } | |
952 | ||
953 | if (!trackAtOuter) | |
954 | { | |
955 | //if we don't reach the outer radius check that the apex is indeed within the outer radius and use apex position | |
956 | Bool_t haveApex = TrackApex(track, fMagField, trackPositionOuter); | |
957 | Float_t apexRadius = TMath::Sqrt(trackPositionOuter[0]*trackPositionOuter[0]+trackPositionOuter[1]*trackPositionOuter[1]); | |
958 | if ( haveApex && apexRadius<=outerRadius && apexRadius>innerRadius) | |
959 | { | |
960 | //printf("pt: %.2f, apexRadius: %.2f(%s), x: %.2f, y: %.2f\n",track->Pt(),apexRadius,(haveApex)?"OK":"BAD",trackPositionOuter[0],trackPositionOuter[1]); | |
961 | } | |
962 | else | |
963 | { | |
964 | inphi=0.0; | |
965 | outphi=0.0; | |
966 | in=0; | |
967 | out=0; | |
968 | return kChamberInvalid; | |
969 | } | |
970 | } | |
971 | ||
972 | ||
1b45e564 | 973 | if (!sectorNumbersInOut(trackPositionInner, |
974 | trackPositionOuter, | |
975 | inphi, | |
976 | outphi, | |
977 | in, | |
f85a3764 | 978 | out)) return kChamberInvalid; |
644666df | 979 | |
f85a3764 | 980 | ///////////////////////////////////////////////////////////////////////////// |
1b45e564 | 981 | //now we have the location of the track we can check |
f85a3764 | 982 | //if it is in a good/bad chamber |
983 | // | |
644666df | 984 | Bool_t sideA = kTRUE; |
1b45e564 | 985 | |
644666df | 986 | if (((in/18)%2==1) && ((out/18)%2==1)) sideA=kFALSE; |
987 | ||
988 | in=in%18; | |
989 | out=out%18; | |
990 | ||
991 | if (TMath::Abs(in-out)>9) | |
992 | { | |
993 | if (TMath::Max(in,out)==out) | |
994 | { | |
995 | Int_t tmp=out; | |
996 | out = in; | |
997 | in = tmp; | |
998 | Float_t tmpphi=outphi; | |
999 | outphi=inphi; | |
1000 | inphi=tmpphi; | |
1001 | } | |
1002 | in-=18; | |
1003 | inphi-=TMath::TwoPi(); | |
1004 | } | |
10d100d4 | 1005 | else |
644666df | 1006 | { |
1007 | if (TMath::Max(in,out)==in) | |
1008 | { | |
1009 | Int_t tmp=out; | |
1010 | out=in; | |
1011 | in=tmp; | |
1012 | Float_t tmpphi=outphi; | |
1013 | outphi=inphi; | |
1014 | inphi=tmpphi; | |
1015 | } | |
1016 | } | |
1017 | ||
1018 | Float_t trackLengthInBad = 0.; | |
1019 | Float_t trackLengthInLowGain = 0.; | |
1020 | Float_t trackLengthTotal = TMath::Abs(outphi-inphi); | |
f85a3764 | 1021 | Float_t lengthFractionInBadSectors = 0.; |
644666df | 1022 | |
1023 | const Float_t sectorWidth = TMath::TwoPi()/18.; | |
1b45e564 | 1024 | |
644666df | 1025 | for (Int_t i=in; i<=out; i++) |
1026 | { | |
1027 | int j=i; | |
1028 | if (i<0) j+=18; //correct for the negative values | |
1029 | if (!sideA) j+=18; //move to the correct side | |
1b45e564 | 1030 | |
644666df | 1031 | Float_t deltaPhi = 0.; |
1032 | Float_t phiEdge=sectorWidth*i; | |
1033 | if (inphi>phiEdge) {deltaPhi=phiEdge+sectorWidth-inphi;} | |
1034 | else if ((outphi>=phiEdge) && (outphi<(phiEdge+sectorWidth))) {deltaPhi=outphi-phiEdge;} | |
1035 | else {deltaPhi=sectorWidth;} | |
1b45e564 | 1036 | |
644666df | 1037 | Float_t v = fVoltageMap[(layer==1)?(j):(j+36)]; |
1b45e564 | 1038 | if (v<=fBadIROCthreshhold) |
1039 | { | |
1040 | trackLengthInBad+=deltaPhi; | |
f85a3764 | 1041 | lengthFractionInBadSectors=1.; |
1042 | } | |
1b45e564 | 1043 | if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold) |
f85a3764 | 1044 | { |
1045 | trackLengthInLowGain+=deltaPhi; | |
1046 | lengthFractionInBadSectors=1.; | |
1047 | } | |
644666df | 1048 | } |
1049 | ||
1050 | //for now low gain and bad (off) chambers are treated equally | |
f85a3764 | 1051 | if (trackLengthTotal>0) |
1052 | lengthFractionInBadSectors = (trackLengthInLowGain+trackLengthInBad)/trackLengthTotal; | |
644666df | 1053 | |
1054 | //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 | 1055 | |
1056 | if (lengthFractionInBadSectors>fMaxBadLengthFraction) | |
644666df | 1057 | { |
1058 | //printf("%%%%%%%% %s kChamberLowGain\n",(layer==1)?"IROC":"OROC"); | |
1059 | return kChamberLowGain; | |
1060 | } | |
1b45e564 | 1061 | |
644666df | 1062 | return kChamberHighGain; |
1063 | } | |
1064 | ||
1b45e564 | 1065 | |
644666df | 1066 | //_____________________________________________________________________________ |
1067 | Float_t AliTPCPIDResponse::MaxClusterRadius(const AliVTrack* track) const | |
1068 | { | |
1069 | //return the radius of the outermost padrow containing a cluster in TPC | |
1070 | //for the track | |
1071 | const TBits* clusterMap=track->GetTPCClusterMapPtr(); | |
1072 | if (!clusterMap) return 0.; | |
1073 | ||
1074 | //from AliTPCParam, radius of first IROC row | |
1b45e564 | 1075 | const Float_t rfirstIROC = 8.52249984741210938e+01; |
644666df | 1076 | const Float_t padrowHeightIROC = 0.75; |
1077 | const Float_t rfirstOROC0 = 1.35100006103515625e+02; | |
1078 | const Float_t padrowHeightOROC0 = 1.0; | |
1079 | const Float_t rfirstOROC1 = 1.99350006103515625e+02; | |
1080 | const Float_t padrowHeightOROC1 = 1.5; | |
1081 | ||
1082 | Int_t maxPadRow=160; | |
1083 | while ((--maxPadRow)>0 && !clusterMap->TestBitNumber(maxPadRow)){} | |
1084 | if (maxPadRow>126) return rfirstOROC1+(maxPadRow-126-1)*padrowHeightOROC1; | |
1085 | if (maxPadRow>62) return rfirstOROC0+(maxPadRow-62-1)*padrowHeightOROC0; | |
1086 | if (maxPadRow>0) return rfirstIROC+(maxPadRow-1)*padrowHeightIROC; | |
1087 | return 0.0; | |
10d100d4 | 1088 | } |
644666df | 1089 | |
1b45e564 | 1090 | |
644666df | 1091 | //_____________________________________________________________________________ |
1092 | Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const | |
1093 | { | |
1094 | //calculate the coordinates of the apex of the track | |
1095 | Double_t x[3]; | |
1096 | track->GetXYZ(x); | |
1097 | Double_t p[3]; | |
1098 | track->GetPxPyPz(p); | |
1099 | Double_t r = 1./track->OneOverPt()/0.0299792458/magField; //signed - will determine the direction of b | |
1b45e564 | 1100 | //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 | 1101 | //find orthogonal vector (Gram-Schmidt) |
1102 | Double_t alpha = (p[0]*x[0] + p[1]*x[1])/(p[0]*p[0] + p[1]*p[1]); | |
1103 | Double_t b[2]; | |
1b45e564 | 1104 | b[0]=x[0]-alpha*p[0]; |
1105 | b[1]=x[1]-alpha*p[1]; | |
1106 | ||
644666df | 1107 | Double_t norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]); |
1108 | if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE; | |
1b45e564 | 1109 | b[0]/=norm; |
1110 | b[1]/=norm; | |
1111 | b[0]*=r; | |
1112 | b[1]*=r; | |
1113 | b[0]+=x[0]; | |
644666df | 1114 | b[1]+=x[1]; |
1115 | //printf("center: x:%.2f, y:%.2f\n",b[0],b[1]); | |
1b45e564 | 1116 | |
644666df | 1117 | norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]); |
1118 | if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE; | |
1119 | ||
1120 | position[0]=b[0]+b[0]*TMath::Abs(r)/norm; | |
1121 | position[1]=b[1]+b[1]*TMath::Abs(r)/norm; | |
1122 | position[2]=0.; | |
1123 | return kTRUE; | |
f85a3764 | 1124 | } |