CMS 3D CMS Logo

HFCherenkov.cc
Go to the documentation of this file.
1 // File: HFCherenkov.cc
3 // Description: Generate Cherenkov photons for HF
5 
7 
8 #include "G4Poisson.hh"
9 #include "G4ParticleDefinition.hh"
10 #include "G4NavigationHistory.hh"
11 #include "TMath.h"
12 
13 #include "Randomize.hh"
14 #include "G4SystemOfUnits.hh"
15 
16 //#define EDM_ML_DEBUG
17 
19  ref_index = m_HF.getParameter<double>("RefIndex");
20  lambda1 = ((m_HF.getParameter<double>("Lambda1")) / 1.e+7) * CLHEP::cm;
21  lambda2 = ((m_HF.getParameter<double>("Lambda2")) / 1.e+7) * CLHEP::cm;
22  aperture = m_HF.getParameter<double>("Aperture");
23  gain = m_HF.getParameter<double>("Gain");
24  checkSurvive = m_HF.getParameter<bool>("CheckSurvive");
25  UseNewPMT = m_HF.getParameter<bool>("UseR7600UPMT");
26  fibreR = m_HF.getParameter<double>("FibreR") * CLHEP::mm;
27 
29  sinPsimax = 1 / ref_index;
30 
31  edm::LogVerbatim("HFShower") << "HFCherenkov:: initialised with ref_index " << ref_index << " lambda1/lambda2 (cm) "
32  << lambda1 / CLHEP::cm << "|" << lambda2 / CLHEP::cm << " aperture(trapped) " << aperture
33  << "|" << aperturetrapped << " sinPsimax " << sinPsimax
34  << " Check photon survival in HF " << checkSurvive << " Gain " << gain << " useNewPMT "
35  << UseNewPMT << " FibreR " << fibreR;
36 
37  clearVectors();
38 }
39 
41 
42 int HFCherenkov::computeNPE(const G4Step* aStep,
43  const G4ParticleDefinition* pDef,
44  double pBeta,
45  double u,
46  double v,
47  double w,
48  double step_length,
49  double zFiber,
50  double dose,
51  int npe_Dose) {
52  clearVectors();
53  if (!isApplicable(pDef)) {
54  return 0;
55  }
56  if (pBeta < (1 / ref_index) || step_length < 0.0001) {
57 #ifdef EDM_ML_DEBUG
58  edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: pBeta " << pBeta << " 1/mu " << (1 / ref_index)
59  << " step_length " << step_length;
60 #endif
61  return 0;
62  }
63 
64  double uv = std::sqrt(u * u + v * v);
65  int nbOfPhotons = computeNbOfPhotons(pBeta, step_length) * aStep->GetTrack()->GetWeight();
66 #ifdef EDM_ML_DEBUG
67  edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: pBeta " << pBeta << " u/v/w " << u << "/" << v << "/" << w
68  << " step_length " << step_length << " zFib " << zFiber << " nbOfPhotons "
69  << nbOfPhotons;
70 #endif
71  if (nbOfPhotons < 0) {
72  return 0;
73  } else if (nbOfPhotons > 0) {
74  G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
75  const G4TouchableHandle& theTouchable = preStepPoint->GetTouchableHandle();
76  G4ThreeVector localprepos =
77  theTouchable->GetHistory()->GetTopTransform().TransformPoint(aStep->GetPreStepPoint()->GetPosition());
78  G4ThreeVector localpostpos =
79  theTouchable->GetHistory()->GetTopTransform().TransformPoint(aStep->GetPostStepPoint()->GetPosition());
80 
81  double length = std::sqrt((localpostpos.x() - localprepos.x()) * (localpostpos.x() - localprepos.x()) +
82  (localpostpos.y() - localprepos.y()) * (localpostpos.y() - localprepos.y()));
83  double yemit = std::sqrt(fibreR * fibreR - length * length / 4.);
84 
85  double u_ph = 0, v_ph = 0, w_ph = 0;
86  for (int i = 0; i < nbOfPhotons; i++) {
87  double rand = G4UniformRand();
88  double theta_C = acos(1. / (pBeta * ref_index));
89  double phi_C = 2 * M_PI * rand;
90  double sinTheta = sin(theta_C);
91  double cosTheta = cos(theta_C);
92  double cosPhi = cos(phi_C);
93  double sinPhi = sin(phi_C);
94  //photon momentum
95  if (uv < 0.001) { // aligned with z-axis
96  u_ph = sinTheta * cosPhi;
97  v_ph = sinTheta * sinPhi;
98  w_ph = cosTheta;
99  } else { // general case
100  u_ph = uv * cosTheta + sinTheta * cosPhi * w;
101  v_ph = sinTheta * sinPhi;
102  w_ph = w * cosTheta - sinTheta * cosPhi * uv;
103  }
104  double r_lambda = G4UniformRand();
105  double lambda0 = (lambda1 * lambda2) / (lambda2 - r_lambda * (lambda2 - lambda1));
106  double lambda = (lambda0 / CLHEP::cm) * 1.e+7; // lambda is in nm
107  wlini.push_back(lambda);
108 #ifdef EDM_ML_DEBUG
109  edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: " << i << " lambda " << lambda << " w_ph " << w_ph;
110 #endif
111  // --------------
112  double xemit = length * (G4UniformRand() - 0.5);
113  double gam = atan2(yemit, xemit);
114  double eps = atan2(v_ph, u_ph);
115  double sinBeta = sin(gam - eps);
116  double rho = std::sqrt(xemit * xemit + yemit * yemit);
117  double sinEta = rho / fibreR * sinBeta;
118  double cosEta = std::sqrt(1. - sinEta * sinEta);
119  double sinPsi = std::sqrt(1. - w_ph * w_ph);
120  double cosKsi = cosEta * sinPsi;
121 
122 #ifdef EDM_ML_DEBUG
123  if (cosKsi < aperturetrapped && w_ph > 0.) {
124  edm::LogVerbatim("HFShower") << "HFCherenkov::Trapped photon : " << u_ph << " " << v_ph << " " << w_ph << " "
125  << xemit << " " << gam << " " << eps << " " << sinBeta << " " << rho << " "
126  << sinEta << " " << cosEta << " "
127  << " " << sinPsi << " " << cosKsi;
128  } else {
129  edm::LogVerbatim("HFShower") << "HFCherenkov::Rejected photon : " << u_ph << " " << v_ph << " " << w_ph << " "
130  << xemit << " " << gam << " " << eps << " " << sinBeta << " " << rho << " "
131  << sinEta << " " << cosEta << " "
132  << " " << sinPsi << " " << cosKsi;
133  }
134 #endif
135  if (cosKsi < aperturetrapped // photon is trapped inside fiber
136  && w_ph > 0. // and moves to PMT
137  && sinPsi < sinPsimax) { // and is not reflected at fiber end
138  wltrap.push_back(lambda);
139  double prob_HF = 1.0;
140  if (checkSurvive) {
141  double a0_inv = 0.1234; //meter^-1
142  double a_inv = a0_inv + 0.14 * pow(dose, 0.30); // ?
143  double z_meters = zFiber;
144  prob_HF = exp(-z_meters * a_inv);
145  }
146  rand = G4UniformRand();
147 #ifdef EDM_ML_DEBUG
148  edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: probHF " << prob_HF << " Random " << rand
149  << " Survive? " << (rand < prob_HF);
150 #endif
151  if (rand < prob_HF) { // survived in HF
152  wlatten.push_back(lambda);
153  rand = G4UniformRand();
154  double effHEM = computeHEMEff(lambda);
155  // compute number of bounces in air guide
156  rand = G4UniformRand();
157  double phi = 0.5 * M_PI * rand;
158  double rad_bundle = 19.; // bundle radius
159  double rad_lg = 25.; // light guide radius
160  rand = G4UniformRand();
161  double rad = rad_bundle * std::sqrt(rand);
162  double rho = rad * sin(phi);
163  double tlength = 2. * std::sqrt(rad_lg * rad_lg - rho * rho);
164  double length_lg = 430;
165  double sin_air = sinPsi * 1.46;
166  double tang = sin_air / std::sqrt(1. - sin_air * sin_air);
167  int nbounce = length_lg / tlength * tang + 0.5;
168  double eff = pow(effHEM, nbounce);
169 #ifdef EDM_ML_DEBUG
170  edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: w_ph " << w_ph << " effHEM " << effHEM << " eff "
171  << eff << " Random " << rand << " Survive? " << (rand < eff);
172 #endif
173  if (rand < eff) { // survived HEM
174  wlhem.push_back(lambda);
175  double qEffic = computeQEff(lambda);
176  rand = G4UniformRand();
177 #ifdef EDM_ML_DEBUG
178  edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: qEffic " << qEffic << " Random " << rand
179  << " Survive? " << (rand < qEffic);
180 #endif
181  if (rand < qEffic) { // made photoelectron
182  npe_Dose += 1;
183  momZ.push_back(w_ph);
184  wl.push_back(lambda);
185  wlqeff.push_back(lambda);
186  } // made pe
187  } // passed HEM
188  } // passed fiber
189  } // end of if(w_ph < w_aperture), trapped inside fiber
190  } // end of ++NbOfPhotons
191  } // end of if(NbOfPhotons)}
192  int npe = npe_Dose; // Nb of photoelectrons
193 #ifdef EDM_ML_DEBUG
194  edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPE: npe " << npe;
195 #endif
196  return npe;
197 }
198 
200  const G4ParticleDefinition* pDef, double pBeta, double u, double v, double w, double step_length) {
201  clearVectors();
202  int npe_ = 0;
203  if (!isApplicable(pDef)) {
204  return 0;
205  }
206  if (pBeta < (1 / ref_index) || step_length < 0.0001) {
207 #ifdef EDM_ML_DEBUG
208  edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPEinPMT: pBeta " << pBeta << " 1/mu " << (1 / ref_index)
209  << " step_length " << step_length;
210 #endif
211  return 0;
212  }
213 
214  double uv = std::sqrt(u * u + v * v);
215  int nbOfPhotons = computeNbOfPhotons(pBeta, step_length);
216 #ifdef EDM_ML_DEBUG
217  edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPEinPMT: pBeta " << pBeta << " u/v/w " << u << "/" << v << "/"
218  << w << " step_length " << step_length << " nbOfPhotons " << nbOfPhotons;
219 #endif
220  if (nbOfPhotons < 0) {
221  return 0;
222  } else if (nbOfPhotons > 0) {
223  double w_ph = 0;
224  for (int i = 0; i < nbOfPhotons; i++) {
225  double rand = G4UniformRand();
226  double theta_C = acos(1. / (pBeta * ref_index));
227  double phi_C = 2 * M_PI * rand;
228  double sinTheta = sin(theta_C);
229  double cosTheta = cos(theta_C);
230  double cosPhi = cos(phi_C);
231  //photon momentum
232  if (uv < 0.001) { // aligned with z-axis
233  w_ph = cosTheta;
234  } else { // general case
235  w_ph = w * cosTheta - sinTheta * cosPhi * uv;
236  }
237  double r_lambda = G4UniformRand();
238  double lambda0 = (lambda1 * lambda2) / (lambda2 - r_lambda * (lambda2 - lambda1));
239  double lambda = (lambda0 / CLHEP::cm) * 1.e+7; // lambda is in nm
240  wlini.push_back(lambda);
241 #ifdef EDM_ML_DEBUG
242  edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPEinPMT: " << i << " lambda " << lambda << " w_ph " << w_ph
243  << " aperture " << aperture;
244 #endif
245  if (w_ph > aperture) { // phton trapped inside PMT glass
246  wltrap.push_back(lambda);
247  rand = G4UniformRand();
248 #ifdef EDM_ML_DEBUG
249  edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPEinPMT: Random " << rand << " Survive? " << (rand < 1.);
250 #endif
251  if (rand < 1.0) { // survived all the times and sent to photo-cathode
252  wlatten.push_back(lambda);
253  double qEffic = computeQEff(lambda); //Quantum efficiency of the PMT
254  rand = G4UniformRand();
255 #ifdef EDM_ML_DEBUG
256  edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPEinPMT: qEffic " << qEffic << " Random " << rand
257  << " Survive? " << (rand < qEffic);
258 #endif
259  if (rand < qEffic) { // made photoelectron
260  npe_ += 1;
261  momZ.push_back(w_ph);
262  wl.push_back(lambda);
263  wlqeff.push_back(lambda);
264  } // made pe
265  } // accepted all Cherenkov photons
266  } // end of if(w_ph < w_aperture), trapped inside glass
267  } // end of ++NbOfPhotons
268  } // end of if(NbOfPhotons)}
269 #ifdef EDM_ML_DEBUG
270  edm::LogVerbatim("HFShower") << "HFCherenkov::computeNPEinPMT: npe " << npe_;
271 #endif
272  return npe_;
273 }
274 
275 std::vector<double> HFCherenkov::getWLIni() {
276  std::vector<double> v = wlini;
277  return v;
278 }
279 
280 std::vector<double> HFCherenkov::getWLTrap() {
281  std::vector<double> v = wltrap;
282  return v;
283 }
284 
285 std::vector<double> HFCherenkov::getWLAtten() {
286  std::vector<double> v = wlatten;
287  return v;
288 }
289 
290 std::vector<double> HFCherenkov::getWLHEM() {
291  std::vector<double> v = wlhem;
292  return v;
293 }
294 
295 std::vector<double> HFCherenkov::getWLQEff() {
296  std::vector<double> v = wlqeff;
297  return v;
298 }
299 
300 std::vector<double> HFCherenkov::getWL() {
301  std::vector<double> v = wl;
302  return v;
303 }
304 
305 std::vector<double> HFCherenkov::getMom() {
306  std::vector<double> v = momZ;
307  return v;
308 }
309 
310 int HFCherenkov::computeNbOfPhotons(double beta, G4double stepL) {
311  double pBeta = beta;
312  double alpha = 0.0073;
313  double step_length = stepL;
314  double theta_C = acos(1. / (pBeta * ref_index));
315  double lambdaDiff = (1. / lambda1 - 1. / lambda2);
316  double cherenPhPerLength = 2 * M_PI * alpha * lambdaDiff * CLHEP::cm;
317  double d_NOfPhotons = cherenPhPerLength * sin(theta_C) * sin(theta_C) * (step_length / CLHEP::cm);
318  int nbOfPhotons = int(d_NOfPhotons);
319 #ifdef EDM_ML_DEBUG
320  edm::LogVerbatim("HFShower") << "HFCherenkov::computeNbOfPhotons: StepLength " << step_length << " theta_C "
321  << theta_C << " lambdaDiff " << lambdaDiff << " cherenPhPerLength " << cherenPhPerLength
322  << " Photons " << d_NOfPhotons << " " << nbOfPhotons;
323 #endif
324  return nbOfPhotons;
325 }
326 
327 double HFCherenkov::computeQEff(double wavelength) {
328  double qeff(0.);
329  if (UseNewPMT) {
330  if (wavelength <= 350) {
331  qeff = 2.45867 * (TMath::Landau(wavelength, 353.820, 59.1324));
332  } else if (wavelength > 350 && wavelength < 500) {
333  qeff = 0.441989 * exp(-pow((wavelength - 358.371), 2) / (2 * pow((138.277), 2)));
334  } else if (wavelength >= 500 && wavelength < 550) {
335  qeff = 0.271862 * exp(-pow((wavelength - 491.505), 2) / (2 * pow((47.0418), 2)));
336  } else if (wavelength >= 550) {
337  qeff = 0.137297 * exp(-pow((wavelength - 520.260), 2) / (2 * pow((75.5023), 2)));
338  }
339 #ifdef EDM_ML_DEBUG
340  edm::LogVerbatim("HFShower") << "HFCherenkov:: for new PMT : wavelength === " << wavelength << "\tqeff ===\t"
341  << qeff;
342 #endif
343  } else {
344  double y = (wavelength - 275.) / 180.;
345  double func = 1. / (1. + 250. * pow((y / 5.), 4));
346  double qE_R7525 = 0.77 * y * exp(-y) * func;
347  qeff = qE_R7525;
348 #ifdef EDM_ML_DEBUG
349  edm::LogVerbatim("HFShower") << "HFCherenkov::computeQEff: wavelength " << wavelength << " y/func " << y << "/"
350  << func << " qeff " << qeff;
351  edm::LogVerbatim("HFShower") << "HFCherenkov:: for old PMT : wavelength === " << wavelength << "; qeff = " << qeff;
352 #endif
353  }
354  return qeff;
355 }
356 
357 double HFCherenkov::computeHEMEff(double wavelength) {
358  double hEMEff = 0;
359  if (wavelength < 400.) {
360  hEMEff = 0.;
361  } else if (wavelength >= 400. && wavelength < 410.) {
362  //hEMEff = .99 * (wavelength - 400.) / 10.;
363  hEMEff = (-1322.453 / wavelength) + 4.2056;
364  } else if (wavelength >= 410.) {
365  hEMEff = 0.99;
366  if (wavelength > 415. && wavelength < 445.) {
367  //abs(wavelength - 430.) < 15.
368  //hEMEff = 0.95;
369  hEMEff = 0.97;
370  } else if (wavelength > 550. && wavelength < 600.) {
371  // abs(wavelength - 575.) < 25.)
372  //hEMEff = 0.96;
373  hEMEff = 0.98;
374  } else if (wavelength > 565. && wavelength <= 635.) { // added later
375  // abs(wavelength - 600.) < 35.)
376  hEMEff = (701.7268 / wavelength) - 0.186;
377  }
378  }
379 #ifdef EDM_ML_DEBUG
380  edm::LogVerbatim("HFShower") << "HFCherenkov::computeHEMEff: wavelength " << wavelength << " hEMEff " << hEMEff;
381 #endif
382  return hEMEff;
383 }
384 
385 double HFCherenkov::smearNPE(int npe) {
386  double pe = 0.;
387  if (npe > 0) {
388  for (int i = 0; i < npe; ++i) {
389  double val = G4Poisson(gain);
390  pe += (val / gain) + 0.001 * G4UniformRand();
391  }
392  }
393 #ifdef EDM_ML_DEBUG
394  edm::LogVerbatim("HFShower") << "HFCherenkov::smearNPE: npe " << npe << " pe " << pe;
395 #endif
396  return pe;
397 }
398 
400  wl.clear();
401  wlini.clear();
402  wltrap.clear();
403  wlatten.clear();
404  wlhem.clear();
405  wlqeff.clear();
406  momZ.clear();
407 }
408 
409 bool HFCherenkov::isApplicable(const G4ParticleDefinition* aParticleType) {
410  bool tmp = (aParticleType->GetPDGCharge() != 0);
411 #ifdef EDM_ML_DEBUG
412  edm::LogVerbatim("HFShower") << "HFCherenkov::isApplicable: aParticleType " << aParticleType->GetParticleName()
413  << " PDGCharge " << aParticleType->GetPDGCharge() << " Result " << tmp;
414 #endif
415  return tmp;
416 }
Log< level::Info, true > LogVerbatim
double computeQEff(double wavelength)
Definition: HFCherenkov.cc:327
float alpha
Definition: AMPTWrapper.h:105
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
virtual ~HFCherenkov()
Definition: HFCherenkov.cc:40
double fibreR
Definition: HFCherenkov.h:61
double aperturetrapped
Definition: HFCherenkov.h:60
T w() const
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
HFCherenkov(edm::ParameterSet const &p)
Definition: HFCherenkov.cc:18
std::vector< double > wlqeff
Definition: HFCherenkov.h:72
void clearVectors()
Definition: HFCherenkov.cc:399
bool UseNewPMT
Definition: HFCherenkov.h:63
std::vector< double > wlini
Definition: HFCherenkov.h:68
double computeHEMEff(double wavelength)
Definition: HFCherenkov.cc:357
double smearNPE(G4int Npe)
Definition: HFCherenkov.cc:385
double sinPsimax
Definition: HFCherenkov.h:61
std::vector< double > getWL()
Definition: HFCherenkov.cc:300
std::vector< double > getWLQEff()
Definition: HFCherenkov.cc:295
T sqrt(T t)
Definition: SSEVec.h:19
std::vector< double > wltrap
Definition: HFCherenkov.h:69
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
double aperture
Definition: HFCherenkov.h:60
std::vector< double > wlhem
Definition: HFCherenkov.h:71
std::vector< double > getMom()
Definition: HFCherenkov.cc:305
#define M_PI
std::vector< double > getWLIni()
Definition: HFCherenkov.cc:275
int computeNbOfPhotons(double pBeta, double step_length)
Definition: HFCherenkov.cc:310
std::vector< double > getWLTrap()
Definition: HFCherenkov.cc:280
std::vector< double > getWLAtten()
Definition: HFCherenkov.cc:285
bool isApplicable(const G4ParticleDefinition *aParticleType)
Definition: HFCherenkov.cc:409
std::vector< double > wlatten
Definition: HFCherenkov.h:70
double gain
Definition: HFCherenkov.h:61
double lambda1
Definition: HFCherenkov.h:59
int computeNPE(const G4Step *step, const G4ParticleDefinition *pDef, double pBeta, double u, double v, double w, double step_length, double zFiber, double Dose, int Npe_Dose)
Definition: HFCherenkov.cc:42
double lambda2
Definition: HFCherenkov.h:59
bool checkSurvive
Definition: HFCherenkov.h:62
tmp
align.sh
Definition: createJobs.py:716
std::vector< double > momZ
Definition: HFCherenkov.h:67
int computeNPEinPMT(const G4ParticleDefinition *pDef, double pBeta, double u, double v, double w, double step_length)
Definition: HFCherenkov.cc:199
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
double ref_index
Definition: HFCherenkov.h:58
std::vector< double > getWLHEM()
Definition: HFCherenkov.cc:290
std::vector< double > wl
Definition: HFCherenkov.h:66