- introduction of gain scenarios (e.g. OROC only - for homogeneous gain in 11h)
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliTPCPIDResponse.cxx
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