- introduction of gain scenarios (e.g. OROC only - for homogeneous gain in 11h)
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliTPCPIDResponse.cxx
CommitLineData
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
10d100d4 25//-----------------------------------------------------------------
26
d2aa6df0 27#include <TGraph.h>
28#include <TObjArray.h>
29#include <TSpline.h>
644666df 30#include <TBits.h>
00a38d07 31#include <TMath.h>
d2aa6df0 32
644666df 33#include <AliLog.h>
10d100d4 34#include "AliExternalTrackParam.h"
644666df 35#include "AliVTrack.h"
d2aa6df0 36#include "AliTPCPIDResponse.h"
644666df 37#include "AliTPCdEdxInfo.h"
38
39ClassImp(AliTPCPIDResponse);
d2aa6df0 40
644666df 41const char* AliTPCPIDResponse::fgkGainScenarioName[fgkNumberOfGainScenarios+1]=
42{
43 "", //default - no name
44 "1", //all high
45 "2", //OROC only
46 "unknown"
47};
10d100d4 48
49//_________________________________________________________________________
50AliTPCPIDResponse::AliTPCPIDResponse():
644666df 51 TNamed(),
d2aa6df0 52 fMIP(50.),
644666df 53 fRes0(),
54 fResN2(),
d2aa6df0 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),
644666df 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.)
10d100d4 73{
74 //
75 // The default constructor
76 //
644666df 77 for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=0.07;fResN2[i]=0.0;}
10d100d4 78}
79
80//_________________________________________________________________________
d2aa6df0 81AliTPCPIDResponse::AliTPCPIDResponse(const Double_t *param):
644666df 82 TNamed(),
d2aa6df0 83 fMIP(param[0]),
644666df 84 fRes0(),
85 fResN2(),
d2aa6df0 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),
644666df 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.)
10d100d4 104{
105 //
106 // The main constructor
107 //
644666df 108 for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=param[1];fResN2[i]=param[2];}
109}
110
111//_________________________________________________________________________
112AliTPCPIDResponse::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//_________________________________________________________________________
141AliTPCPIDResponse& 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;
10d100d4 163}
164
d2aa6df0 165//_________________________________________________________________________
10d100d4 166Double_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 //
d2aa6df0 178
179// const Float_t kmeanCorrection =0.1;
10d100d4 180 Double_t bb=
181 AliExternalTrackParam::BetheBlochAleph(betaGamma,fKp1,fKp2,fKp3,fKp4,fKp5);
182 return bb*fMIP;
183}
184
185//_________________________________________________________________________
186void 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}
644666df 200
10d100d4 201//_________________________________________________________________________
202void AliTPCPIDResponse::SetSigma(Float_t res0, Float_t resN2) {
203 //
204 // Set the relative resolution sigma_rel = res0 * sqrt(1+resN2/npoint)
205 //
644666df 206 for (Int_t i=0; i<fgkNumberOfGainScenarios; i++) {fRes0[i]=res0;fResN2[i]=resN2;}
10d100d4 207}
208
209//_________________________________________________________________________
210Double_t AliTPCPIDResponse::GetExpectedSignal(const Float_t mom,
d2aa6df0 211 AliPID::EParticleType n) const {
10d100d4 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
539a5a59 222 Double_t mass=AliPID::ParticleMassZ(n);
644666df 223 if (!fUseDatabase) return Bethe(mom/mass);
d2aa6df0 224 //
644666df 225 const TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n);
00a38d07 226
d2aa6df0 227 if (!responseFunction) return Bethe(mom/mass);
00a38d07 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;
d2aa6df0 232
10d100d4 233}
234
235//_________________________________________________________________________
236Double_t AliTPCPIDResponse::GetExpectedSigma(const Float_t mom,
644666df 237 const Int_t nPoints,
238 AliPID::EParticleType n) const {
10d100d4 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 //
644666df 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//_________________________________________________________________________
253void 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//_________________________________________________________________________
263Double_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.
10d100d4 274 //
275
644666df 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//_________________________________________________________________________
286Double_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//_________________________________________________________________________
311TSpline3* 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//_________________________________________________________________________
319TSpline3* 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//_________________________________________________________________________
334void 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//_________________________________________________________________________
343Int_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//_________________________________________________________________________
351void AliTPCPIDResponse::SetResponseFunction( TObject* o,
352 AliPID::EParticleType species,
353 ETPCgainScenario gainScenario )
354{
355 fResponseFunctions.AddAtAndExpand(o,ResponseFunctionIndex(species,gainScenario));
356}
357
358//_________________________________________________________________________
359Double_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;
10d100d4 371 if (nPoints != 0)
644666df 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//_________________________________________________________________________
379Double_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//_________________________________________________________________________
400Float_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//_________________________________________________________________________
423Bool_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//_________________________________________________________________________
504Bool_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//_____________________________________________________________________________
581Int_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//_____________________________________________________________________________
589void AliTPCPIDResponse::Print(Option_t* /*option*/) const
590{
591 //Print info
592 fResponseFunctions.Print();
593}
594
595//_____________________________________________________________________________
596AliTPCPIDResponse::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 }
10d100d4 636 else
644666df 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//_____________________________________________________________________________
687Float_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;
10d100d4 708}
644666df 709
710//_____________________________________________________________________________
711void 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//_____________________________________________________________________________
721Bool_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