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   //TODO Remove? 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   //Better prevent user from setting fMIP != 50. because fMIP set fix to 50 for much other code:
  71   void SetMip(Float_t mip) { fMIP = mip; } // Set overall normalisation; mean dE/dx for MIP
  72   Double_t Bethe(Double_t bg) const;
  73   void SetUseDatabase(Bool_t useDatabase) { fUseDatabase = useDatabase;}
  74   Bool_t GetUseDatabase() const { return fUseDatabase;}
  75
  76   void SetResponseFunction(AliPID::EParticleType type, TObject * const o) { fResponseFunctions.AddAt(o,(Int_t)type); }
  77   const TObject * GetResponseFunction(AliPID::EParticleType type) { return fResponseFunctions.At((Int_t)type); }
  78   void SetVoltage(Int_t n, Float_t v) {fVoltageMap[n]=v;}
  79   void SetVoltageMap(const TVectorF& a) {fVoltageMap=a;} //resets ownership, ~ will not delete contents
  80   Float_t GetVoltage(Int_t n) const {return fVoltageMap[n];}
  81   void SetLowGainIROCthreshold(Float_t v) {fLowGainIROCthreshold=v;}
  82   void SetBadIROCthreshold(Float_t v) {fBadIROCthreshhold=v;}
  83   void SetLowGainOROCthreshold(Float_t v) {fLowGainOROCthreshold=v;}
  84   void SetBadOROCthreshold(Float_t v) {fBadOROCthreshhold=v;}
  85   void SetMaxBadLengthFraction(Float_t f) {fMaxBadLengthFraction=f;}
  86
  87   void SetMagField(Double_t mf) { fMagField=mf; }
  88
  89   const TH2D* GetEtaCorrMap() const { return fhEtaCorr; };
  90   Bool_t SetEtaCorrMap(TH2D* hMap);
  91
  92   Double_t GetEtaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
  93
  94   Double_t GetEtaCorrectedTrackdEdx(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
  95
  96   const TH2D* GetSigmaPar1Map() const { return fhEtaSigmaPar1; };
  97   Double_t GetSigmaPar0() const { return fSigmaPar0; };
  98   Bool_t SetSigmaParams(TH2D* hSigmaPar1Map, Double_t sigmaPar0);
  99
 100   Double_t GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
 101
 102   //NEW
 103   void SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario );
 104   Double_t GetExpectedSignal( const AliVTrack* track,
 105                               AliPID::EParticleType species,
 106                               ETPCdEdxSource dedxSource = kdEdxDefault,
 107                               Bool_t correctEta = kFALSE) const;//TODO: In future, default kTRUE
 108   Double_t GetExpectedSigma( const AliVTrack* track,
 109                              AliPID::EParticleType species,
 110                              ETPCdEdxSource dedxSource = kdEdxDefault,
 111                              Bool_t correctEta = kFALSE) const;//TODO: In future, default kTRUE
 112   Float_t GetNumberOfSigmas( const AliVTrack* track,
 113                              AliPID::EParticleType species,
 114                              ETPCdEdxSource dedxSource = kdEdxDefault,
 115                              Bool_t correctEta = kFALSE) const;//TODO: In future, default kTRUE
 116
 117   Float_t GetSignalDelta( const AliVTrack* track,
 118                           AliPID::EParticleType species,
 119                           ETPCdEdxSource dedxSource = kdEdxDefault,
 120                           Bool_t correctEta = kFALSE, Bool_t ratio=kFALSE) const;
 121
 122   void SetResponseFunction(TObject* o,
 123                            AliPID::EParticleType type,
 124                            ETPCgainScenario gainScenario);
 125   void Print(Option_t* option="") const;
 126   TSpline3* GetResponseFunction( AliPID::EParticleType species,
 127                                  ETPCgainScenario gainScenario ) const;
 128   TSpline3* GetResponseFunction( const AliVTrack* track,
 129                                  AliPID::EParticleType species,
 130                                  ETPCdEdxSource dedxSource = kdEdxDefault) const;
 131   Bool_t ResponseFunctiondEdxN(const AliVTrack* track,
 132                                AliPID::EParticleType species,
 133                                ETPCdEdxSource dedxSource,
 134                                Double_t& dEdx, Int_t& nPoints, ETPCgainScenario& gainScenario, TSpline3** responseFunction) const;
 135   Bool_t sectorNumbersInOut(Double_t* trackPositionInner,
 136                             Double_t* trackPositionOuter,
 137                             Float_t& phiIn, Float_t& phiOut,
 138                             Int_t& in, Int_t& out ) const;
 139   AliTPCPIDResponse::EChamberStatus TrackStatus(const AliVTrack* track, Int_t layer) const;
 140   Float_t MaxClusterRadius(const AliVTrack* track) const;
 141   Bool_t TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const;
 142   static const char* GainScenarioName(Int_t n) {return fgkGainScenarioName[(n>fgkNumberOfdEdxSourceScenarios)?fgkNumberOfdEdxSourceScenarios+1:n];}
 143   Int_t ResponseFunctionIndex( AliPID::EParticleType species,
 144                                ETPCgainScenario gainScenario ) const;
 145   void ResetSplines();
 146
 147   //OLD
 148   Double_t GetExpectedSignal(const Float_t mom,
 149                      AliPID::EParticleType n=AliPID::kKaon) const;
 150   Double_t GetExpectedSigma(const Float_t mom, const Int_t nPoints,
 151                             AliPID::EParticleType n=AliPID::kKaon) const;
 152   Float_t  GetNumberOfSigmas(const Float_t mom,
 153                              const Float_t dEdx,
 154                                                const Int_t nPoints,
 155                              AliPID::EParticleType n=AliPID::kKaon) const {
 156     //
 157     // Deprecated function (for backward compatibility). Please use
 158     // GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource )
 159     // instead!TODO
 160     //
 161
 162     Double_t bethe=GetExpectedSignal(mom,n);
 163     Double_t sigma=GetExpectedSigma(mom,nPoints,n);
 164     return (dEdx-bethe)/sigma;
 165   }
 166
 167   Double_t GetMIP() const { return fMIP;}
 168   Float_t  GetRes0()  const { return fRes0[0];  }
 169   Float_t  GetResN2() const { return fResN2[0]; }
 170   Float_t  GetRes0(ETPCgainScenario s)  const { return fRes0[s];  }
 171   Float_t  GetResN2(ETPCgainScenario s) const { return fResN2[s]; }
 172
 173 protected:
 174   Double_t GetExpectedSignal(const AliVTrack* track,
 175                              AliPID::EParticleType species,
 176                              Double_t dEdx,
 177                              const TSpline3* responseFunction,
 178                              Bool_t correctEta) const;
 179
 180   Double_t GetExpectedSigma(const AliVTrack* track,
 181                             AliPID::EParticleType species,
 182                             ETPCgainScenario gainScenario,
 183                             Double_t dEdx,
 184                             Int_t nPoints,
 185                             const TSpline3* responseFunction,
 186                             Bool_t correctEta) const;
 187
 188   Double_t GetEtaCorrection(const AliVTrack *track, Double_t dEdxSplines) const;
 189
 190   Double_t GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species,
 191                         Double_t dEdx, const TSpline3* responseFunction) const;
 192
 193 private:
 194   Float_t fMIP;          // dEdx for MIP
 195   Float_t fRes0[fgkNumberOfGainScenarios];  // relative dEdx resolution  rel sigma = fRes0*sqrt(1+fResN2/npoint)
 196   Float_t fResN2[fgkNumberOfGainScenarios]; // relative Npoint dependence rel  sigma = fRes0*sqrt(1+fResN2/npoint)
 197
 198   Double_t fKp1;   // Parameters
 199   Double_t fKp2;   //    of
 200   Double_t fKp3;   // the ALEPH
 201   Double_t fKp4;   // Bethe-Bloch
 202   Double_t fKp5;   // formula
 203
 204   Bool_t fUseDatabase; // flag if fine-tuned database-response or simple ALEPH BB should be used
 205
 206   TObjArray fResponseFunctions; //! ObjArray of response functions individually for each particle
 207   TVectorF fVoltageMap; //!stores a map of voltages wrt nominal for all chambers
 208   Float_t fLowGainIROCthreshold;  //voltage threshold below which the IROC is considered low gain
 209   Float_t fBadIROCthreshhold;     //voltage threshold for bad IROCS
 210   Float_t fLowGainOROCthreshold;  //voltage threshold below which the OROC is considered low gain
 211   Float_t fBadOROCthreshhold;     //voltage threshold for bad OROCS
 212   Float_t fMaxBadLengthFraction;  //the maximum allowed fraction of track length in a bad sector.
 213
 214   Int_t sectorNumber(Double_t phi) const;
 215
 216   Double_t fMagField;  //! Magnetic field
 217
 218   static const char* fgkGainScenarioName[fgkNumberOfGainScenarios+1];
 219
 220   TH2D* fhEtaCorr; //! Map for TPC eta correction
 221   TH2D* fhEtaSigmaPar1; //! Map for parameter 1 of the dEdx sigma parametrisation
 222
 223   Double_t fSigmaPar0; // Parameter 0 of the dEdx sigma parametrisation
 224
 225   ClassDef(AliTPCPIDResponse,5)   // TPC PID class
 226 };
 227
 228 #endif
 229
 230