ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ECCEmRICHFastPIDMap.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ECCEmRICHFastPIDMap.cc
1 // $Id: $
2 
11 #include "ECCEmRICHFastPIDMap.h"
12 
13 #include <TF1.h>
14 #include <TFile.h>
15 #include <TH2F.h>
16 #include <TMath.h>
17 #include <TRandom.h>
18 #include <TVector3.h>
19 
20 #include <cassert>
21 #include <cmath>
22 #include <iostream>
23 
24 ECCEmRICHFastPIDMap::ECCEmRICHFastPIDMap(double trackResolution,
25  double incidentAngle, double pixS) {
26  fTrackResolution = trackResolution;
27  pLow = 0.6;
28  pHigh = 20;
29  c = 0.0299792458; // cm/picosecond
30  n = 1.03; // Aerogel
31  a = pixS; // pixel size 3.0; // mm -- one side
32  f = 152.4; // focal length mm = 6"
33  N_gam = 10;
34  mElectron = 0.00051099895; //GeV/c^2
35  mPion = 0.13957018; // GeV/c^2
36  mKaon = 0.493677; // GeV/c^2
37  mProton = 0.93827208816; // GeV/c^2
38  pi = 3.14159;
39  alpha = 0.0072973525693; // hyperfine const
40  L = 3.0; // Aerogel block thickness in cm
41 
42  //===============
43  th0 = incidentAngle; // incidence angle in radians
44 }
45 
46 ECCEmRICHFastPIDMap::~ECCEmRICHFastPIDMap() {}
47 
48 ECCEmRICHFastPIDMap::PIDCandidate_LogLikelihood_map
49 ECCEmRICHFastPIDMap::getFastSmearLogLikelihood(int truth_pid,
50  const double momentum,
51  const double theta_rad) const {
52  const int abs_truth_pid = abs(truth_pid);
53 
54  PIDCandidate_LogLikelihood_map ll_map;
55 
56  if (abs_truth_pid != EICPIDDefs::ElectronCandiate and
57  abs_truth_pid != EICPIDDefs::PionCandiate and
58  abs_truth_pid != EICPIDDefs::KaonCandiate and
59  abs_truth_pid != EICPIDDefs::ProtonCandiate) {
60  // not processing non-hadronic tracks
61  if (Verbosity())
62  std::cout << __PRETTY_FUNCTION__
63  << ": not processing non-hadronic tracks " << truth_pid
64  << std::endl;
65 
66  return ll_map;
67  }
68  if (theta_rad < m_acceptanceThetaMin or theta_rad > m_acceptanceThetaMax) {
69  // not processing out of acceptance tracks
70  if (Verbosity())
71  std::cout << __PRETTY_FUNCTION__
72  << ": not processing out of acceptance tracks theta_rad = "
73  << theta_rad << std::endl;
74 
75  return ll_map;
76  }
77  if (momentum < pLow or momentum > pHigh) {
78  // not processing out of acceptance tracks
79  if (Verbosity())
80  std::cout << __PRETTY_FUNCTION__
81  << ": not processing out of acceptance tracks momentum = "
82  << momentum << std::endl;
83 
84  return ll_map;
85  }
86 
87  double Nsigma_epi = 0;
88  double Nsigma_piK = 0;
89  double Nsigma_Kp = 0;
90 
91  //electron-pion case
92  if(momentum>0.59){ // pion cerenkov's threhsold at n=1.03 is ~0.59 GeV/c
93  // Angle difference
94  double dth = getAng(mElectron, momentum) - getAng(mPion, momentum);
95  // Detector uncertainty
96  double sigTh = sqrt(pow(getdAng(mPion, momentum), 2) +
97  pow(getdAng(mElectron, momentum), 2));
98  // Global uncertainty
99  double sigThTrk = getdAngTrk(mElectron, momentum);
100  double sigThc = sqrt(pow(sigTh / sqrt(getNgamma(L, mElectron, momentum)), 2) +
101  pow(sigThTrk, 2));
102  Nsigma_epi = dth / sigThc;
103  if (isnan(Nsigma_epi)) Nsigma_epi = 0;
104  }
105  //pion-Kaon case
106  if(momentum>2.0){ // kaon cerenkov's threhsold at n=1.03 is ~2 GeV/c
107  // Angle difference
108  double dth = getAng(mPion, momentum) - getAng(mKaon, momentum);
109  // Detector uncertainty
110  double sigTh = sqrt(pow(getdAng(mPion, momentum), 2) +
111  pow(getdAng(mKaon, momentum), 2));
112  // Global uncertainty
113  double sigThTrk = getdAngTrk(mPion, momentum);
114  double sigThc = sqrt(pow(sigTh / sqrt(getNgamma(L, mPion, momentum)), 2) +
115  pow(sigThTrk, 2));
116  Nsigma_piK = dth / sigThc;
117  if (isnan(Nsigma_piK)) Nsigma_piK = 0;
118  }
119  //kaon-proton case
120  if(momentum>3.8){ // Proton cerenkov's threhsold at n=1.03 is ~3.8 GeV/c
121  // Angle difference
122  double dth = getAng(mKaon, momentum) - getAng(mProton, momentum);
123  // Detector uncertainty
124  double sigTh = sqrt(pow(getdAng(mKaon, momentum), 2) +
125  pow(getdAng(mProton, momentum), 2));
126  // Global uncertainty
127  double sigThTrk = getdAngTrk(mKaon, momentum);
128  double sigThc = sqrt(pow(sigTh / sqrt(getNgamma(L, mKaon, momentum)), 2) +
129  pow(sigThTrk, 2));
130  Nsigma_Kp = dth / sigThc;
131  if (isnan(Nsigma_Kp)) Nsigma_Kp = 0;
132  }
133 
134  if (Verbosity())
135  std::cout << __PRETTY_FUNCTION__
136  << ": processing tracks momentum = " << momentum
137  << " Nsigma_epi = " << Nsigma_epi
138  << " Nsigma_piK = " << Nsigma_piK
139  << " Nsigma_Kp = " << Nsigma_Kp
140  << std::endl;
141 
142  const double electron_sigma_space_ring_radius = +Nsigma_piK + Nsigma_epi;
143  const double pion_sigma_space_ring_radius = +Nsigma_piK;
144  const double kaon_sigma_space_ring_radius = 0;
145  const double proton_sigma_space_ring_radius = -Nsigma_Kp;
146 
147  double sigma_space_ring_radius = gRandom->Gaus(0, 1);
148  if (abs_truth_pid == EICPIDDefs::ElectronCandiate)
149  sigma_space_ring_radius += electron_sigma_space_ring_radius;
150  else if (abs_truth_pid == EICPIDDefs::PionCandiate)
151  sigma_space_ring_radius += pion_sigma_space_ring_radius;
152  else if (abs_truth_pid == EICPIDDefs::KaonCandiate)
153  sigma_space_ring_radius += kaon_sigma_space_ring_radius;
154  else if (abs_truth_pid == EICPIDDefs::ProtonCandiate)
155  sigma_space_ring_radius += proton_sigma_space_ring_radius;
156 
157  ll_map[EICPIDDefs::ElectronCandiate] =
158  -0.5 * pow(sigma_space_ring_radius - electron_sigma_space_ring_radius, 2);
159  ll_map[EICPIDDefs::PionCandiate] =
160  -0.5 * pow(sigma_space_ring_radius - pion_sigma_space_ring_radius, 2);
161  ll_map[EICPIDDefs::KaonCandiate] =
162  -0.5 * pow(sigma_space_ring_radius - kaon_sigma_space_ring_radius, 2);
163  ll_map[EICPIDDefs::ProtonCandiate] =
164  -0.5 * pow(sigma_space_ring_radius - proton_sigma_space_ring_radius, 2);
165 
166  return ll_map;
167 }
168 
169 // Angle exiting the Aerogel
170 double ECCEmRICHFastPIDMap::getAng(double mass, double mom) const {
171  double beta = mom / sqrt(pow(mom, 2) + pow(mass, 2));
172  double thc = acos(1. / n / beta);
173  double th0p = asin(sin(th0) / n);
174  double dthc = thc - th0p;
175  double theta = asin(n * sin(dthc));
176 
177  return theta;
178 }
179 
180 // Uncertainty due to detector effects
181 double ECCEmRICHFastPIDMap::getdAng(double mass, double mom) const {
182  double beta = mom / sqrt(pow(mom, 2) + pow(mass, 2));
183  double thc = acos(1. / n / beta);
184  double th0p = asin(sin(th0) / n);
185  double dthc = thc - th0p;
186  double theta = asin(n * sin(dthc));
187 
188  double sig_ep = 0; // Emission point error
189  double sig_chro = 0; // Chromatic dispersion error
190  double sig_pix = a * pow(cos(theta), 2) / f / sqrt(12.);
191 
192  double sigTh = sqrt(pow(sig_ep, 2) + pow(sig_chro, 2) + pow(sig_pix, 2));
193 
194  return sigTh;
195 }
196 
197 // Uncertainty due to tracking resolution
198 double ECCEmRICHFastPIDMap::getdAngTrk(double mass, double mom) const {
199  double beta = mom / sqrt(pow(mom, 2) + pow(mass, 2));
200  double thc = acos(1. / n / beta);
201  double th0p = asin(sin(th0) / n);
202  double dthc = thc - th0p;
203  double theta = asin(n * sin(dthc));
204 
205  double sig_trk =
206  (cos(dthc) / cos(theta)) * (cos(th0) / cos(th0p)) * fTrackResolution;
207 
208  return sig_trk;
209 }
210 
211 // no. of gamms
212 double ECCEmRICHFastPIDMap::getNgamma(double t, double mass, double mom) const {
213  int tot = 10000;
214  double beta = mom / sqrt(pow(mom, 2) + pow(mass, 2));
215  double fact = 2. * pi * alpha * t * (1. - 1. / pow(n * beta, 2));
216  double T_lensWin = 0.92 * 0.92;
217  double xmin = 300.e-7;
218  double xmax = 650.e-7;
219  double dx = (xmax - xmin) / tot;
220  double sum = 0;
221  for (int j = 0; j < tot; j++) {
222  double x = xmin + j * dx + dx / 2.;
223  sum += T_QE(x) * T_Aer(t, x) / pow(x, 2);
224  }
225  return fact * T_lensWin * sum * dx;
226 }
227 
228 // Quantum efficiency
229 double ECCEmRICHFastPIDMap::T_QE(double lam) const {
230  return 0.34 * exp(-1. * pow(lam - 345.e-7, 2) / (2. * pow(119.e-7, 2)));
231 }
232 
233 // Transmissions of the radiator block
234 double ECCEmRICHFastPIDMap::T_Aer(double t, double lam) const {
235  return 0.83 * exp(-1. * t * 56.29e-20 / pow(lam, 4));
236 }