ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
mRICH.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file mRICH.C
1 #include "mRICH.h"
2 
3 using namespace std;
4 
5 
6 mRICH::mRICH(double trackResolution, double timePrecision, double pix, double p)
7 {
8 
9  // In this arameterization we assume that the particle enters the mRICH perpendicular to its front face. More realistic case will be implemented in the subsequent updates.
10 
11 
12  fTrackResolution = trackResolution;
13  fTimePrecision = timePrecision;
14  mom = p;
15  pLow = 3;
16  pHigh = 10.1;
17  c = 0.0299792458; // cm/picosecond
18  n = 1.03; //Aerogel
19  a = pix;// pixel size 3.0; // mm -- one side
20  f = 152.4; //focal length mm = 6"
21  N_gam = 10;
22  mPion = 0.13957018; //GeV/c^2
23  mKaon = 0.493677; //GeV/c^2
24  mProton = 0.93827208816; //GeV/c^2
25  pi = 3.14159;
26  alpha = 0.0072973525693; // hyperfine const
27  L = 3.0;//Aerogel block thickness in cm
28 
29  //===============
30  th0 = 0.; //incidence angle in radians
31 
32  //Angle difference
33  double dth = getAng(mPion) - getAng(mKaon);
34  //Detector uncertainty
35  double sigTh = sqrt(pow(getdAng(mPion),2)+pow(getdAng(mKaon),2));
36  //Global uncertainty
37  double sigThTrk = sqrt(pow(getdAngTrk(mPion),2)+pow(getdAngTrk(mKaon),2));
38  double sigThc = sqrt(pow(sigTh/sqrt(getNgamma(L,mKaon)),2)+pow(sigThTrk,2));
39  Nsigma = dth/sigThc;
40 
41 #if 0
42  //From simulations
43  ReadMap("mRICHresol.root");
44  Nsigma = fpixMap->GetBinContent(p,pix);
45 #endif
46 
47  //cout<<getNgamma(mPion)<<"\t"<<getNgamma(mKaon)<<endl;
48 
49  char nameit[500];
50  sprintf(nameit,"mRICH Pixel Size =%dx%d (mm^{2}) / p =%.f (GeV/c)",pix,pix,p);
51  myName = nameit;
52 
53  //cout<<".................\t"<<Nsigma<<endl;
54 }
55 
56 double mRICH::maxP(double eta, double p, double numSigma, PID::type PID)
57 {
58  if (valid(eta,p))
59  {
60  return Nsigma;
61  }
62  else
63  {
64  cout << "mRICH.C: Invalid (eta) for this detector." <<endl;
65  return 0.0;
66  }
67 }
68 
69 //Angle exiting the Aerogel
70 double mRICH::getAng(double mass)
71 {
72 
73  double beta = mom/sqrt(pow(mom,2)+pow(mass,2));
74  double thc = acos(1./n/beta);
75  double th0p = asin(sin(th0)/n);
76  double dthc = thc - th0p;
77  double theta = asin(n*sin(dthc));
78 
79  return theta;
80 }
81 
82 //Uncertainty due to detector effects
83 double mRICH::getdAng(double mass)
84 {
85 
86  double beta = mom/sqrt(pow(mom,2)+pow(mass,2));
87  double thc = acos(1./n/beta);
88  double th0p = asin(sin(th0)/n);
89  double dthc = thc - th0p;
90  double theta = asin(n*sin(dthc));
91 
92  double sig_ep = 0;//Emission point error
93  double sig_chro = 0; //Chromatic dispersion error
94  double sig_pix = a*pow(cos(theta),2)/f/sqrt(12.);;
95 
96  double sigTh = sqrt(pow(sig_ep,2)+pow(sig_chro,2)+pow(sig_pix,2));
97 
98  return sigTh;
99 }
100 
101 //Uncertainty due to tracking resolution
102 double mRICH::getdAngTrk(double mass)
103 {
104 
105  double beta = mom/sqrt(pow(mom,2)+pow(mass,2));
106  double thc = acos(1./n/beta);
107  double th0p = asin(sin(th0)/n);
108  double dthc = thc - th0p;
109  double theta = asin(n*sin(dthc));
110 
111  double sig_trk = (cos(dthc)/cos(theta))*(cos(th0)/cos(th0p))*fTrackResolution;;
112 
113  return sig_trk;
114 }
115 
116 //no. of gamms
117 double mRICH::getNgamma(double t, double mass)
118 {
119  int tot = 10000;
120  double beta = mom/sqrt(pow(mom,2)+pow(mass,2));
121  double fact = 2.*pi*alpha*t*(1.-1./pow(n*beta,2));
122  double T_lensWin = 0.92*0.92;
123  double xmin = 300.e-7;
124  double xmax = 650.e-7;
125  double dx = (xmax-xmin)/tot;
126  double sum = 0;
127  for(int j=0; j<tot; j++){
128  double x = xmin + j*dx + dx/2.;
129  sum += T_QE(x)*T_Aer(t,x)/pow(x,2);
130  }
131  return fact*T_lensWin*sum*dx;
132 }
133 
134 //Quantum efficiency
135 double mRICH::T_QE(double lam)
136 {
137  return 0.34*exp(-1.*pow(lam-345.e-7,2)/(2.*pow(119.e-7,2)));
138 }
139 
140 //Transmissions of the radiator block
141 double mRICH::T_Aer(double t, double lam)
142 {
143  return 0.83*exp(-1.*t*56.29e-20/pow(lam,4));
144 }
145 
146 
147 //Read histogram from MC parametrization
148 void mRICH::ReadMap(TString name){
149  TFile* file = TFile::Open(name);
150  fpixMap = (TH2F *) file->Get("h2pix");
151 }
152 
153 
155 {
156  // Here one should describe what this detector is and what assumptions are
157  // made in calculating the performance specifications.
158 
159  cout << "My name is \"" << myName << "\" and I am described as follows:" <<endl;
160  cout << " Momentum coverage for k/pi separation= [" << 3.0 << " to " << 10.0 <<" (GeV/c) ]"<<endl;
161  cout << " Since mRICH detector system consists of array of identical shoebox-like modules, in principle, " << endl;
162  cout << " the mRICH pid performance is invariant relative to the tracking polar angle but the momentum of the " << endl;
163  cout << " charge particle and the entrance point location on the front face of a given mRICH module." << endl;
164  cout << " Assumed time precision = " << fTimePrecision << " ns" <<endl;
165  cout << " Assumed track resolution = "<< fTrackResolution << " mrad" <<endl;
166  //cout << " Assumed quantum efficiency of the photosensors = "<< fSensorQMefficiency << "%" <<endl;
167  cout << endl;
168 
169 }