8 #include "G4Poisson.hh"
9 #include "G4ParticleDefinition.hh"
10 #include "G4NavigationHistory.hh"
13 #include "Randomize.hh"
14 #include "G4SystemOfUnits.hh"
34 <<
" Check photon survival in HF " <<
checkSurvive <<
" Gain " <<
gain <<
" useNewPMT "
43 const G4ParticleDefinition* pDef,
56 if (pBeta < (1 /
ref_index) || step_length < 0.0001) {
59 <<
" step_length " << step_length;
65 int nbOfPhotons =
computeNbOfPhotons(pBeta, step_length) * aStep->GetTrack()->GetWeight();
67 edm::LogVerbatim(
"HFShower") <<
"HFCherenkov::computeNPE: pBeta " << pBeta <<
" u/v/w " << u <<
"/" <<
v <<
"/" <<
w
68 <<
" step_length " << step_length <<
" zFib " << zFiber <<
" nbOfPhotons "
71 if (nbOfPhotons < 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());
81 double length =
std::sqrt((localpostpos.x() - localprepos.x()) * (localpostpos.x() - localprepos.x()) +
82 (localpostpos.y() - localprepos.y()) * (localpostpos.y() - localprepos.y()));
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);
100 u_ph = uv * cosTheta + sinTheta *
cosPhi *
w;
102 w_ph =
w * cosTheta - sinTheta *
cosPhi * uv;
104 double r_lambda = G4UniformRand();
106 double lambda = (lambda0 / CLHEP::cm) * 1.
e+7;
107 wlini.push_back(lambda);
109 edm::LogVerbatim(
"HFShower") <<
"HFCherenkov::computeNPE: " <<
i <<
" lambda " << lambda <<
" w_ph " << w_ph;
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);
118 double cosEta =
std::sqrt(1. - sinEta * sinEta);
119 double sinPsi =
std::sqrt(1. - w_ph * w_ph);
120 double cosKsi = cosEta * sinPsi;
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;
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;
139 double prob_HF = 1.0;
141 double a0_inv = 0.1234;
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);
146 rand = G4UniformRand();
148 edm::LogVerbatim(
"HFShower") <<
"HFCherenkov::computeNPE: probHF " << prob_HF <<
" Random " << rand
149 <<
" Survive? " << (rand < prob_HF);
151 if (rand < prob_HF) {
153 rand = G4UniformRand();
156 rand = G4UniformRand();
157 double phi = 0.5 *
M_PI * rand;
158 double rad_bundle = 19.;
160 rand = G4UniformRand();
161 double rad = rad_bundle *
std::sqrt(rand);
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);
170 edm::LogVerbatim(
"HFShower") <<
"HFCherenkov::computeNPE: w_ph " << w_ph <<
" effHEM " << effHEM <<
" eff "
171 << eff <<
" Random " << rand <<
" Survive? " << (rand < eff);
174 wlhem.push_back(lambda);
176 rand = G4UniformRand();
178 edm::LogVerbatim(
"HFShower") <<
"HFCherenkov::computeNPE: qEffic " << qEffic <<
" Random " << rand
179 <<
" Survive? " << (rand < qEffic);
183 momZ.push_back(w_ph);
184 wl.push_back(lambda);
200 const G4ParticleDefinition* pDef,
double pBeta,
double u,
double v,
double w,
double step_length) {
206 if (pBeta < (1 /
ref_index) || step_length < 0.0001) {
209 <<
" step_length " << step_length;
217 edm::LogVerbatim(
"HFShower") <<
"HFCherenkov::computeNPEinPMT: pBeta " << pBeta <<
" u/v/w " << u <<
"/" <<
v <<
"/"
218 <<
w <<
" step_length " << step_length <<
" nbOfPhotons " << nbOfPhotons;
220 if (nbOfPhotons < 0) {
222 }
else if (nbOfPhotons > 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);
235 w_ph =
w * cosTheta - sinTheta *
cosPhi * uv;
237 double r_lambda = G4UniformRand();
239 double lambda = (lambda0 / CLHEP::cm) * 1.
e+7;
240 wlini.push_back(lambda);
242 edm::LogVerbatim(
"HFShower") <<
"HFCherenkov::computeNPEinPMT: " <<
i <<
" lambda " << lambda <<
" w_ph " << w_ph
247 rand = G4UniformRand();
249 edm::LogVerbatim(
"HFShower") <<
"HFCherenkov::computeNPEinPMT: Random " << rand <<
" Survive? " << (rand < 1.);
254 rand = G4UniformRand();
256 edm::LogVerbatim(
"HFShower") <<
"HFCherenkov::computeNPEinPMT: qEffic " << qEffic <<
" Random " << rand
257 <<
" Survive? " << (rand < qEffic);
261 momZ.push_back(w_ph);
262 wl.push_back(lambda);
270 edm::LogVerbatim(
"HFShower") <<
"HFCherenkov::computeNPEinPMT: npe " << npe_;
276 std::vector<double>
v =
wlini;
281 std::vector<double>
v =
wltrap;
291 std::vector<double>
v =
wlhem;
296 std::vector<double>
v =
wlqeff;
301 std::vector<double>
v =
wl;
306 std::vector<double>
v =
momZ;
312 double alpha = 0.0073;
313 double step_length = stepL;
314 double theta_C = acos(1. / (pBeta *
ref_index));
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);
320 edm::LogVerbatim(
"HFShower") <<
"HFCherenkov::computeNbOfPhotons: StepLength " << step_length <<
" theta_C "
321 << theta_C <<
" lambdaDiff " << lambdaDiff <<
" cherenPhPerLength " << cherenPhPerLength
322 <<
" Photons " << d_NOfPhotons <<
" " << nbOfPhotons;
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)));
340 edm::LogVerbatim(
"HFShower") <<
"HFCherenkov:: for new PMT : wavelength === " << wavelength <<
"\tqeff ===\t"
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;
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;
359 if (wavelength < 400.) {
361 }
else if (wavelength >= 400. && wavelength < 410.) {
363 hEMEff = (-1322.453 / wavelength) + 4.2056;
364 }
else if (wavelength >= 410.) {
366 if (wavelength > 415. && wavelength < 445.) {
370 }
else if (wavelength > 550. && wavelength < 600.) {
374 }
else if (wavelength > 565. && wavelength <= 635.) {
376 hEMEff = (701.7268 / wavelength) - 0.186;
380 edm::LogVerbatim(
"HFShower") <<
"HFCherenkov::computeHEMEff: wavelength " << wavelength <<
" hEMEff " << hEMEff;
388 for (
int i = 0;
i < npe; ++
i) {
390 pe += (
val /
gain) + 0.001 * G4UniformRand();
394 edm::LogVerbatim(
"HFShower") <<
"HFCherenkov::smearNPE: npe " << npe <<
" pe " << pe;
410 bool tmp = (aParticleType->GetPDGCharge() != 0);
412 edm::LogVerbatim(
"HFShower") <<
"HFCherenkov::isApplicable: aParticleType " << aParticleType->GetParticleName()
413 <<
" PDGCharge " << aParticleType->GetPDGCharge() <<
" Result " <<
tmp;