[u/mrichter/AliRoot.git] / STEER / STEERBase / AliTPCPIDResponse.h

#ifndef ALITPCPIDRESPONSE_H
#define ALITPCPIDRESPONSE_H
/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 * See cxx source for full Copyright notice                               */

//-------------------------------------------------------
//                    TPC PID class
// A very naive design... Should be made better by the detector experts...
//   Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch 
// With many additions and modifications suggested by
//      Alexander Kalweit, GSI, alexander.philipp.kalweit@cern.ch
//      Dariusz Miskowiec, GSI, D.Miskowiec@gsi.de
//-------------------------------------------------------
#include <Rtypes.h>

#include <TNamed.h>
#include <TVectorF.h>
#include <TObjArray.h>

#include "AliPID.h"

class AliVTrack;
class TSpline3;

class AliTPCPIDResponse: public TNamed {
public:
  AliTPCPIDResponse();
  AliTPCPIDResponse(const Double_t *param);
  AliTPCPIDResponse(const AliTPCPIDResponse&);
  AliTPCPIDResponse& operator=(const AliTPCPIDResponse&);
  virtual ~AliTPCPIDResponse() {}

  enum EChamberStatus {
    kChamberOff=0,
    kChamberHighGain=1,
    kChamberLowGain=2,
    kChamberInvalid=3
  };
  
  enum ETPCgainScenario {
    kDefault= 0,
    kALLhigh = 1,
    kOROChigh = 2,
    kGainScenarioInvalid = 3
  };

  static const Int_t fgkNumberOfParticleSpecies=AliPID::kSPECIESC;
  static const Int_t fgkNumberOfGainScenarios=3;
  static const Int_t fgkNumberOfdEdxSourceScenarios=3;

  enum ETPCdEdxSource {
    kdEdxDefault=0,        // use combined dEdx from IROC+OROC (assumes ideal detector)
    kdEdxOROC=1,       // use only OROC
    kdEdxHybrid=2,   // Use IROC+OROC dEdx only where IROCS are good (high gain), otherwise fall back to OROC only
    kdEdxInvalid=3     //invalid
  };

  void SetSigma(Float_t res0, Float_t resN2);
  void SetBetheBlochParameters(Double_t kp1,
                               Double_t kp2,
                               Double_t kp3,
                               Double_t kp4,
                               Double_t kp5
                               );
  void SetMip(Float_t mip) { fMIP = mip; } // Set overall normalisation; mean dE/dx for MIP
  Double_t Bethe(Double_t bg) const;
  void SetUseDatabase(Bool_t useDatabase) { fUseDatabase = useDatabase;}
  Bool_t GetUseDatabase() const { return fUseDatabase;}
  
  void SetResponseFunction(AliPID::EParticleType type, TObject * const o) { fResponseFunctions.AddAt(o,(Int_t)type); }
  const TObject * GetResponseFunction(AliPID::EParticleType type) { return fResponseFunctions.At((Int_t)type); }
  void SetVoltage(Int_t n, Float_t v) {fVoltageMap[n]=v;}
  void SetVoltageMap(const TVectorF& a) {fVoltageMap=a;} //resets ownership, ~ will not delete contents
  Float_t GetVoltage(Int_t n) const {return fVoltageMap[n];}
  void SetLowGainIROCthreshold(Float_t v) {fLowGainIROCthreshold=v;}
  void SetBadIROCthreshold(Float_t v) {fBadIROCthreshhold=v;}
  void SetLowGainOROCthreshold(Float_t v) {fLowGainOROCthreshold=v;}
  void SetBadOROCthreshold(Float_t v) {fBadOROCthreshhold=v;}
  void SetMaxBadLengthFraction(Float_t f) {fMaxBadLengthFraction=f;}

  void SetMagField(Double_t mf) { fMagField=mf; }
  
  //NEW
  void SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario );
  Double_t GetExpectedSignal( Double_t momentum,
                              AliPID::EParticleType species,
                              const TSpline3* responseFunction ) const;
  Double_t GetExpectedSignal( const AliVTrack* track,
                              AliPID::EParticleType species,
                              ETPCdEdxSource dedxSource );
  Double_t GetExpectedSigma( const AliVTrack* track, 
                             AliPID::EParticleType species,
                             ETPCdEdxSource dedxSource );
  Double_t GetExpectedSigma( Double_t mom,
                             Int_t nPoints,
                             AliPID::EParticleType species,
                             ETPCgainScenario gainScenario,
                             const TSpline3* responseFunction) const;
  Float_t GetNumberOfSigmas( const AliVTrack* track,
                             AliPID::EParticleType species,
                             ETPCdEdxSource dedxSource );

  void SetResponseFunction(TObject* o,
                           AliPID::EParticleType type,
                           ETPCgainScenario gainScenario);
  void Print(Option_t* option="") const;
  TSpline3* GetResponseFunction( AliPID::EParticleType species,
                                 ETPCgainScenario gainScenario ) const;
  TSpline3* GetResponseFunction( const AliVTrack* track,
                                 AliPID::EParticleType species,
                                 ETPCdEdxSource dedxSource );
  Bool_t ResponseFunctiondEdxN(const AliVTrack* track, 
                               AliPID::EParticleType species,
                               ETPCdEdxSource dedxSource);
  Bool_t sectorNumbersInOut(const AliVTrack* track, 
                            Double_t innerRadius, Double_t outerRadius, 
                            Float_t& phiIn, Float_t& phiOut, 
                            Int_t& in, Int_t& out ) const;
  AliTPCPIDResponse::EChamberStatus TrackStatus(const AliVTrack* track, Int_t layer) const;
  Float_t MaxClusterRadius(const AliVTrack* track) const;
  Bool_t TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const;
  static const char* GainScenarioName(Int_t n) {return fgkGainScenarioName[(n>fgkNumberOfdEdxSourceScenarios)?fgkNumberOfdEdxSourceScenarios+1:n];}
  Int_t ResponseFunctionIndex( AliPID::EParticleType species,
                               ETPCgainScenario gainScenario ) const;
  void ResetSplines();

  void InvalidateCurrentValues();
  TSpline3* GetCurrentResponseFunction() const {return fCurrentResponseFunction;}
  Double_t GetCurrentdEdx() const {return fCurrentdEdx;}
  Int_t GetCurrentNPoints() const {return fCurrentNPoints;}
  ETPCgainScenario GetCurrentGainScenario() const {return fCurrentGainScenario;}

  //OLD
  Double_t GetExpectedSignal(const Float_t mom,
                     AliPID::EParticleType n=AliPID::kKaon) const;
  Double_t GetExpectedSigma(const Float_t mom, const Int_t nPoints,
                            AliPID::EParticleType n=AliPID::kKaon) const;
  Float_t  GetNumberOfSigmas(const Float_t mom, 
                             const Float_t dEdx, 
			                       const Int_t nPoints,
                             AliPID::EParticleType n=AliPID::kKaon) const {

    Double_t bethe=GetExpectedSignal(mom,n);
    Double_t sigma=GetExpectedSigma(mom,nPoints,n);
    return (dEdx-bethe)/sigma;
  }

  Double_t GetMIP() const { return fMIP;} 
  Float_t  GetRes0()  const { return fRes0[0];  }
  Float_t  GetResN2() const { return fResN2[0]; }
  Float_t  GetRes0(ETPCgainScenario s)  const { return fRes0[s];  }
  Float_t  GetResN2(ETPCgainScenario s) const { return fResN2[s]; }

private:
  Float_t fMIP;          // dEdx for MIP
  Float_t fRes0[fgkNumberOfGainScenarios];  // relative dEdx resolution  rel sigma = fRes0*sqrt(1+fResN2/npoint)
  Float_t fResN2[fgkNumberOfGainScenarios]; // relative Npoint dependence rel  sigma = fRes0*sqrt(1+fResN2/npoint)

  Double_t fKp1;   // Parameters
  Double_t fKp2;   //    of
  Double_t fKp3;   // the ALEPH
  Double_t fKp4;   // Bethe-Bloch
  Double_t fKp5;   // formula

  Bool_t fUseDatabase; // flag if fine-tuned database-response or simple ALEPH BB should be used
  
  TObjArray fResponseFunctions; //! ObjArray of response functions individually for each particle
  TVectorF fVoltageMap; //!stores a map of voltages wrt nominal for all chambers
  Float_t fLowGainIROCthreshold;  //voltage threshold below which the IROC is considered low gain
  Float_t fBadIROCthreshhold;     //voltage threshold for bad IROCS
  Float_t fLowGainOROCthreshold;  //voltage threshold below which the OROC is considered low gain
  Float_t fBadOROCthreshhold;     //voltage threshold for bad OROCS
  Float_t fMaxBadLengthFraction;  //the maximum allowed fraction of track length in a bad sector.

  TSpline3* fCurrentResponseFunction;      //!response function for current track
  Double_t fCurrentdEdx;                  //!dEdx for currently processed track
  Int_t fCurrentNPoints;                  //!number of points used for dEdx calculation for current track
  ETPCgainScenario fCurrentGainScenario;  //!gain scenario used for current track
  Int_t sectorNumber(Double_t phi) const;

  Double_t fMagField;  //! Magnetic field

  static const char* fgkGainScenarioName[fgkNumberOfGainScenarios+1];

  ClassDef(AliTPCPIDResponse,4)   // TPC PID class
};

#endif
Commit	Line	Data
10d100d4	1	#ifndef ALITPCPIDRESPONSE_H
10d100d4	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
8c6a71ab	13	//-------------------------------------------------------
8c6a71ab	14	#include <Rtypes.h>
644666df	15
	16	#include <TNamed.h>
	17	#include <TVectorF.h>
d2aa6df0	18	#include <TObjArray.h>
8c6a71ab	19
aea7a46d	20	#include "AliPID.h"
aea7a46d	21
644666df	22	class AliVTrack;
	23	class TSpline3;
	24
	25	class AliTPCPIDResponse: public TNamed {
8c6a71ab	26	public:
10d100d4	27	AliTPCPIDResponse();
d2aa6df0	28	AliTPCPIDResponse(const Double_t *param);
644666df	29	AliTPCPIDResponse(const AliTPCPIDResponse&);
644666df	30	AliTPCPIDResponse& operator=(const AliTPCPIDResponse&);
10d100d4	31	virtual ~AliTPCPIDResponse() {}
644666df	32
	33	enum EChamberStatus {
	34	kChamberOff=0,
	35	kChamberHighGain=1,
	36	kChamberLowGain=2,
	37	kChamberInvalid=3
	38	};
	39
	40	enum ETPCgainScenario {
	41	kDefault= 0,
	42	kALLhigh = 1,
	43	kOROChigh = 2,
	44	kGainScenarioInvalid = 3
	45	};
	46
	47	static const Int_t fgkNumberOfParticleSpecies=AliPID::kSPECIESC;
	48	static const Int_t fgkNumberOfGainScenarios=3;
	49	static const Int_t fgkNumberOfdEdxSourceScenarios=3;
	50
	51	enum ETPCdEdxSource {
	52	kdEdxDefault=0, // use combined dEdx from IROC+OROC (assumes ideal detector)
	53	kdEdxOROC=1, // use only OROC
	54	kdEdxHybrid=2, // Use IROC+OROC dEdx only where IROCS are good (high gain), otherwise fall back to OROC only
	55	kdEdxInvalid=3 //invalid
	56	};
	57
10d100d4	58	void SetSigma(Float_t res0, Float_t resN2);
aea7a46d	59	void SetBetheBlochParameters(Double_t kp1,
	60	Double_t kp2,
	61	Double_t kp3,
	62	Double_t kp4,
	63	Double_t kp5
	64	);
10d100d4	65	void SetMip(Float_t mip) { fMIP = mip; } // Set overall normalisation; mean dE/dx for MIP
aea7a46d	66	Double_t Bethe(Double_t bg) const;
d2aa6df0	67	void SetUseDatabase(Bool_t useDatabase) { fUseDatabase = useDatabase;}
644666df	68	Bool_t GetUseDatabase() const { return fUseDatabase;}
d2aa6df0	69
d2aa6df0	70	void SetResponseFunction(AliPID::EParticleType type, TObject * const o) { fResponseFunctions.AddAt(o,(Int_t)type); }
deae51a8	71	const TObject * GetResponseFunction(AliPID::EParticleType type) { return fResponseFunctions.At((Int_t)type); }
644666df	72	void SetVoltage(Int_t n, Float_t v) {fVoltageMap[n]=v;}
	73	void SetVoltageMap(const TVectorF& a) {fVoltageMap=a;} //resets ownership, ~ will not delete contents
	74	Float_t GetVoltage(Int_t n) const {return fVoltageMap[n];}
	75	void SetLowGainIROCthreshold(Float_t v) {fLowGainIROCthreshold=v;}
	76	void SetBadIROCthreshold(Float_t v) {fBadIROCthreshhold=v;}
	77	void SetLowGainOROCthreshold(Float_t v) {fLowGainOROCthreshold=v;}
	78	void SetBadOROCthreshold(Float_t v) {fBadOROCthreshhold=v;}
	79	void SetMaxBadLengthFraction(Float_t f) {fMaxBadLengthFraction=f;}
	80
	81	void SetMagField(Double_t mf) { fMagField=mf; }
d2aa6df0	82
644666df	83	//NEW
	84	void SetSigma(Float_t res0, Float_t resN2, ETPCgainScenario gainScenario );
	85	Double_t GetExpectedSignal( Double_t momentum,
	86	AliPID::EParticleType species,
	87	const TSpline3* responseFunction ) const;
	88	Double_t GetExpectedSignal( const AliVTrack* track,
	89	AliPID::EParticleType species,
	90	ETPCdEdxSource dedxSource );
	91	Double_t GetExpectedSigma( const AliVTrack* track,
	92	AliPID::EParticleType species,
	93	ETPCdEdxSource dedxSource );
	94	Double_t GetExpectedSigma( Double_t mom,
	95	Int_t nPoints,
	96	AliPID::EParticleType species,
	97	ETPCgainScenario gainScenario,
	98	const TSpline3* responseFunction) const;
	99	Float_t GetNumberOfSigmas( const AliVTrack* track,
	100	AliPID::EParticleType species,
	101	ETPCdEdxSource dedxSource );
	102
	103	void SetResponseFunction(TObject* o,
	104	AliPID::EParticleType type,
	105	ETPCgainScenario gainScenario);
	106	void Print(Option_t* option="") const;
	107	TSpline3* GetResponseFunction( AliPID::EParticleType species,
	108	ETPCgainScenario gainScenario ) const;
	109	TSpline3* GetResponseFunction( const AliVTrack* track,
	110	AliPID::EParticleType species,
	111	ETPCdEdxSource dedxSource );
	112	Bool_t ResponseFunctiondEdxN(const AliVTrack* track,
	113	AliPID::EParticleType species,
	114	ETPCdEdxSource dedxSource);
	115	Bool_t sectorNumbersInOut(const AliVTrack* track,
	116	Double_t innerRadius, Double_t outerRadius,
	117	Float_t& phiIn, Float_t& phiOut,
	118	Int_t& in, Int_t& out ) const;
	119	AliTPCPIDResponse::EChamberStatus TrackStatus(const AliVTrack* track, Int_t layer) const;
	120	Float_t MaxClusterRadius(const AliVTrack* track) const;
	121	Bool_t TrackApex(const AliVTrack* track, Float_t magField, Double_t position[3]) const;
	122	static const char* GainScenarioName(Int_t n) {return fgkGainScenarioName[(n>fgkNumberOfdEdxSourceScenarios)?fgkNumberOfdEdxSourceScenarios+1:n];}
	123	Int_t ResponseFunctionIndex( AliPID::EParticleType species,
	124	ETPCgainScenario gainScenario ) const;
	125	void ResetSplines();
	126
	127	void InvalidateCurrentValues();
	128	TSpline3* GetCurrentResponseFunction() const {return fCurrentResponseFunction;}
	129	Double_t GetCurrentdEdx() const {return fCurrentdEdx;}
	130	Int_t GetCurrentNPoints() const {return fCurrentNPoints;}
	131	ETPCgainScenario GetCurrentGainScenario() const {return fCurrentGainScenario;}
	132
	133	//OLD
10d100d4	134	Double_t GetExpectedSignal(const Float_t mom,
aea7a46d	135	AliPID::EParticleType n=AliPID::kKaon) const;
10d100d4	136	Double_t GetExpectedSigma(const Float_t mom, const Int_t nPoints,
644666df	137	AliPID::EParticleType n=AliPID::kKaon) const;
	138	Float_t GetNumberOfSigmas(const Float_t mom,
	139	const Float_t dEdx,
	140	const Int_t nPoints,
	141	AliPID::EParticleType n=AliPID::kKaon) const {
10d100d4	142
	143	Double_t bethe=GetExpectedSignal(mom,n);
	144	Double_t sigma=GetExpectedSigma(mom,nPoints,n);
	145	return (dEdx-bethe)/sigma;
	146	}
aea7a46d	147
10d100d4	148	Double_t GetMIP() const { return fMIP;}
644666df	149	Float_t GetRes0() const { return fRes0[0]; }
	150	Float_t GetResN2() const { return fResN2[0]; }
	151	Float_t GetRes0(ETPCgainScenario s) const { return fRes0[s]; }
	152	Float_t GetResN2(ETPCgainScenario s) const { return fResN2[s]; }
b981edcd	153
8c6a71ab	154	private:
10d100d4	155	Float_t fMIP; // dEdx for MIP
644666df	156	Float_t fRes0[fgkNumberOfGainScenarios]; // relative dEdx resolution rel sigma = fRes0*sqrt(1+fResN2/npoint)
644666df	157	Float_t fResN2[fgkNumberOfGainScenarios]; // relative Npoint dependence rel sigma = fRes0*sqrt(1+fResN2/npoint)
aea7a46d	158
	159	Double_t fKp1; // Parameters
	160	Double_t fKp2; // of
	161	Double_t fKp3; // the ALEPH
	162	Double_t fKp4; // Bethe-Bloch
	163	Double_t fKp5; // formula
	164
d2aa6df0	165	Bool_t fUseDatabase; // flag if fine-tuned database-response or simple ALEPH BB should be used
644666df	166
d2aa6df0	167	TObjArray fResponseFunctions; //! ObjArray of response functions individually for each particle
644666df	168	TVectorF fVoltageMap; //!stores a map of voltages wrt nominal for all chambers
	169	Float_t fLowGainIROCthreshold; //voltage threshold below which the IROC is considered low gain
	170	Float_t fBadIROCthreshhold; //voltage threshold for bad IROCS
	171	Float_t fLowGainOROCthreshold; //voltage threshold below which the OROC is considered low gain
	172	Float_t fBadOROCthreshhold; //voltage threshold for bad OROCS
	173	Float_t fMaxBadLengthFraction; //the maximum allowed fraction of track length in a bad sector.
	174
	175	TSpline3* fCurrentResponseFunction; //!response function for current track
	176	Double_t fCurrentdEdx; //!dEdx for currently processed track
	177	Int_t fCurrentNPoints; //!number of points used for dEdx calculation for current track
	178	ETPCgainScenario fCurrentGainScenario; //!gain scenario used for current track
	179	Int_t sectorNumber(Double_t phi) const;
	180
	181	Double_t fMagField; //! Magnetic field
	182
	183	static const char* fgkGainScenarioName[fgkNumberOfGainScenarios+1];
d2aa6df0	184
644666df	185	ClassDef(AliTPCPIDResponse,4) // TPC PID class
8c6a71ab	186	};
	187
	188	#endif
	189
	190