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 #include <TF1.h>
  24
  25 #include "AliPID.h"
  26 #include "AliVTrack.h"
  27
  28 class TH2D;
  29 class TSpline3;
  30
  31 class AliTPCPIDResponse: public TNamed {
  32 public:
  33   AliTPCPIDResponse();
  34   //TODO Remove? AliTPCPIDResponse(const Double_t *param);
  35   AliTPCPIDResponse(const AliTPCPIDResponse&);
  36   AliTPCPIDResponse& operator=(const AliTPCPIDResponse&);
  37   virtual ~AliTPCPIDResponse();
  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
  64   void SetSigma(Float_t res0, Float_t resN2);
  65   void SetBetheBlochParameters(Double_t kp1,
  66                                Double_t kp2,
  67                                Double_t kp3,
  68                                Double_t kp4,
  69                                Double_t kp5
  70                                );
  71   //Better prevent user from setting fMIP != 50. because fMIP set fix to 50 for much other code:
  72   void SetMip(Float_t mip) { fMIP = mip; } // Set overall normalisation; mean dE/dx for MIP
  73   Double_t Bethe(Double_t bg) const;
  74   void SetUseDatabase(Bool_t useDatabase) { fUseDatabase = useDatabase;}
  75   Bool_t GetUseDatabase() const { return fUseDatabase;}
  76
  77   void SetResponseFunction(AliPID::EParticleType type, TObject * const o) { fResponseFunctions.AddAt(o,(Int_t)type); }
  78   const TObject * GetResponseFunction(AliPID::EParticleType type) { return fResponseFunctions.At((Int_t)type); }
  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; }
  89
  90   const TH2D* GetEtaCorrMap() const { return fhEtaCorr; };
  91   Bool_t SetEtaCorrMap(TH2D* hMap);
  92
  93   Double_t GetTrackTanTheta(const AliVTrack *track) const;
  94
  95   Double_t GetEtaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
  96
  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
 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
 123   Double_t GetMultiplicityCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
 124
 125   Double_t GetMultiplicitySigmaCorrection(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource = kdEdxDefault) const;
 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;
 131
 132   // Fast functions for expert use only
 133   Double_t GetEtaCorrectionFast(const AliVTrack *track, Double_t dEdxSplines) const;
 134
 135   Double_t GetMultiplicityCorrectionFast(const AliVTrack *track, Double_t dEdxExpected, Int_t multiplicity) const;
 136
 137   Double_t GetMultiplicitySigmaCorrectionFast(Double_t dEdxExpected, Int_t multiplicity) const;
 138
 139   Double_t GetSigmaPar1Fast(const AliVTrack *track, AliPID::EParticleType species,
 140                             Double_t dEdx, const TSpline3* responseFunction) const;
 141
 142   //NEW
 143   void SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario );
 144   Double_t GetExpectedSignal( const AliVTrack* track,
 145                               AliPID::EParticleType species,
 146                               ETPCdEdxSource dedxSource = kdEdxDefault,
 147                               Bool_t correctEta = kFALSE,
 148                               Bool_t correctMultiplicity = kFALSE) const;
 149   Double_t GetExpectedSigma( const AliVTrack* track,
 150                              AliPID::EParticleType species,
 151                              ETPCdEdxSource dedxSource = kdEdxDefault,
 152                              Bool_t correctEta = kFALSE,
 153                              Bool_t correctMultiplicity = kFALSE) const;
 154   Float_t GetNumberOfSigmas( const AliVTrack* track,
 155                              AliPID::EParticleType species,
 156                              ETPCdEdxSource dedxSource = kdEdxDefault,
 157                              Bool_t correctEta = kFALSE,
 158                              Bool_t correctMultiplicity = kFALSE) const;
 159
 160   Float_t GetSignalDelta( const AliVTrack* track,
 161                           AliPID::EParticleType species,
 162                           ETPCdEdxSource dedxSource = kdEdxDefault,
 163                           Bool_t correctEta = kFALSE,
 164                           Bool_t correctMultiplicity = kFALSE,
 165                           Bool_t ratio = kFALSE) const;
 166
 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,
 175                                  ETPCdEdxSource dedxSource = kdEdxDefault) const;
 176   Bool_t ResponseFunctiondEdxN(const AliVTrack* track,
 177                                AliPID::EParticleType species,
 178                                ETPCdEdxSource dedxSource,
 179                                Double_t& dEdx, Int_t& nPoints, ETPCgainScenario& gainScenario, TSpline3** responseFunction) const;
 180   Bool_t sectorNumbersInOut(Double_t* trackPositionInner,
 181                             Double_t* trackPositionOuter,
 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;
 187   static const char* GainScenarioName(Int_t n) {return fgkGainScenarioName[(n>fgkNumberOfGainScenarios)?fgkNumberOfGainScenarios:n];}
 188   Int_t ResponseFunctionIndex( AliPID::EParticleType species,
 189                                ETPCgainScenario gainScenario ) const;
 190   void ResetSplines();
 191
 192   //OLD
 193   Double_t GetExpectedSignal(Float_t mom,
 194                      AliPID::EParticleType n=AliPID::kKaon) const;
 195   Double_t GetExpectedSigma(Float_t mom, Int_t nPoints,
 196                             AliPID::EParticleType n=AliPID::kKaon) const;
 197   Float_t  GetNumberOfSigmas(Float_t mom,
 198                              Float_t dEdx,
 199                              Int_t nPoints,
 200                              AliPID::EParticleType n=AliPID::kKaon) const {
 201     //
 202     // Deprecated function (for backward compatibility). Please use
 203     // GetNumberOfSigmas(const AliVTrack *track, AliPID::EParticleType species, ETPCdEdxSource dedxSource,
 204     // Bool_t correctEta, Bool_t correctMultiplicity)
 205     // instead!TODO
 206     //
 207
 208     Double_t bethe=GetExpectedSignal(mom,n);
 209     Double_t sigma=GetExpectedSigma(mom,nPoints,n);
 210     return (dEdx-bethe)/sigma;
 211   }
 212
 213   Double_t GetMIP() const { return fMIP;}
 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]; }
 218
 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
 223 protected:
 224   Double_t GetExpectedSignal(const AliVTrack* track,
 225                              AliPID::EParticleType species,
 226                              Double_t dEdx,
 227                              const TSpline3* responseFunction,
 228                              Bool_t correctEta,
 229                              Bool_t correctMultiplicity) const;
 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,
 237                             Bool_t correctEta,
 238                             Bool_t correctMultiplicity) const;
 239
 240   Double_t GetMultiplicityCorrection(const AliVTrack *track, const Double_t dEdxExpected, const Int_t multiplicity) const;
 241
 242   Double_t GetMultiplicitySigmaCorrection(const Double_t dEdxExpected, const Int_t multiplicity) const;
 243
 244   Double_t GetSigmaPar1(const AliVTrack *track, AliPID::EParticleType species,
 245                         Double_t dEdx, const TSpline3* responseFunction) const;
 246   //
 247   // function for numberical debugging 0 registed splines can be used in the TFormula and tree visualizations
 248   //
 249 private:
 250   Float_t fMIP;          // dEdx for MIP
 251   Float_t fRes0[fgkNumberOfGainScenarios];  // relative dEdx resolution  rel sigma = fRes0*sqrt(1+fResN2/npoint)
 252   Float_t fResN2[fgkNumberOfGainScenarios]; // relative Npoint dependence rel  sigma = fRes0*sqrt(1+fResN2/npoint)
 253
 254   Double_t fKp1;   // Parameters
 255   Double_t fKp2;   //    of
 256   Double_t fKp3;   // the ALEPH
 257   Double_t fKp4;   // Bethe-Bloch
 258   Double_t fKp5;   // formula
 259
 260   Bool_t fUseDatabase; // flag if fine-tuned database-response or simple ALEPH BB should be used
 261
 262   TObjArray fResponseFunctions; //! ObjArray of response functions individually for each particle
 263   TVectorF fVoltageMap; //!stores a map of voltages wrt nominal for all chambers
 264   Float_t fLowGainIROCthreshold;  //voltage threshold below which the IROC is considered low gain
 265   Float_t fBadIROCthreshhold;     //voltage threshold for bad IROCS
 266   Float_t fLowGainOROCthreshold;  //voltage threshold below which the OROC is considered low gain
 267   Float_t fBadOROCthreshhold;     //voltage threshold for bad OROCS
 268   Float_t fMaxBadLengthFraction;  //the maximum allowed fraction of track length in a bad sector.
 269
 270   Int_t sectorNumber(Double_t phi) const;
 271
 272   Double_t fMagField;  //! Magnetic field
 273
 274   static const char* fgkGainScenarioName[fgkNumberOfGainScenarios+1];
 275
 276   TH2D* fhEtaCorr; //! Map for TPC eta correction
 277   TH2D* fhEtaSigmaPar1; //! Map for parameter 1 of the dEdx sigma parametrisation
 278
 279   Double_t fSigmaPar0; // Parameter 0 of the dEdx sigma parametrisation
 280
 281   Int_t fCurrentEventMultiplicity; // Multiplicity of the current event
 282   TF1* fCorrFuncMultiplicity; //! Function to correct for the multiplicity dependence of the TPC dEdx
 283   TF1* fCorrFuncMultiplicityTanTheta; //! Function to correct the additional tanTheta dependence of the multiplicity dependence of the TPC dEdx
 284   TF1* fCorrFuncSigmaMultiplicity; //! Function to correct for the multiplicity dependence of the TPC dEdx resolution
 285
 286   //
 287   //
 288   static AliTPCPIDResponse*   fgInstance;     //! Instance of this class (singleton implementation)
 289   TObjArray                   fSplineArray;   //array of registered splines
 290   ClassDef(AliTPCPIDResponse,6)   // TPC PID class
 291 };
 292
 293
 294 #endif
 295
 296