ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
genericRICH.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file genericRICH.h
1 
2 
3 
4 
5 #ifndef __GENERICRICH_H__
6 #define __GENERICRICH_H__
7 
8 #include "genericDetector.h"
9 
10 class genericRICH : public genericDetector
11 {
12  public:
13  genericRICH() = default;
14  virtual ~genericRICH() = default;
15 
17  void setIndex(double val) { mIndex = val; };
18  void setEfficiency(double val) { mEfficiency = val; };
19  void setMinPhotons(double val) { mMinPhotons = val; };
20  void setThresholdMode(bool val) { mThresholdMode = val; };
21 
22  void setChromaticSigma(int n, double *valx, double *valy) {
23  if (mChromaticSigma) delete mChromaticSigma;
24  mChromaticSigma = new TGraph(n, valx, valy);
25  }
26  void setPositionSigma(int n, double *valx, double *valy) {
27  if (mPositionSigma) delete mPositionSigma;
28  mPositionSigma = new TGraph(n, valx, valy);
29  }
30  void setEmissionSigma(int n, double *valx, double *valy) {
31  if (mEmissionSigma) delete mEmissionSigma;
32  mEmissionSigma = new TGraph(n, valx, valy);
33  }
34  void setFieldSigma(int n, double *valx, double *valy) {
35  if (mFieldSigma) delete mFieldSigma;
36  mFieldSigma = new TGraph(n, valx, valy);
37  }
38  void setTrackingSigma(int n, double *valx, double *valy) {
39  if (mTrackingSigma) delete mTrackingSigma;
40  mTrackingSigma = new TGraph(n, valx, valy);
41  }
42 
44  double numSigma (double eta, double p, PID::type PID) override;
45  double maxP (double eta, double nsigma, PID::type PID) override;
46  double minP (double eta, double nsigma, PID::type PID) override;
47 
48  double cherenkovAngle(double p, double m) const { return acos( sqrt( m * m + p * p ) / ( mIndex * p ) ); };
49  double cherenkovThreshold(double m) const { return m / sqrt(mIndex * mIndex - 1.); };
50  double numberOfPhotons(double angle) const { return 490. * sin(angle) * sin(angle) * mLength; };
51  double numberOfDetectedPhotons(double angle) const { return numberOfPhotons(angle) * mEfficiency; };
52  double cherenkovAngleSigma(double eta, double p, double m) const;
53 
54  protected:
55 
56  // RICH parameters
57  double mIndex = 1.0014; // refractive index
58  double mEfficiency = 0.25; // overall photon detection efficiency
59  double mMinPhotons = 3.; // minimum number of detected photons
60 
61  // contributions to resolution
62  TGraph *mChromaticSigma = nullptr; // chromatic resolution vs. polar angle [rad]
63  TGraph *mPositionSigma = nullptr; // position resolution vs. polar angle [rad]
64  TGraph *mEmissionSigma = nullptr; // emission resolution vs. polar angle [rad]
65  TGraph *mFieldSigma = nullptr; // field resolution vs. polar angle [rad]
66  TGraph *mTrackingSigma = nullptr; // tracking resolution vs. polar angle [rad]
67 
68  // threshold mode
69  bool mThresholdMode = true;
70 
71 };
72 
73 double genericRICH::cherenkovAngleSigma(double eta, double p, double m) const
74 {
75 
76  auto theta = 2. * atan(exp(-eta)) * 57.295780;
77  auto chromatic = mChromaticSigma ? mChromaticSigma->Eval(theta) : 0.;
78  auto position = mPositionSigma ? mPositionSigma->Eval(theta) : 0.;
79  auto emission = mEmissionSigma ? mEmissionSigma->Eval(theta) : 0.;
80  auto field = mFieldSigma ? mFieldSigma->Eval(theta) : 0.;
81  auto tracking = mTrackingSigma ? mTrackingSigma->Eval(theta) : 0.;
82 
83  // contributions that scale with number of detected photons
84  auto ndet = numberOfDetectedPhotons(cherenkovAngle(p, m));
85  auto sigma1 = sqrt(chromatic * chromatic +
86  position * position +
87  emission * emission +
88  field * field +
89  tracking * tracking);
90  // contributions that do not
91  auto sigma2 = 0.;
92  //
93  return sqrt(sigma1 * sigma1 / ndet + sigma2 * sigma2);
94 };
95 
96 
97 double genericRICH::numSigma(double eta, double p, PID::type PID)
98 {
99  double mass1, mass2;
100  switch (PID) {
101  case pi_k:
102  mass1 = mMassPion;
103  mass2 = mMassKaon;
104  break;
105  case k_p:
106  mass1 = mMassKaon;
107  mass2 = mMassProton;
108  break;
109  }
110 
111  double thr1 = cherenkovThreshold(mass1);
112  double thr2 = cherenkovThreshold(mass2);
113 
115  if (p > thr1 && p > thr2)
116  return (cherenkovAngle(p, mass1) - cherenkovAngle(p, mass2)) / cherenkovAngleSigma(eta, p, mass1);
117 
119  if (mThresholdMode && p > thr1)
120  return (cherenkovAngle(thr2 + 0.001, mass1) - cherenkovAngle(thr2 + 0.001, mass2)) / cherenkovAngleSigma(eta, thr2 + 0.001, mass1);
121 
123  return 0.;
124 }
125 
126 double genericRICH::maxP(double eta, double nsigma, PID::type PID)
127 {
128  double mass1, mass2;
129  switch (PID) {
130  case pi_k:
131  mass1 = mMassPion;
132  mass2 = mMassKaon;
133  break;
134  case k_p:
135  mass1 = mMassKaon;
136  mass2 = mMassProton;
137  break;
138  }
139 
141  double p = minP(eta, nsigma, PID) + 0.001;
142  while (numSigma(eta, p, PID) > nsigma) p += 0.001;
143  return p;
144 }
145 
146 double genericRICH::minP(double eta, double nsigma, PID::type PID)
147 {
148  double mass;
149  switch (PID) {
150  case pi_k:
151  mass = mThresholdMode ? mMassPion : mMassKaon;
152  break;
153  case k_p:
154  mass = mThresholdMode ? mMassKaon : mMassProton;
155  break;
156  }
157 
159  double p = cherenkovThreshold(mass);
160  while (numberOfPhotons(cherenkovAngle(p, mass)) * mEfficiency < mMinPhotons) p += 0.001;
161  return std::max(p, pMin(eta));
162 }
163 
164 #endif /* __GENERICRICH_H__ */