]>
Commit | Line | Data |
---|---|---|
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 | |
23 | // ...and some modifications by | |
24 | // Mikolaj Krzewicki, GSI, mikolaj.krzewicki@cern.ch | |
25 | //----------------------------------------------------------------- | |
26 | ||
27 | #include <TGraph.h> | |
28 | #include <TObjArray.h> | |
29 | #include <TSpline.h> | |
30 | #include <TBits.h> | |
31 | #include <TMath.h> | |
32 | ||
33 | #include <AliLog.h> | |
34 | #include "AliExternalTrackParam.h" | |
35 | #include "AliVTrack.h" | |
36 | #include "AliTPCPIDResponse.h" | |
37 | #include "AliTPCdEdxInfo.h" | |
38 | ||
39 | ClassImp(AliTPCPIDResponse); | |
40 | ||
41 | const char* AliTPCPIDResponse::fgkGainScenarioName[fgkNumberOfGainScenarios+1]= | |
42 | { | |
43 | "", //default - no name | |
44 | "1", //all high | |
45 | "2", //OROC only | |
46 | "unknown" | |
47 | }; | |
48 | ||
49 | //_________________________________________________________________________ | |
50 | AliTPCPIDResponse::AliTPCPIDResponse(): | |
51 | TNamed(), | |
52 | fMIP(50.), | |
53 | fRes0(), | |
54 | fResN2(), | |
55 | fKp1(0.0283086), | |
56 | fKp2(2.63394e+01), | |
57 | fKp3(5.04114e-11), | |
58 | fKp4(2.12543), | |
59 | fKp5(4.88663), | |
60 | fUseDatabase(kFALSE), | |
61 | fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios), | |
62 | fVoltageMap(72), | |
63 | fLowGainIROCthreshold(-40), | |
64 | fBadIROCthreshhold(-70), | |
65 | fLowGainOROCthreshold(-40), | |
66 | fBadOROCthreshhold(-40), | |
67 | fMaxBadLengthFraction(0.5), | |
68 | fCurrentResponseFunction(NULL), | |
69 | fCurrentdEdx(0.0), | |
70 | fCurrentNPoints(0), | |
71 | fCurrentGainScenario(kGainScenarioInvalid), | |
72 | fMagField(0.) | |
73 | { | |
74 | // | |
75 | // The default constructor | |
76 | // | |
77 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=0.07;fResN2[i]=0.0;} | |
78 | } | |
79 | ||
80 | //_________________________________________________________________________ | |
81 | AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param): | |
82 | TNamed(), | |
83 | fMIP(param[0]), | |
84 | fRes0(), | |
85 | fResN2(), | |
86 | fKp1(0.0283086), | |
87 | fKp2(2.63394e+01), | |
88 | fKp3(5.04114e-11), | |
89 | fKp4(2.12543), | |
90 | fKp5(4.88663), | |
91 | fUseDatabase(kFALSE), | |
92 | fResponseFunctions(fgkNumberOfParticleSpecies*fgkNumberOfGainScenarios), | |
93 | fVoltageMap(72), | |
94 | fLowGainIROCthreshold(-40), | |
95 | fBadIROCthreshhold(-70), | |
96 | fLowGainOROCthreshold(-40), | |
97 | fBadOROCthreshhold(-40), | |
98 | fMaxBadLengthFraction(0.5), | |
99 | fCurrentResponseFunction(NULL), | |
100 | fCurrentdEdx(0.0), | |
101 | fCurrentNPoints(0), | |
102 | fCurrentGainScenario(kGainScenarioInvalid), | |
103 | fMagField(0.) | |
104 | { | |
105 | // | |
106 | // The main constructor | |
107 | // | |
108 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=param[1];fResN2[i]=param[2];} | |
109 | } | |
110 | ||
111 | //_________________________________________________________________________ | |
112 | AliTPCPIDResponse::AliTPCPIDResponse(const AliTPCPIDResponse& that): | |
113 | TNamed(that), | |
114 | fMIP(that.fMIP), | |
115 | fRes0(), | |
116 | fResN2(), | |
117 | fKp1(that.fKp1), | |
118 | fKp2(that.fKp2), | |
119 | fKp3(that.fKp3), | |
120 | fKp4(that.fKp4), | |
121 | fKp5(that.fKp5), | |
122 | fUseDatabase(that.fUseDatabase), | |
123 | fResponseFunctions(that.fResponseFunctions), | |
124 | fVoltageMap(that.fVoltageMap), | |
125 | fLowGainIROCthreshold(that.fLowGainIROCthreshold), | |
126 | fBadIROCthreshhold(that.fBadIROCthreshhold), | |
127 | fLowGainOROCthreshold(that.fLowGainOROCthreshold), | |
128 | fBadOROCthreshhold(that.fBadOROCthreshhold), | |
129 | fMaxBadLengthFraction(that.fMaxBadLengthFraction), | |
130 | fCurrentResponseFunction(NULL), | |
131 | fCurrentdEdx(0.0), | |
132 | fCurrentNPoints(0), | |
133 | fCurrentGainScenario(kGainScenarioInvalid), | |
134 | fMagField(that.fMagField) | |
135 | { | |
136 | //copy ctor | |
137 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];} | |
138 | } | |
139 | ||
140 | //_________________________________________________________________________ | |
141 | AliTPCPIDResponse& AliTPCPIDResponse::operator=(const AliTPCPIDResponse& that) | |
142 | { | |
143 | //assignment | |
144 | if (&that==this) return *this; | |
145 | TNamed::operator=(that); | |
146 | fMIP=that.fMIP; | |
147 | fKp1=that.fKp1; | |
148 | fKp2=that.fKp2; | |
149 | fKp3=that.fKp3; | |
150 | fKp4=that.fKp4; | |
151 | fKp5=that.fKp5; | |
152 | fUseDatabase=that.fUseDatabase; | |
153 | fResponseFunctions=that.fResponseFunctions; | |
154 | fVoltageMap=that.fVoltageMap; | |
155 | fLowGainIROCthreshold=that.fLowGainIROCthreshold; | |
156 | fBadIROCthreshhold=that.fBadIROCthreshhold; | |
157 | fLowGainOROCthreshold=that.fLowGainOROCthreshold; | |
158 | fBadOROCthreshhold=that.fBadOROCthreshhold; | |
159 | fMaxBadLengthFraction=that.fMaxBadLengthFraction; | |
160 | fMagField=that.fMagField; | |
161 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=that.fRes0[i];fResN2[i]=that.fResN2[i];} | |
162 | return *this; | |
163 | } | |
164 | ||
165 | //_________________________________________________________________________ | |
166 | Double_t AliTPCPIDResponse::Bethe(Double_t betaGamma) const { | |
167 | // | |
168 | // This is the Bethe-Bloch function normalised to 1 at the minimum | |
169 | // WARNING | |
170 | // Simulated and reconstructed Bethe-Bloch differs | |
171 | // Simulated curve is the dNprim/dx | |
172 | // Reconstructed is proportianal dNtot/dx | |
173 | // Temporary fix for production - Simple linear correction function | |
174 | // Future 2 Bethe Bloch formulas needed | |
175 | // 1. for simulation | |
176 | // 2. for reconstructed PID | |
177 | // | |
178 | ||
179 | // const Float_t kmeanCorrection =0.1; | |
180 | Double_t bb= | |
181 | AliExternalTrackParam::BetheBlochAleph(betaGamma,fKp1,fKp2,fKp3,fKp4,fKp5); | |
182 | return bb*fMIP; | |
183 | } | |
184 | ||
185 | //_________________________________________________________________________ | |
186 | void AliTPCPIDResponse::SetBetheBlochParameters(Double_t kp1, | |
187 | Double_t kp2, | |
188 | Double_t kp3, | |
189 | Double_t kp4, | |
190 | Double_t kp5) { | |
191 | // | |
192 | // Set the parameters of the ALEPH Bethe-Bloch formula | |
193 | // | |
194 | fKp1=kp1; | |
195 | fKp2=kp2; | |
196 | fKp3=kp3; | |
197 | fKp4=kp4; | |
198 | fKp5=kp5; | |
199 | } | |
200 | ||
201 | //_________________________________________________________________________ | |
202 | void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2) { | |
203 | // | |
204 | // Set the relative resolution sigma_rel = res0 * sqrt(1+resN2/npoint) | |
205 | // | |
206 | for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=res0;fResN2[i]=resN2;} | |
207 | } | |
208 | ||
209 | //_________________________________________________________________________ | |
210 | Double_t AliTPCPIDResponse::GetExpectedSignal(const Float_t mom, | |
211 | AliPID::EParticleType n) const { | |
212 | // | |
213 | // Calculates the expected PID signal as the function of | |
214 | // the information stored in the track, for the specified particle type | |
215 | // | |
216 | // At the moment, these signals are just the results of calling the | |
217 | // Bethe-Bloch formula. | |
218 | // This can be improved. By taking into account the number of | |
219 | // assigned clusters and/or the track dip angle, for example. | |
220 | // | |
221 | ||
222 | Double_t mass=AliPID::ParticleMassZ(n); | |
223 | if (!fUseDatabase) return Bethe(mom/mass); | |
224 | // | |
225 | const TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n); | |
226 | ||
227 | if (!responseFunction) return Bethe(mom/mass); | |
228 | //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...) | |
229 | // !!! Splines for light nuclei need to be normalised to this factor !!! | |
230 | const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(n),2.3); | |
231 | return fMIP*responseFunction->Eval(mom/mass)*chargeFactor; | |
232 | ||
233 | } | |
234 | ||
235 | //_________________________________________________________________________ | |
236 | Double_t AliTPCPIDResponse::GetExpectedSigma(const Float_t mom, | |
237 | const Int_t nPoints, | |
238 | AliPID::EParticleType n) const { | |
239 | // | |
240 | // Calculates the expected sigma of the PID signal as the function of | |
241 | // the information stored in the track, for the specified particle type | |
242 | // | |
243 | ||
244 | if (nPoints != 0) | |
245 | return GetExpectedSignal(mom,n)*fRes0[0]*sqrt(1. + fResN2[0]/nPoints); | |
246 | else | |
247 | return GetExpectedSignal(mom,n)*fRes0[0]; | |
248 | } | |
249 | ||
250 | ////////////////////////////////////////////////////NEW////////////////////////////// | |
251 | ||
252 | //_________________________________________________________________________ | |
253 | void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario) { | |
254 | // | |
255 | // Set the relative resolution sigma_rel = res0 * sqrt(1+resN2/npoint) | |
256 | // | |
257 | if ((Int_t)gainScenario>(Int_t)fgkNumberOfGainScenarios) return; //TODO: better handling! | |
258 | fRes0[gainScenario]=res0; | |
259 | fResN2[gainScenario]=resN2; | |
260 | } | |
261 | ||
262 | //_________________________________________________________________________ | |
263 | Double_t AliTPCPIDResponse::GetExpectedSignal(Double_t mom, | |
264 | AliPID::EParticleType species, | |
265 | const TSpline3* responseFunction) const | |
266 | { | |
267 | // Calculates the expected PID signal as the function of | |
268 | // the information stored in the track, for the specified particle type | |
269 | // | |
270 | // At the moment, these signals are just the results of calling the | |
271 | // Bethe-Bloch formula. | |
272 | // This can be improved. By taking into account the number of | |
273 | // assigned clusters and/or the track dip angle, for example. | |
274 | // | |
275 | ||
276 | ||
277 | Double_t mass=AliPID::ParticleMass(species); | |
278 | ||
279 | if (!fUseDatabase) return Bethe(mom/mass); | |
280 | ||
281 | if (!responseFunction) return Bethe(mom/mass); | |
282 | return fMIP*responseFunction->Eval(mom/mass); | |
283 | } | |
284 | ||
285 | //_________________________________________________________________________ | |
286 | Double_t AliTPCPIDResponse::GetExpectedSignal(const AliVTrack* track, | |
287 | AliPID::EParticleType species, | |
288 | ETPCdEdxSource dedxSource) | |
289 | { | |
290 | // Calculates the expected PID signal as the function of | |
291 | // the information stored in the track, for the specified particle type | |
292 | // | |
293 | // At the moment, these signals are just the results of calling the | |
294 | // Bethe-Bloch formula. | |
295 | // This can be improved. By taking into account the number of | |
296 | // assigned clusters and/or the track dip angle, for example. | |
297 | // | |
298 | ||
299 | Double_t mom = track->GetTPCmomentum(); | |
300 | Double_t mass=AliPID::ParticleMass(species); | |
301 | ||
302 | if (!fUseDatabase) return Bethe(mom/mass); | |
303 | ||
304 | const TSpline3* responseFunction = GetResponseFunction(track,species,dedxSource); | |
305 | if (!responseFunction) return Bethe(mom/mass); | |
306 | return fMIP*responseFunction->Eval(mom/mass); | |
307 | ||
308 | } | |
309 | ||
310 | //_________________________________________________________________________ | |
311 | TSpline3* AliTPCPIDResponse::GetResponseFunction( AliPID::EParticleType type, | |
312 | AliTPCPIDResponse::ETPCgainScenario gainScenario ) const | |
313 | { | |
314 | //get response function | |
315 | return dynamic_cast<TSpline3*>(fResponseFunctions.At(ResponseFunctionIndex(type,gainScenario))); | |
316 | } | |
317 | ||
318 | //_________________________________________________________________________ | |
319 | TSpline3* AliTPCPIDResponse::GetResponseFunction( const AliVTrack* track, | |
320 | AliPID::EParticleType species, | |
321 | ETPCdEdxSource dedxSource ) | |
322 | { | |
323 | //the splines are stored in an array, different scenarios | |
324 | ||
325 | if (ResponseFunctiondEdxN(track, | |
326 | species, | |
327 | dedxSource)) | |
328 | return fCurrentResponseFunction; | |
329 | ||
330 | return NULL; | |
331 | } | |
332 | ||
333 | //_________________________________________________________________________ | |
334 | void AliTPCPIDResponse::ResetSplines() | |
335 | { | |
336 | //reset the array with splines | |
337 | for (Int_t i=0;i<fResponseFunctions.GetEntriesFast();i++) | |
338 | { | |
339 | fResponseFunctions.AddAt(NULL,i); | |
340 | } | |
341 | } | |
342 | //_________________________________________________________________________ | |
343 | Int_t AliTPCPIDResponse::ResponseFunctionIndex( AliPID::EParticleType species, | |
344 | ETPCgainScenario gainScenario ) const | |
345 | { | |
346 | //get the index in fResponseFunctions given type and scenario | |
347 | return Int_t(species)+Int_t(gainScenario)*fgkNumberOfParticleSpecies; | |
348 | } | |
349 | ||
350 | //_________________________________________________________________________ | |
351 | void AliTPCPIDResponse::SetResponseFunction( TObject* o, | |
352 | AliPID::EParticleType species, | |
353 | ETPCgainScenario gainScenario ) | |
354 | { | |
355 | fResponseFunctions.AddAtAndExpand(o,ResponseFunctionIndex(species,gainScenario)); | |
356 | } | |
357 | ||
358 | //_________________________________________________________________________ | |
359 | Double_t AliTPCPIDResponse::GetExpectedSigma( Double_t mom, | |
360 | Int_t nPoints, | |
361 | AliPID::EParticleType species, | |
362 | ETPCgainScenario gainScenario, | |
363 | const TSpline3* responseFunction) const | |
364 | { | |
365 | // Calculates the expected sigma of the PID signal as the function of | |
366 | // the information stored in the track, for the specified particle type | |
367 | // and dedx scenario | |
368 | // | |
369 | ||
370 | if (!responseFunction) return 999; | |
371 | if (nPoints != 0) | |
372 | return GetExpectedSignal(mom,species,responseFunction) * | |
373 | fRes0[gainScenario] * sqrt(1. + fResN2[gainScenario]/nPoints); | |
374 | else | |
375 | return GetExpectedSignal(mom,species,responseFunction)*fRes0[gainScenario]; | |
376 | } | |
377 | ||
378 | //_________________________________________________________________________ | |
379 | Double_t AliTPCPIDResponse::GetExpectedSigma(const AliVTrack* track, | |
380 | AliPID::EParticleType species, | |
381 | ETPCdEdxSource dedxSource) | |
382 | { | |
383 | // Calculates the expected sigma of the PID signal as the function of | |
384 | // the information stored in the track, for the specified particle type | |
385 | // and dedx scenario | |
386 | // | |
387 | ||
388 | if (!ResponseFunctiondEdxN(track, | |
389 | species, | |
390 | dedxSource)) return 998; //TODO: better handling! | |
391 | ||
392 | return GetExpectedSigma( track->GetTPCmomentum(), | |
393 | fCurrentNPoints, | |
394 | species, | |
395 | fCurrentGainScenario, | |
396 | fCurrentResponseFunction ); | |
397 | } | |
398 | ||
399 | //_________________________________________________________________________ | |
400 | Float_t AliTPCPIDResponse::GetNumberOfSigmas(const AliVTrack* track, | |
401 | AliPID::EParticleType species, | |
402 | ETPCdEdxSource dedxSource) | |
403 | { | |
404 | //calculates the number of sigmas of the PID signal from the expected value | |
405 | //for a gicven particle species in the presence of multiple gain scenarios | |
406 | //inside the TPC | |
407 | ||
408 | if (!ResponseFunctiondEdxN(track, | |
409 | species, | |
410 | dedxSource)) return 997; //TODO: better handling! | |
411 | ||
412 | Double_t mom = track->GetTPCmomentum(); | |
413 | Double_t bethe=GetExpectedSignal(mom,species,fCurrentResponseFunction); | |
414 | Double_t sigma=GetExpectedSigma( mom, | |
415 | fCurrentNPoints, | |
416 | species, | |
417 | fCurrentGainScenario, | |
418 | fCurrentResponseFunction ); | |
419 | return (fCurrentdEdx-bethe)/sigma; | |
420 | } | |
421 | ||
422 | //_________________________________________________________________________ | |
423 | Bool_t AliTPCPIDResponse::ResponseFunctiondEdxN( const AliVTrack* track, | |
424 | AliPID::EParticleType species, | |
425 | ETPCdEdxSource dedxSource ) | |
426 | { | |
427 | // Calculates the right parameters for PID | |
428 | // dEdx parametrization for the proper gain scenario, dEdx | |
429 | // and NPoints used for dEdx | |
430 | // based on the track geometry (which chambers it crosses) for the specified particle type | |
431 | // and preferred source of dedx. | |
432 | // returns true on success | |
433 | ||
434 | Double32_t signal[4]; //0: IROC, 1: OROC medium, 2:OROC long, 3: OROC all (def. truncation used) | |
435 | Char_t ncl[3]; //same | |
436 | Char_t nrows[3]; //same | |
437 | AliTPCdEdxInfo* dEdxInfo = track->GetTPCdEdxInfo(); | |
438 | ||
439 | if (!dEdxInfo && dedxSource!=kdEdxDefault) //in one case its ok if we dont have the info | |
440 | { | |
441 | AliError("AliTPCdEdxInfo not available"); | |
442 | InvalidateCurrentValues(); | |
443 | return kFALSE; | |
444 | } | |
445 | ||
446 | if (dEdxInfo) dEdxInfo->GetTPCSignalRegionInfo(signal,ncl,nrows); | |
447 | ||
448 | //check if we cross a bad OROC in which case we reject | |
449 | EChamberStatus trackStatus = TrackStatus(track,2); | |
450 | if (trackStatus==kChamberOff || trackStatus==kChamberLowGain) | |
451 | { | |
452 | InvalidateCurrentValues(); | |
453 | return kFALSE; | |
454 | } | |
455 | ||
456 | switch (dedxSource) | |
457 | { | |
458 | case kdEdxDefault: | |
459 | { | |
460 | fCurrentdEdx = track->GetTPCsignal(); | |
461 | fCurrentNPoints = track->GetTPCsignalN(); | |
462 | fCurrentGainScenario = kDefault; | |
463 | break; | |
464 | } | |
465 | case kdEdxOROC: | |
466 | { | |
467 | fCurrentdEdx = signal[3]; | |
468 | fCurrentNPoints = ncl[2]+ncl[1]; | |
469 | fCurrentGainScenario = kOROChigh; | |
470 | break; | |
471 | } | |
472 | case kdEdxHybrid: | |
473 | { | |
474 | //if we cross a bad IROC we use OROC dedx, if we dont we use combined | |
475 | EChamberStatus status = TrackStatus(track,1); | |
476 | if (status!=kChamberHighGain) | |
477 | { | |
478 | fCurrentdEdx = signal[3]; | |
479 | fCurrentNPoints = ncl[2]+ncl[1]; | |
480 | fCurrentGainScenario = kOROChigh; | |
481 | } | |
482 | else | |
483 | { | |
484 | fCurrentdEdx = track->GetTPCsignal(); | |
485 | fCurrentNPoints = track->GetTPCsignalN(); | |
486 | fCurrentGainScenario = kALLhigh; | |
487 | } | |
488 | break; | |
489 | } | |
490 | default: | |
491 | { | |
492 | fCurrentdEdx = 0.; | |
493 | fCurrentNPoints = 0; | |
494 | fCurrentGainScenario = kGainScenarioInvalid; | |
495 | return kFALSE; | |
496 | } | |
497 | } | |
498 | TObject* obj = fResponseFunctions.UncheckedAt(ResponseFunctionIndex(species,fCurrentGainScenario)); | |
499 | fCurrentResponseFunction = dynamic_cast<TSpline3*>(obj); //TODO:maybe static cast? | |
500 | return kTRUE; | |
501 | } | |
502 | ||
503 | //_________________________________________________________________________ | |
504 | Bool_t AliTPCPIDResponse::sectorNumbersInOut(const AliVTrack* track, | |
505 | Double_t innerRadius, | |
506 | Double_t outerRadius, | |
507 | Float_t& inphi, | |
508 | Float_t& outphi, | |
509 | Int_t& in, | |
510 | Int_t& out ) const | |
511 | { | |
512 | //calculate the sector numbers (equivalent to IROC chamber numbers) a track crosses | |
513 | //for OROC chamber numbers add 36 | |
514 | //returned angles are between (0,2pi) | |
515 | ||
516 | Double_t trackPositionInner[3]; | |
517 | Double_t trackPositionOuter[3]; | |
518 | ||
519 | Bool_t trackAtInner=kTRUE; | |
520 | Bool_t trackAtOuter=kTRUE; | |
521 | const AliExternalTrackParam* ip = track->GetInnerParam(); | |
522 | ||
523 | //if there is no inner param this could mean we're using the AOD track, | |
524 | //we still can extrapolate from the vertex - so use those params. | |
525 | if (ip) track=ip; | |
526 | ||
527 | if (!track->GetXYZAt(innerRadius, fMagField, trackPositionInner)) | |
528 | trackAtInner=kFALSE; | |
529 | ||
530 | if (!track->GetXYZAt(outerRadius, fMagField, trackPositionOuter)) | |
531 | trackAtOuter=kFALSE; | |
532 | ||
533 | if (!trackAtInner) | |
534 | { | |
535 | //if we dont even enter inner radius we do nothing | |
536 | inphi=0.0; | |
537 | outphi=0.0; | |
538 | in=0; | |
539 | out=0; | |
540 | return kFALSE; | |
541 | } | |
542 | ||
543 | if (!trackAtOuter) | |
544 | { | |
545 | //if we don't reach the outer radius check that the apex is indeed within the outer radius and use apex position | |
546 | Bool_t haveApex = TrackApex(track, fMagField, trackPositionOuter); | |
547 | Float_t apexRadius = TMath::Sqrt(trackPositionOuter[0]*trackPositionOuter[0]+trackPositionOuter[1]*trackPositionOuter[1]); | |
548 | if ( haveApex && apexRadius<=outerRadius && apexRadius>innerRadius) | |
549 | { | |
550 | //printf("pt: %.2f, apexRadius: %.2f(%s), x: %.2f, y: %.2f\n",track->Pt(),apexRadius,(haveApex)?"OK":"BAD",trackPositionOuter[0],trackPositionOuter[1]); | |
551 | } | |
552 | else | |
553 | { | |
554 | inphi=0.0; | |
555 | outphi=0.0; | |
556 | in=0; | |
557 | out=0; | |
558 | return kFALSE; | |
559 | } | |
560 | } | |
561 | ||
562 | inphi = TMath::ATan2(trackPositionInner[1],trackPositionInner[0]); | |
563 | outphi = TMath::ATan2(trackPositionOuter[1], trackPositionOuter[0]); | |
564 | ||
565 | if (inphi<0) {inphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi | |
566 | if (outphi<0) {outphi+=TMath::TwoPi();} //because ATan2 gives angles between -pi,pi | |
567 | ||
568 | in = sectorNumber(inphi); | |
569 | out = sectorNumber(outphi); | |
570 | ||
571 | //for the C side (positive z) offset by 18 | |
572 | if (trackPositionInner[2]>0.0) | |
573 | { | |
574 | in+=18; | |
575 | out+=18; | |
576 | } | |
577 | return kTRUE; | |
578 | } | |
579 | ||
580 | //_____________________________________________________________________________ | |
581 | Int_t AliTPCPIDResponse::sectorNumber(Double_t phi) const | |
582 | { | |
583 | //calculate sector number | |
584 | const Float_t width=TMath::TwoPi()/18.; | |
585 | return TMath::Floor(phi/width); | |
586 | } | |
587 | ||
588 | //_____________________________________________________________________________ | |
589 | void AliTPCPIDResponse::Print(Option_t* /*option*/) const | |
590 | { | |
591 | //Print info | |
592 | fResponseFunctions.Print(); | |
593 | } | |
594 | ||
595 | //_____________________________________________________________________________ | |
596 | AliTPCPIDResponse::EChamberStatus AliTPCPIDResponse::TrackStatus(const AliVTrack* track, Int_t layer) const | |
597 | { | |
598 | //status of the track: if it crosses any bad chambers return kChamberOff; | |
599 | //IROC:layer=1, OROC:layer=2 | |
600 | if (layer<1 || layer>2) layer=1; | |
601 | Int_t in=0; | |
602 | Int_t out=0; | |
603 | Float_t inphi=0.; | |
604 | Float_t outphi=0.; | |
605 | Float_t innerRadius = (layer==1)?83.0:133.7; | |
606 | Float_t outerRadius = (layer==1)?133.5:247.7; | |
607 | if (!sectorNumbersInOut(track, | |
608 | innerRadius, | |
609 | outerRadius, | |
610 | inphi, | |
611 | outphi, | |
612 | in, | |
613 | out)) return kChamberOff; | |
614 | ||
615 | Bool_t sideA = kTRUE; | |
616 | ||
617 | if (((in/18)%2==1) && ((out/18)%2==1)) sideA=kFALSE; | |
618 | ||
619 | in=in%18; | |
620 | out=out%18; | |
621 | ||
622 | if (TMath::Abs(in-out)>9) | |
623 | { | |
624 | if (TMath::Max(in,out)==out) | |
625 | { | |
626 | Int_t tmp=out; | |
627 | out = in; | |
628 | in = tmp; | |
629 | Float_t tmpphi=outphi; | |
630 | outphi=inphi; | |
631 | inphi=tmpphi; | |
632 | } | |
633 | in-=18; | |
634 | inphi-=TMath::TwoPi(); | |
635 | } | |
636 | else | |
637 | { | |
638 | if (TMath::Max(in,out)==in) | |
639 | { | |
640 | Int_t tmp=out; | |
641 | out=in; | |
642 | in=tmp; | |
643 | Float_t tmpphi=outphi; | |
644 | outphi=inphi; | |
645 | inphi=tmpphi; | |
646 | } | |
647 | } | |
648 | ||
649 | Float_t trackLengthInBad = 0.; | |
650 | Float_t trackLengthInLowGain = 0.; | |
651 | Float_t trackLengthTotal = TMath::Abs(outphi-inphi); | |
652 | ||
653 | const Float_t sectorWidth = TMath::TwoPi()/18.; | |
654 | ||
655 | for (Int_t i=in; i<=out; i++) | |
656 | { | |
657 | int j=i; | |
658 | if (i<0) j+=18; //correct for the negative values | |
659 | if (!sideA) j+=18; //move to the correct side | |
660 | ||
661 | Float_t deltaPhi = 0.; | |
662 | Float_t phiEdge=sectorWidth*i; | |
663 | if (inphi>phiEdge) {deltaPhi=phiEdge+sectorWidth-inphi;} | |
664 | else if ((outphi>=phiEdge) && (outphi<(phiEdge+sectorWidth))) {deltaPhi=outphi-phiEdge;} | |
665 | else {deltaPhi=sectorWidth;} | |
666 | ||
667 | Float_t v = fVoltageMap[(layer==1)?(j):(j+36)]; | |
668 | if (v<=fBadIROCthreshhold) trackLengthInBad+=deltaPhi; | |
669 | if (v<=fLowGainIROCthreshold && v>fBadIROCthreshhold) trackLengthInLowGain+=deltaPhi; | |
670 | } | |
671 | ||
672 | //for now low gain and bad (off) chambers are treated equally | |
673 | Float_t lengthFractionInBadSectors = (trackLengthInLowGain+trackLengthInBad)/trackLengthTotal; | |
674 | ||
675 | //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); | |
676 | ||
677 | if (lengthFractionInBadSectors>fMaxBadLengthFraction) | |
678 | { | |
679 | //printf("%%%%%%%% %s kChamberLowGain\n",(layer==1)?"IROC":"OROC"); | |
680 | return kChamberLowGain; | |
681 | } | |
682 | ||
683 | return kChamberHighGain; | |
684 | } | |
685 | ||
686 | //_____________________________________________________________________________ | |
687 | Float_t AliTPCPIDResponse::MaxClusterRadius(const AliVTrack* track) const | |
688 | { | |
689 | //return the radius of the outermost padrow containing a cluster in TPC | |
690 | //for the track | |
691 | const TBits* clusterMap=track->GetTPCClusterMapPtr(); | |
692 | if (!clusterMap) return 0.; | |
693 | ||
694 | //from AliTPCParam, radius of first IROC row | |
695 | const Float_t rfirstIROC = 8.52249984741210938e+01; | |
696 | const Float_t padrowHeightIROC = 0.75; | |
697 | const Float_t rfirstOROC0 = 1.35100006103515625e+02; | |
698 | const Float_t padrowHeightOROC0 = 1.0; | |
699 | const Float_t rfirstOROC1 = 1.99350006103515625e+02; | |
700 | const Float_t padrowHeightOROC1 = 1.5; | |
701 | ||
702 | Int_t maxPadRow=160; | |
703 | while ((--maxPadRow)>0 && !clusterMap->TestBitNumber(maxPadRow)){} | |
704 | if (maxPadRow>126) return rfirstOROC1+(maxPadRow-126-1)*padrowHeightOROC1; | |
705 | if (maxPadRow>62) return rfirstOROC0+(maxPadRow-62-1)*padrowHeightOROC0; | |
706 | if (maxPadRow>0) return rfirstIROC+(maxPadRow-1)*padrowHeightIROC; | |
707 | return 0.0; | |
708 | } | |
709 | ||
710 | //_____________________________________________________________________________ | |
711 | void AliTPCPIDResponse::InvalidateCurrentValues() | |
712 | { | |
713 | //make the "current" stored values from the last processed track invalid | |
714 | fCurrentResponseFunction=NULL; | |
715 | fCurrentdEdx=0.; | |
716 | fCurrentNPoints=0; | |
717 | fCurrentGainScenario=kGainScenarioInvalid; | |
718 | } | |
719 | ||
720 | //_____________________________________________________________________________ | |
721 | Bool_t AliTPCPIDResponse::TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const | |
722 | { | |
723 | //calculate the coordinates of the apex of the track | |
724 | Double_t x[3]; | |
725 | track->GetXYZ(x); | |
726 | Double_t p[3]; | |
727 | track->GetPxPyPz(p); | |
728 | Double_t r = 1./track->OneOverPt()/0.0299792458/magField; //signed - will determine the direction of b | |
729 | //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); | |
730 | //find orthogonal vector (Gram-Schmidt) | |
731 | Double_t alpha = (p[0]*x[0] + p[1]*x[1])/(p[0]*p[0] + p[1]*p[1]); | |
732 | Double_t b[2]; | |
733 | b[0]=x[0]-alpha*p[0]; | |
734 | b[1]=x[1]-alpha*p[1]; | |
735 | ||
736 | Double_t norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]); | |
737 | if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE; | |
738 | b[0]/=norm; | |
739 | b[1]/=norm; | |
740 | b[0]*=r; | |
741 | b[1]*=r; | |
742 | b[0]+=x[0]; | |
743 | b[1]+=x[1]; | |
744 | //printf("center: x:%.2f, y:%.2f\n",b[0],b[1]); | |
745 | ||
746 | norm = TMath::Sqrt(b[0]*b[0]+b[1]*b[1]); | |
747 | if (TMath::AreEqualAbs(norm,0.0,1e-10)) return kFALSE; | |
748 | ||
749 | position[0]=b[0]+b[0]*TMath::Abs(r)/norm; | |
750 | position[1]=b[1]+b[1]*TMath::Abs(r)/norm; | |
751 | position[2]=0.; | |
752 | return kTRUE; | |
753 | } | |
754 |