8 #include "G4Poisson.hh"
9 #include "G4ParticleDefinition.hh"
10 #include "G4NavigationHistory.hh"
29 edm::LogInfo(
"HFShower") <<
"HFCherenkov:: initialised with ref_index "
31 << lambda1/cm <<
"|" <<
lambda2/cm
32 <<
" aperture(total/trapped) " <<
aperture <<
"|"
34 <<
" Check photon survival in HF " << checkSurvive
35 <<
" Gain " << gain <<
" useNewPMT " << UseNewPMT
44 double u,
double v,
double w,
45 double step_length,
double zFiber,
46 double dose,
int npe_Dose) {
48 if (pBeta < (1/
ref_index) || step_length < 0.0001) {
return 0;}
50 double uv =
sqrt(u*u + v*v);
53 if (nbOfPhotons < 0) {
55 }
else if (nbOfPhotons > 0) {
57 for (
int i = 0;
i < nbOfPhotons;
i++) {
58 double rand = G4UniformRand();
59 double theta_C = acos(1./(pBeta*
ref_index));
61 double sinTheta =
sin(theta_C);
62 double cosTheta =
cos(theta_C);
63 double cosPhi =
cos(phi_C);
68 w_ph = w * cosTheta - sinTheta * cosPhi * uv;
75 int n_photons = npe_Dose;
80 double pBeta,
double u,
double v,
double w,
81 double step_length,
double zFiber,
82 double dose,
int npe_Dose) {
86 if (pBeta < (1/
ref_index) || step_length < 0.0001) {
88 LogDebug(
"HFShower") <<
"HFCherenkov::computeNPE: pBeta " << pBeta
89 <<
" 1/mu " << (1/
ref_index) <<
" step_length "
95 double uv =
sqrt(u*u + v*v);
97 *aStep->GetTrack()->GetWeight();
99 LogDebug(
"HFShower") <<
"HFCherenkov::computeNPE: pBeta " << pBeta
100 <<
" u/v/w " << u <<
"/" << v <<
"/" << w
101 <<
" step_length " << step_length <<
" zFib " << zFiber
102 <<
" nbOfPhotons " << nbOfPhotons;
104 if (nbOfPhotons < 0) {
106 }
else if (nbOfPhotons > 0) {
107 G4StepPoint * preStepPoint = aStep->GetPreStepPoint();
108 G4TouchableHandle theTouchable = preStepPoint->GetTouchableHandle();
109 G4ThreeVector localprepos = theTouchable->GetHistory()->
110 GetTopTransform().TransformPoint(aStep->GetPreStepPoint()->GetPosition());
111 G4ThreeVector localpostpos = theTouchable->GetHistory()->
112 GetTopTransform().TransformPoint(aStep->GetPostStepPoint()->GetPosition());
114 double length=
sqrt((localpostpos.x()-localprepos.x())*(localpostpos.x()-localprepos.x())
115 +(localpostpos.y()-localprepos.y())*(localpostpos.y()-localprepos.y()));
118 double u_ph=0,v_ph=0, w_ph=0;
119 for (
int i = 0;
i < nbOfPhotons;
i++) {
120 double rand = G4UniformRand();
121 double theta_C = acos(1./(pBeta*
ref_index));
123 double sinTheta =
sin(theta_C);
124 double cosTheta =
cos(theta_C);
125 double cosPhi =
cos(phi_C);
126 double sinPhi =
sin(phi_C);
129 u_ph = sinTheta * cosPhi ;
130 v_ph = sinTheta * sinPhi;
133 u_ph = uv * cosTheta + sinTheta * cosPhi *
w;
134 v_ph = sinTheta * sinPhi;
135 w_ph = w * cosTheta - sinTheta * cosPhi * uv;
137 double r_lambda = G4UniformRand();
140 double lambda = (lambda0/cm) *
pow(
double(10),7);
141 wlini.push_back(lambda);
143 LogDebug(
"HFShower") <<
"HFCherenkov::computeNPE: " <<
i <<
" lambda "
144 << lambda <<
" w_ph " << w_ph <<
" aperture "
148 double xemit=length*(G4UniformRand()-0.5);
149 double gam=atan2(yemit,xemit);
150 double eps=atan2(v_ph,u_ph);
151 double sinBeta=
sin(gam-eps);
152 double rho=
sqrt(xemit*xemit+yemit*yemit);
153 double sinEta=rho/
fibreR*sinBeta;
154 double cosEta=
sqrt(1.-sinEta*sinEta);
155 double sinPsi=
sqrt(1.-w_ph*w_ph);
156 double cosKsi=cosEta*sinPsi;
158 if (cosKsi < aperturetrapped && w_ph>0.) {
159 LogDebug(
"HFShower") <<
"HFCherenkov::Trapped photon : " << u_ph <<
" "
160 << v_ph <<
" " << w_ph <<
" " << xemit <<
" "
161 << gam <<
" " << eps <<
" " << sinBeta <<
" "
162 << rho <<
" " << sinEta <<
" " << cosEta <<
" "
163 <<
" " << sinPsi <<
" " << cosKsi;
165 LogDebug(
"HFShower") <<
"HFCherenkov::Rejected photon : " << u_ph <<
" "
166 << v_ph <<
" " << w_ph <<
" " << xemit <<
" "
167 << gam <<
" " << eps <<
" " << sinBeta <<
" "
168 << rho <<
" " << sinEta <<
" " << cosEta <<
" "
169 <<
" " << sinPsi <<
" " << cosKsi;
176 double prob_HF = 1.0;
177 double a0_inv = 0.1234;
178 double prob_MX =
exp( - 0.5 * a0_inv );
180 double a_inv = a0_inv + 0.14 *
pow(dose,0.30);
181 double z_meters = zFiber;
182 prob_HF =
exp(-z_meters * a_inv );
184 rand = G4UniformRand();
186 LogDebug(
"HFShower") <<
"HFCherenkov::computeNPE: probHF " << prob_HF
187 <<
" prob_MX " << prob_MX <<
" Random " << rand
188 <<
" Survive? " << (rand < (prob_HF * prob_MX));
190 if (rand < (prob_HF * prob_MX)) {
192 rand = G4UniformRand();
195 LogDebug(
"HFShower") <<
"HFCherenkov::computeNPE: w_ph " << w_ph
196 <<
" effHEM " << effHEM <<
" Random " << rand
197 <<
" Survive? " << (w_ph>0.997||(rand<effHEM));
199 if (w_ph>0.997 || (rand<effHEM)) {
200 wlhem.push_back(lambda);
202 rand = G4UniformRand();
204 LogDebug(
"HFShower") <<
"HFCherenkov::computeNPE: qEffic "
205 << qEffic <<
" Random " << rand
206 <<
" Survive? " <<(rand < qEffic);
210 momZ.push_back(w_ph);
211 wl.push_back(lambda);
221 LogDebug(
"HFShower") <<
"HFCherenkov::computeNPE: npe " << npe;
227 double u,
double v,
double w,
232 if (pBeta < (1/
ref_index) || step_length < 0.0001) {
234 LogDebug(
"HFShower") <<
"HFCherenkov::computeNPEinPMT: pBeta " << pBeta
235 <<
" 1/mu " << (1/
ref_index) <<
" step_length "
241 double uv =
sqrt(u*u + v*v);
244 LogDebug(
"HFShower") <<
"HFCherenkov::computeNPEinPMT: pBeta " << pBeta
245 <<
" u/v/w " << u <<
"/" << v <<
"/" << w
246 <<
" step_length " << step_length
247 <<
" nbOfPhotons " << nbOfPhotons;
249 if (nbOfPhotons < 0) {
251 }
else if (nbOfPhotons > 0) {
253 for (
int i = 0;
i < nbOfPhotons;
i++) {
254 double rand = G4UniformRand();
255 double theta_C = acos(1./(pBeta*
ref_index));
257 double sinTheta =
sin(theta_C);
258 double cosTheta =
cos(theta_C);
259 double cosPhi =
cos(phi_C);
264 w_ph = w * cosTheta - sinTheta * cosPhi * uv;
266 double r_lambda = G4UniformRand();
269 double lambda = (lambda0/cm) *
pow(
double(10),7);
270 wlini.push_back(lambda);
272 LogDebug(
"HFShower") <<
"HFCherenkov::computeNPEinPMT: " <<
i <<
" lambda "
273 << lambda <<
" w_ph " << w_ph <<
" aperture "
276 if (w_ph > aperture) {
278 rand = G4UniformRand();
280 LogDebug(
"HFShower") <<
"HFCherenkov::computeNPEinPMT: Random " << rand
281 <<
" Survive? " << (rand < 1.);
285 rand = G4UniformRand();
287 rand = G4UniformRand();
289 LogDebug(
"HFShower") <<
"HFCherenkov::computeNPEinPMT: qEffic "
290 << qEffic <<
" Random " << rand
291 <<
" Survive? " <<(rand < qEffic);
295 momZ.push_back(w_ph);
296 wl.push_back(lambda);
304 LogDebug(
"HFShower") <<
"HFCherenkov::computeNPEinPMT: npe " << npe_;
310 std::vector<double>
v =
wlini;
315 std::vector<double>
v =
wltrap;
325 std::vector<double>
v =
wlhem;
330 std::vector<double>
v =
wlqeff;
335 std::vector<double>
v =
wl;
340 std::vector<double>
v =
momZ;
347 double alpha = 0.0073;
348 double step_length = stepL;
349 double theta_C = acos(1./(pBeta*
ref_index));
351 double cherenPhPerLength = 2 *
M_PI * alpha * lambdaDiff*cm;
352 double d_NOfPhotons = cherenPhPerLength *
sin(theta_C)*
sin(theta_C) * (step_length/cm);
353 int nbOfPhotons = int(d_NOfPhotons);
355 LogDebug(
"HFShower") <<
"HFCherenkov::computeNbOfPhotons: StepLength "
356 << step_length <<
" theta_C " << theta_C
357 <<
" lambdaDiff " << lambdaDiff
358 <<
" cherenPhPerLength " << cherenPhPerLength
359 <<
" Photons " << d_NOfPhotons <<
" " << nbOfPhotons;
368 if (wavelength<=350) {
369 qeff=2.45867*(TMath::Landau(wavelength,353.820,59.1324));
370 }
else if (wavelength>350 && wavelength<500) {
371 qeff= 0.441989*
exp(-
pow((wavelength-358.371),2)/(2*
pow((138.277),2)));
372 }
else if (wavelength>=500 && wavelength<550) {
373 qeff= 0.271862*
exp(-
pow((wavelength-491.505),2)/(2*
pow((47.0418),2)));
374 }
else if (wavelength>=550) {
375 qeff= 0.137297*
exp(-
pow((wavelength-520.260),2)/(2*
pow((75.5023),2)));
378 LogDebug(
"HFShower") <<
"HFCherenkov:: for new PMT : wavelength === "
379 << wavelength <<
"\tqeff ===\t" << qeff;
382 double y = (wavelength - 275.) /180.;
383 double func = 1. / (1. + 250.*
pow((y/5.),4));
384 double qE_R7525 = 0.77 * y *
exp(-y) * func ;
387 LogDebug(
"HFShower") <<
"HFCherenkov:: for old PMT : wavelength === "
388 << wavelength <<
"; qeff = " << qeff;
393 LogDebug(
"HFShower") <<
"HFCherenkov::computeQEff: wavelength " << wavelength
394 <<
" y/func " <<
y <<
"/" << func <<
" qeff " << qeff;
402 if (wavelength < 400.) {
404 }
else if (wavelength >= 400. && wavelength < 410.) {
406 hEMEff = (-1322.453 / wavelength) + 4.2056;
407 }
else if (wavelength >= 410.) {
409 if (wavelength > 415. && wavelength < 445.) {
413 }
else if (wavelength > 550. && wavelength < 600.) {
417 }
else if (wavelength > 565. && wavelength <= 635.) {
419 hEMEff = (701.7268 / wavelength) - 0.186;
423 LogDebug(
"HFShower") <<
"HFCherenkov::computeHEMEff: wavelength "
424 << wavelength <<
" hEMEff " << hEMEff;
433 for (
int i = 0;
i < npe; ++
i) {
434 double val = G4Poisson(
gain);
435 pe += (val/
gain) + 0.001*G4UniformRand();
439 LogDebug(
"HFShower") <<
"HFCherenkov::smearNPE: npe " << npe <<
" pe " << pe;
456 bool tmp = (aParticleType->GetPDGCharge() != 0);
458 LogDebug(
"HFShower") <<
"HFCherenkov::isApplicable: aParticleType "
459 << aParticleType->GetParticleName() <<
" PDGCharge "
460 << aParticleType->GetPDGCharge() <<
" Result " <<
tmp;
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
double computeQEff(double wavelength)
Sin< T >::type sin(const T &t)
HFCherenkov(edm::ParameterSet const &p)
std::vector< double > wlqeff
int computeNPEinPMT(G4ParticleDefinition *pDef, double pBeta, double u, double v, double w, double step_length)
std::vector< double > wlini
double computeHEMEff(double wavelength)
double smearNPE(G4int Npe)
std::vector< double > getWL()
std::vector< double > getWLQEff()
std::vector< double > wltrap
Cos< T >::type cos(const T &t)
std::vector< double > wlhem
std::vector< double > getMom()
std::vector< double > getWLIni()
std::vector< std::vector< double > > tmp
int computeNbOfPhotons(double pBeta, double step_length)
std::vector< double > getWLTrap()
std::vector< double > getWLAtten()
bool isApplicable(const G4ParticleDefinition *aParticleType)
std::vector< double > wlatten
std::vector< double > momZ
int computeNPhTrapped(double pBeta, double u, double v, double w, double step_length, double zFiber, double Dose, int Npe_Dose)
Power< A, B >::type pow(const A &a, const B &b)
int computeNPE(G4Step *step, G4ParticleDefinition *pDef, double pBeta, double u, double v, double w, double step_length, double zFiber, double Dose, int Npe_Dose)
std::vector< double > getWLHEM()