Update master to aliroot
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliTPCPIDResponse.h
CommitLineData
10d100d4 1#ifndef ALITPCPIDRESPONSE_H
2#define ALITPCPIDRESPONSE_H
8c6a71ab 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
aea7a46d 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
f84b18dd 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
8c6a71ab 17//-------------------------------------------------------
18#include <Rtypes.h>
644666df 19
20#include <TNamed.h>
21#include <TVectorF.h>
d2aa6df0 22#include <TObjArray.h>
87da0205 23#include <TF1.h>
8c6a71ab 24
aea7a46d 25#include "AliPID.h"
f84b18dd 26#include "AliVTrack.h"
aea7a46d 27
f84b18dd 28class TH2D;
644666df 29class TSpline3;
30
31class AliTPCPIDResponse: public TNamed {
8c6a71ab 32public:
10d100d4 33 AliTPCPIDResponse();
1b45e564 34 //TODO Remove? AliTPCPIDResponse(const Double_t *param);
644666df 35 AliTPCPIDResponse(const AliTPCPIDResponse&);
36 AliTPCPIDResponse& operator=(const AliTPCPIDResponse&);
f84b18dd 37 virtual ~AliTPCPIDResponse();
644666df 38
39 enum EChamberStatus {
40 kChamberOff=0,
41 kChamberHighGain=1,
42 kChamberLowGain=2,
43 kChamberInvalid=3
44 };
45
46 enum ETPCgainScenario {
47 kDefault= 0,
48 kALLhigh = 1,
49 kOROChigh = 2,
50 kGainScenarioInvalid = 3
51 };
52
53 static const Int_t fgkNumberOfParticleSpecies=AliPID::kSPECIESC;
54 static const Int_t fgkNumberOfGainScenarios=3;
55 static const Int_t fgkNumberOfdEdxSourceScenarios=3;
56
57 enum ETPCdEdxSource {
58 kdEdxDefault=0, // use combined dEdx from IROC+OROC (assumes ideal detector)
59 kdEdxOROC=1, // use only OROC
60 kdEdxHybrid=2, // Use IROC+OROC dEdx only where IROCS are good (high gain), otherwise fall back to OROC only
61 kdEdxInvalid=3 //invalid
62 };
63
10d100d4 64 void SetSigma(Float_t res0, Float_t resN2);
aea7a46d 65 void SetBetheBlochParameters(Double_t kp1,
66 Double_t kp2,
67 Double_t kp3,
68 Double_t kp4,
69 Double_t kp5
70 );
1b45e564 71 //Better prevent user from setting fMIP != 50. because fMIP set fix to 50 for much other code:
10d100d4 72 void SetMip(Float_t mip) { fMIP = mip; } // Set overall normalisation; mean dE/dx for MIP
aea7a46d 73 Double_t Bethe(Double_t bg) const;
d2aa6df0 74 void SetUseDatabase(Bool_t useDatabase) { fUseDatabase = useDatabase;}
644666df 75 Bool_t GetUseDatabase() const { return fUseDatabase;}
d2aa6df0 76
77 void SetResponseFunction(AliPID::EParticleType type, TObject * const o) { fResponseFunctions.AddAt(o,(Int_t)type); }
deae51a8 78 const TObject * GetResponseFunction(AliPID::EParticleType type) { return fResponseFunctions.At((Int_t)type); }
644666df 79 void SetVoltage(Int_t n, Float_t v) {fVoltageMap[n]=v;}
80 void SetVoltageMap(const TVectorF& a) {fVoltageMap=a;} //resets ownership, ~ will not delete contents
81 Float_t GetVoltage(Int_t n) const {return fVoltageMap[n];}
82 void SetLowGainIROCthreshold(Float_t v) {fLowGainIROCthreshold=v;}
83 void SetBadIROCthreshold(Float_t v) {fBadIROCthreshhold=v;}
84 void SetLowGainOROCthreshold(Float_t v) {fLowGainOROCthreshold=v;}
85 void SetBadOROCthreshold(Float_t v) {fBadOROCthreshhold=v;}
86 void SetMaxBadLengthFraction(Float_t f) {fMaxBadLengthFraction=f;}
87
88 void SetMagField(Double_t mf) { fMagField=mf; }
d2aa6df0 89
f84b18dd 90 const TH2D* GetEtaCorrMap() const { return fhEtaCorr; };
91 Bool_t SetEtaCorrMap(TH2D* hMap);
92
87da0205 93 Double_t GetTrackTanTheta(const AliVTrack *track) const;
f84b18dd 94
87da0205 95 Double_t GetEtaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
96
f84b18dd 97 Double_t GetEtaCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
98
99 const TH2D* GetSigmaPar1Map() const { return fhEtaSigmaPar1; };
100 Double_t GetSigmaPar0() const { return fSigmaPar0; };
101 Bool_t SetSigmaParams(TH2D* hSigmaPar1Map, Double_t sigmaPar0);
102
103 Double_t GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
104
87da0205 105
106 const TF1* GetMultiplicityCorrectionFunction() const { return fCorrFuncMultiplicity; };
107 void SetParameterMultiplicityCorrection(Int_t parIndex, Double_t parValue)
108 { if (fCorrFuncMultiplicity) fCorrFuncMultiplicity->SetParameter(parIndex, parValue); };
109
110 const TF1* GetMultiplicityCorrectionFunctionTanTheta() const { return fCorrFuncMultiplicityTanTheta; };
111 void SetParameterMultiplicityCorrectionTanTheta(Int_t parIndex, Double_t parValue)
112 { if (fCorrFuncMultiplicityTanTheta) fCorrFuncMultiplicityTanTheta->SetParameter(parIndex, parValue); };
113
114 const TF1* GetMultiplicitySigmaCorrectionFunction() const { return fCorrFuncSigmaMultiplicity; };
115 void SetParameterMultiplicitySigmaCorrection(Int_t parIndex, Double_t parValue)
116 { if (fCorrFuncSigmaMultiplicity) fCorrFuncSigmaMultiplicity->SetParameter(parIndex, parValue); };
117
118 void ResetMultiplicityCorrectionFunctions();
119
120 void SetCurrentEventMultiplicity(Int_t value) { fCurrentEventMultiplicity = value; };
121 Int_t GetCurrentEventMultiplicity() const { return fCurrentEventMultiplicity; };
122
854cc597 123 Double_t GetMultiplicityCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
87da0205 124
854cc597 125 Double_t GetMultiplicitySigmaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
87da0205 126
127 Double_t GetMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
128
129 Double_t GetEtaAndMultiplicityCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species,
130 ETPCdEdxSource dedxSource = kdEdxDefault) const;
854cc597 131
132 // Fast functions for expert use only
133 Double_t GetEtaCorrectionFast(const AliVTrack *track, Double_t dEdxSplines) const;
134
205d83d8 135 Double_t GetMultiplicityCorrectionFast(const AliVTrack *track, Double_t dEdxExpected, Int_t multiplicity) const;
854cc597 136
205d83d8 137 Double_t GetMultiplicitySigmaCorrectionFast(Double_t dEdxExpected, Int_t multiplicity) const;
854cc597 138
139 Double_t GetSigmaPar1Fast(const AliVTrack *track, AliPID::EParticleType species,
140 Double_t dEdx, const TSpline3* responseFunction) const;
141
644666df 142 //NEW
143 void SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario );
644666df 144 Double_t GetExpectedSignal( const AliVTrack* track,
145 AliPID::EParticleType species,
f84b18dd 146 ETPCdEdxSource dedxSource = kdEdxDefault,
87da0205 147 Bool_t correctEta = kFALSE,
148 Bool_t correctMultiplicity = kFALSE) const;
644666df 149 Double_t GetExpectedSigma( const AliVTrack* track,
150 AliPID::EParticleType species,
f84b18dd 151 ETPCdEdxSource dedxSource = kdEdxDefault,
87da0205 152 Bool_t correctEta = kFALSE,
153 Bool_t correctMultiplicity = kFALSE) const;
644666df 154 Float_t GetNumberOfSigmas( const AliVTrack* track,
155 AliPID::EParticleType species,
f84b18dd 156 ETPCdEdxSource dedxSource = kdEdxDefault,
87da0205 157 Bool_t correctEta = kFALSE,
158 Bool_t correctMultiplicity = kFALSE) const;
567624b5 159
160 Float_t GetSignalDelta( const AliVTrack* track,
161 AliPID::EParticleType species,
162 ETPCdEdxSource dedxSource = kdEdxDefault,
87da0205 163 Bool_t correctEta = kFALSE,
164 Bool_t correctMultiplicity = kFALSE,
165 Bool_t ratio = kFALSE) const;
567624b5 166
644666df 167 void SetResponseFunction(TObject* o,
168 AliPID::EParticleType type,
169 ETPCgainScenario gainScenario);
170 void Print(Option_t* option="") const;
171 TSpline3* GetResponseFunction( AliPID::EParticleType species,
172 ETPCgainScenario gainScenario ) const;
173 TSpline3* GetResponseFunction( const AliVTrack* track,
174 AliPID::EParticleType species,
f84b18dd 175 ETPCdEdxSource dedxSource = kdEdxDefault) const;
644666df 176 Bool_t ResponseFunctiondEdxN(const AliVTrack* track,
177 AliPID::EParticleType species,
f84b18dd 178 ETPCdEdxSource dedxSource,
179 Double_t& dEdx, Int_t& nPoints, ETPCgainScenario& gainScenario, TSpline3** responseFunction) const;
f85a3764 180 Bool_t sectorNumbersInOut(Double_t* trackPositionInner,
181 Double_t* trackPositionOuter,
644666df 182 Float_t& phiIn, Float_t& phiOut,
183 Int_t& in, Int_t& out ) const;
184 AliTPCPIDResponse::EChamberStatus TrackStatus(const AliVTrack* track, Int_t layer) const;
185 Float_t MaxClusterRadius(const AliVTrack* track) const;
186 Bool_t TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const;
854cc597 187 static const char* GainScenarioName(Int_t n) {return fgkGainScenarioName[(n>fgkNumberOfGainScenarios)?fgkNumberOfGainScenarios:n];}
644666df 188 Int_t ResponseFunctionIndex( AliPID::EParticleType species,
189 ETPCgainScenario gainScenario ) const;
190 void ResetSplines();
191
644666df 192 //OLD
205d83d8 193 Double_t GetExpectedSignal(Float_t mom,
aea7a46d 194 AliPID::EParticleType n=AliPID::kKaon) const;
205d83d8 195 Double_t GetExpectedSigma(Float_t mom, Int_t nPoints,
644666df 196 AliPID::EParticleType n=AliPID::kKaon) const;
205d83d8 197 Float_t GetNumberOfSigmas(Float_t mom,
198 Float_t dEdx,
199 Int_t nPoints,
644666df 200 AliPID::EParticleType n=AliPID::kKaon) const {
f84b18dd 201 //
202 // Deprecated function (for backward compatibility). Please use
87da0205 203 // GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource,
204 // Bool_t correctEta, Bool_t correctMultiplicity)
f84b18dd 205 // instead!TODO
206 //
207
10d100d4 208 Double_t bethe=GetExpectedSignal(mom,n);
209 Double_t sigma=GetExpectedSigma(mom,nPoints,n);
210 return (dEdx-bethe)/sigma;
211 }
aea7a46d 212
10d100d4 213 Double_t GetMIP() const { return fMIP;}
644666df 214 Float_t GetRes0() const { return fRes0[0]; }
215 Float_t GetResN2() const { return fResN2[0]; }
216 Float_t GetRes0(ETPCgainScenario s) const { return fRes0[s]; }
217 Float_t GetResN2(ETPCgainScenario s) const { return fResN2[s]; }
b981edcd 218
912f7fdc 219 Bool_t RegisterSpline(const char * name, Int_t index);
220 Double_t EvaldEdxSpline(Double_t bg,Int_t entry);
221 static Double_t SEvaldEdx(Double_t bg,Int_t entry){ return (fgInstance!=0)? fgInstance->EvaldEdxSpline(bg,entry):0;};
222
f84b18dd 223protected:
224 Double_t GetExpectedSignal(const AliVTrack* track,
225 AliPID::EParticleType species,
226 Double_t dEdx,
227 const TSpline3* responseFunction,
87da0205 228 Bool_t correctEta,
229 Bool_t correctMultiplicity) const;
f84b18dd 230
231 Double_t GetExpectedSigma(const AliVTrack* track,
232 AliPID::EParticleType species,
233 ETPCgainScenario gainScenario,
234 Double_t dEdx,
235 Int_t nPoints,
236 const TSpline3* responseFunction,
87da0205 237 Bool_t correctEta,
238 Bool_t correctMultiplicity) const;
912f7fdc 239 //
240 // function for numberical debugging 0 registed splines can be used in the TFormula and tree visualizations
241 //
8c6a71ab 242private:
10d100d4 243 Float_t fMIP; // dEdx for MIP
644666df 244 Float_t fRes0[fgkNumberOfGainScenarios]; // relative dEdx resolution rel sigma = fRes0*sqrt(1+fResN2/npoint)
245 Float_t fResN2[fgkNumberOfGainScenarios]; // relative Npoint dependence rel sigma = fRes0*sqrt(1+fResN2/npoint)
aea7a46d 246
247 Double_t fKp1; // Parameters
248 Double_t fKp2; // of
249 Double_t fKp3; // the ALEPH
250 Double_t fKp4; // Bethe-Bloch
251 Double_t fKp5; // formula
252
d2aa6df0 253 Bool_t fUseDatabase; // flag if fine-tuned database-response or simple ALEPH BB should be used
644666df 254
d2aa6df0 255 TObjArray fResponseFunctions; //! ObjArray of response functions individually for each particle
644666df 256 TVectorF fVoltageMap; //!stores a map of voltages wrt nominal for all chambers
257 Float_t fLowGainIROCthreshold; //voltage threshold below which the IROC is considered low gain
258 Float_t fBadIROCthreshhold; //voltage threshold for bad IROCS
259 Float_t fLowGainOROCthreshold; //voltage threshold below which the OROC is considered low gain
260 Float_t fBadOROCthreshhold; //voltage threshold for bad OROCS
261 Float_t fMaxBadLengthFraction; //the maximum allowed fraction of track length in a bad sector.
262
644666df 263 Int_t sectorNumber(Double_t phi) const;
264
265 Double_t fMagField; //! Magnetic field
266
267 static const char* fgkGainScenarioName[fgkNumberOfGainScenarios+1];
d2aa6df0 268
f84b18dd 269 TH2D* fhEtaCorr; //! Map for TPC eta correction
270 TH2D* fhEtaSigmaPar1; //! Map for parameter 1 of the dEdx sigma parametrisation
271
272 Double_t fSigmaPar0; // Parameter 0 of the dEdx sigma parametrisation
87da0205 273
274 Int_t fCurrentEventMultiplicity; // Multiplicity of the current event
275 TF1* fCorrFuncMultiplicity; //! Function to correct for the multiplicity dependence of the TPC dEdx
276 TF1* fCorrFuncMultiplicityTanTheta; //! Function to correct the additional tanTheta dependence of the multiplicity dependence of the TPC dEdx
277 TF1* fCorrFuncSigmaMultiplicity; //! Function to correct for the multiplicity dependence of the TPC dEdx resolution
f84b18dd 278
912f7fdc 279 //
280 //
281 static AliTPCPIDResponse* fgInstance; //! Instance of this class (singleton implementation)
282 TObjArray fSplineArray; //array of registered splines
87da0205 283 ClassDef(AliTPCPIDResponse,6) // TPC PID class
8c6a71ab 284};
285
912f7fdc 286
8c6a71ab 287#endif
288
289