TPC PID update:
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliTPCPIDResponse.h
1 #ifndef ALITPCPIDRESPONSE_H
2 #define ALITPCPIDRESPONSE_H
3 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4  * See cxx source for full Copyright notice                               */
5
6 //-------------------------------------------------------
7 //                    TPC PID class
8 // A very naive design... Should be made better by the detector experts...
9 //   Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch 
10 // With many additions and modifications suggested by
11 //      Alexander Kalweit, GSI, alexander.philipp.kalweit@cern.ch
12 //      Dariusz Miskowiec, GSI, D.Miskowiec@gsi.de
13 // ...and some modifications by
14 //      Mikolaj Krzewicki, GSI, mikolaj.krzewicki@cern.ch
15 // ...and some modifications plus eta correction functions by
16 //      Benjamin Hess, University of Tuebingen, bhess@cern.ch
17 //-------------------------------------------------------
18 #include <Rtypes.h>
19
20 #include <TNamed.h>
21 #include <TVectorF.h>
22 #include <TObjArray.h>
23
24 #include "AliPID.h"
25 #include "AliVTrack.h"
26
27 class TH2D;
28 class TSpline3;
29
30 class AliTPCPIDResponse: public TNamed {
31 public:
32   AliTPCPIDResponse();
33   AliTPCPIDResponse(const Double_t *param);
34   AliTPCPIDResponse(const AliTPCPIDResponse&);
35   AliTPCPIDResponse& operator=(const AliTPCPIDResponse&);
36   virtual ~AliTPCPIDResponse();
37
38   enum EChamberStatus {
39     kChamberOff=0,
40     kChamberHighGain=1,
41     kChamberLowGain=2,
42     kChamberInvalid=3
43   };
44   
45   enum ETPCgainScenario {
46     kDefault= 0,
47     kALLhigh = 1,
48     kOROChigh = 2,
49     kGainScenarioInvalid = 3
50   };
51
52   static const Int_t fgkNumberOfParticleSpecies=AliPID::kSPECIESC;
53   static const Int_t fgkNumberOfGainScenarios=3;
54   static const Int_t fgkNumberOfdEdxSourceScenarios=3;
55
56   enum ETPCdEdxSource {
57     kdEdxDefault=0,        // use combined dEdx from IROC+OROC (assumes ideal detector)
58     kdEdxOROC=1,       // use only OROC
59     kdEdxHybrid=2,   // Use IROC+OROC dEdx only where IROCS are good (high gain), otherwise fall back to OROC only
60     kdEdxInvalid=3     //invalid
61   };
62
63   void SetSigma(Float_t res0, Float_t resN2);
64   void SetBetheBlochParameters(Double_t kp1,
65                                Double_t kp2,
66                                Double_t kp3,
67                                Double_t kp4,
68                                Double_t kp5
69                                );
70   void SetMip(Float_t mip) { fMIP = mip; } // Set overall normalisation; mean dE/dx for MIP
71   Double_t Bethe(Double_t bg) const;
72   void SetUseDatabase(Bool_t useDatabase) { fUseDatabase = useDatabase;}
73   Bool_t GetUseDatabase() const { return fUseDatabase;}
74   
75   void SetResponseFunction(AliPID::EParticleType type, TObject * const o) { fResponseFunctions.AddAt(o,(Int_t)type); }
76   const TObject * GetResponseFunction(AliPID::EParticleType type) { return fResponseFunctions.At((Int_t)type); }
77   void SetVoltage(Int_t n, Float_t v) {fVoltageMap[n]=v;}
78   void SetVoltageMap(const TVectorF& a) {fVoltageMap=a;} //resets ownership, ~ will not delete contents
79   Float_t GetVoltage(Int_t n) const {return fVoltageMap[n];}
80   void SetLowGainIROCthreshold(Float_t v) {fLowGainIROCthreshold=v;}
81   void SetBadIROCthreshold(Float_t v) {fBadIROCthreshhold=v;}
82   void SetLowGainOROCthreshold(Float_t v) {fLowGainOROCthreshold=v;}
83   void SetBadOROCthreshold(Float_t v) {fBadOROCthreshhold=v;}
84   void SetMaxBadLengthFraction(Float_t f) {fMaxBadLengthFraction=f;}
85
86   void SetMagField(Double_t mf) { fMagField=mf; }
87   
88   const TH2D* GetEtaCorrMap() const { return fhEtaCorr; };
89   Bool_t SetEtaCorrMap(TH2D* hMap);
90   
91   Double_t GetEtaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
92   
93   Double_t GetEtaCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
94
95   const TH2D* GetSigmaPar1Map() const { return fhEtaSigmaPar1; };
96   Double_t GetSigmaPar0() const { return fSigmaPar0; };
97   Bool_t SetSigmaParams(TH2D* hSigmaPar1Map, Double_t sigmaPar0);
98   
99   Double_t GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
100
101   //NEW
102   void SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario );
103   Double_t GetExpectedSignal( const AliVTrack* track,
104                               AliPID::EParticleType species,
105                               ETPCdEdxSource dedxSource = kdEdxDefault,
106                               Bool_t correctEta = kFALSE) const;//TODO: In future, default kTRUE
107   Double_t GetExpectedSigma( const AliVTrack* track, 
108                              AliPID::EParticleType species,
109                              ETPCdEdxSource dedxSource = kdEdxDefault,
110                              Bool_t correctEta = kFALSE) const;//TODO: In future, default kTRUE
111   Float_t GetNumberOfSigmas( const AliVTrack* track,
112                              AliPID::EParticleType species,
113                              ETPCdEdxSource dedxSource = kdEdxDefault,
114                              Bool_t correctEta = kFALSE) const;//TODO: In future, default kTRUE
115
116   void SetResponseFunction(TObject* o,
117                            AliPID::EParticleType type,
118                            ETPCgainScenario gainScenario);
119   void Print(Option_t* option="") const;
120   TSpline3* GetResponseFunction( AliPID::EParticleType species,
121                                  ETPCgainScenario gainScenario ) const;
122   TSpline3* GetResponseFunction( const AliVTrack* track,
123                                  AliPID::EParticleType species,
124                                  ETPCdEdxSource dedxSource = kdEdxDefault) const;
125   Bool_t ResponseFunctiondEdxN(const AliVTrack* track, 
126                                AliPID::EParticleType species,
127                                ETPCdEdxSource dedxSource,
128                                Double_t& dEdx, Int_t& nPoints, ETPCgainScenario& gainScenario, TSpline3** responseFunction) const;
129   Bool_t sectorNumbersInOut(const AliVTrack* track, 
130                             Double_t innerRadius, Double_t outerRadius, 
131                             Float_t& phiIn, Float_t& phiOut, 
132                             Int_t& in, Int_t& out ) const;
133   AliTPCPIDResponse::EChamberStatus TrackStatus(const AliVTrack* track, Int_t layer) const;
134   Float_t MaxClusterRadius(const AliVTrack* track) const;
135   Bool_t TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const;
136   static const char* GainScenarioName(Int_t n) {return fgkGainScenarioName[(n>fgkNumberOfdEdxSourceScenarios)?fgkNumberOfdEdxSourceScenarios+1:n];}
137   Int_t ResponseFunctionIndex( AliPID::EParticleType species,
138                                ETPCgainScenario gainScenario ) const;
139   void ResetSplines();
140
141   //OLD
142   Double_t GetExpectedSignal(const Float_t mom,
143                      AliPID::EParticleType n=AliPID::kKaon) const;
144   Double_t GetExpectedSigma(const Float_t mom, const Int_t nPoints,
145                             AliPID::EParticleType n=AliPID::kKaon) const;
146   Float_t  GetNumberOfSigmas(const Float_t mom, 
147                              const Float_t dEdx, 
148                                                const Int_t nPoints,
149                              AliPID::EParticleType n=AliPID::kKaon) const {
150     //
151     // Deprecated function (for backward compatibility). Please use 
152     // GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource )
153     // instead!TODO
154     //
155     
156     Double_t bethe=GetExpectedSignal(mom,n);
157     Double_t sigma=GetExpectedSigma(mom,nPoints,n);
158     return (dEdx-bethe)/sigma;
159   }
160
161   Double_t GetMIP() const { return fMIP;} 
162   Float_t  GetRes0()  const { return fRes0[0];  }
163   Float_t  GetResN2() const { return fResN2[0]; }
164   Float_t  GetRes0(ETPCgainScenario s)  const { return fRes0[s];  }
165   Float_t  GetResN2(ETPCgainScenario s) const { return fResN2[s]; }
166
167 protected:
168   Double_t GetExpectedSignal(const AliVTrack* track,
169                              AliPID::EParticleType species,
170                              Double_t dEdx,
171                              const TSpline3* responseFunction,
172                              Bool_t correctEta) const; 
173   
174   Double_t GetExpectedSigma(const AliVTrack* track, 
175                             AliPID::EParticleType species,
176                             ETPCgainScenario gainScenario,
177                             Double_t dEdx,
178                             Int_t nPoints,
179                             const TSpline3* responseFunction,
180                             Bool_t correctEta) const;
181                              
182   Double_t GetEtaCorrection(const AliVTrack *track, Double_t dEdxSplines) const;
183   
184   Double_t GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species,
185                         Double_t dEdx, const TSpline3* responseFunction) const;
186   
187 private:
188   Float_t fMIP;          // dEdx for MIP
189   Float_t fRes0[fgkNumberOfGainScenarios];  // relative dEdx resolution  rel sigma = fRes0*sqrt(1+fResN2/npoint)
190   Float_t fResN2[fgkNumberOfGainScenarios]; // relative Npoint dependence rel  sigma = fRes0*sqrt(1+fResN2/npoint)
191
192   Double_t fKp1;   // Parameters
193   Double_t fKp2;   //    of
194   Double_t fKp3;   // the ALEPH
195   Double_t fKp4;   // Bethe-Bloch
196   Double_t fKp5;   // formula
197
198   Bool_t fUseDatabase; // flag if fine-tuned database-response or simple ALEPH BB should be used
199   
200   TObjArray fResponseFunctions; //! ObjArray of response functions individually for each particle
201   TVectorF fVoltageMap; //!stores a map of voltages wrt nominal for all chambers
202   Float_t fLowGainIROCthreshold;  //voltage threshold below which the IROC is considered low gain
203   Float_t fBadIROCthreshhold;     //voltage threshold for bad IROCS
204   Float_t fLowGainOROCthreshold;  //voltage threshold below which the OROC is considered low gain
205   Float_t fBadOROCthreshhold;     //voltage threshold for bad OROCS
206   Float_t fMaxBadLengthFraction;  //the maximum allowed fraction of track length in a bad sector.
207
208   Int_t sectorNumber(Double_t phi) const;
209
210   Double_t fMagField;  //! Magnetic field
211
212   static const char* fgkGainScenarioName[fgkNumberOfGainScenarios+1];
213
214   TH2D* fhEtaCorr; //! Map for TPC eta correction
215   TH2D* fhEtaSigmaPar1; //! Map for parameter 1 of the dEdx sigma parametrisation
216   
217   Double_t fSigmaPar0; // Parameter 0 of the dEdx sigma parametrisation
218
219   ClassDef(AliTPCPIDResponse,5)   // TPC PID class
220 };
221
222 #endif
223
224