79 return 4*costheta2/(2*costheta2-klk2+2*std::sqrt(costheta2*(costheta2-klk2)));
81 return 4*costheta2/klk2;
91 (2*costhetas2+klks2+2*std::sqrt(costhetas2*(costhetas2+klks2)));
105 if ((std::fabs(thetai-thetao)<AngCut) && (std::fabs(phio)<AngCut))
111 G4double sinthetai=std::sin(thetai);
112 G4double sinthetao=std::sin(thetao);
114 (sinthetai*sinthetai+sinthetao*sinthetao-
115 2.*sinthetai*sinthetao*std::cos(phio));
120 return b2*w2/
twopi*std::exp(-mu_squared*w2/2);
136 if ((std::fabs(thetarefract-thetaSo)<AngCut) && (std::fabs(phiSo)<AngCut))
140 G4double sinthetai=std::sin(thetai);
141 G4double sinthetaSo=std::sin(thetaSo);
144 mu_squared=k*k*sinthetai*sinthetai+kS*kS*sinthetaSo*sinthetaSo-
145 2.*k*kS*sinthetai*sinthetaSo*std::cos(phiSo);
150 return b2*w2/
twopi*std::exp(-mu_squared*w2/2);
166 G4double a_max_theta_o, max_theta_o=theta_i, a_max_phi_o, max_phi_o=0.;
178 G4double costheta_i=std::cos(theta_i);
179 G4double costheta_i_squared=costheta_i*costheta_i;
182 G4double kl4d4=neutron_mass_c2/hbarc_squared*neutron_mass_c2/
183 hbarc_squared*fermipot*fermipot;
191 G4double k2=2*neutron_mass_c2*E/hbarc_squared;
197 for (theta_o=0.*
degree; theta_o<=90.*
degree+1e-6; theta_o+=ang_steptheta)
199 costheta_o_squared=std::cos(theta_o)*std::cos(theta_o);
203 for (phi_o=-180.*
degree; phi_o<=180.*
degree+1e-6; phi_o+=ang_stepphi)
207 Intens=kl4d4/costheta_i*
S2(costheta_i_squared,klk2)*
208 S2(costheta_o_squared,klk2)*
209 Fmu(k2,theta_i,theta_o,phi_o,b2,w2,AngCut)*std::sin(theta_o);
226 wkeit+=Intens*ang_steptheta*ang_stepphi;
239 while ((ang_stepphi>=AngCut*AngCut) || (ang_steptheta>=AngCut*AngCut))
241 a_max_theta_o=max_theta_o;
242 a_max_phi_o=max_phi_o;
250 for (theta_o=a_max_theta_o-ang_steptheta;
251 theta_o<=a_max_theta_o-ang_steptheta+1e-6;
252 theta_o+=ang_steptheta)
255 costheta_o_squared=std::cos(theta_o)*std::cos(theta_o);
256 for (phi_o=a_max_phi_o-ang_stepphi;
257 phi_o<=a_max_phi_o+ang_stepphi+1e-6;
261 Intens=kl4d4/costheta_i*
S2(costheta_i_squared, klk2)*
262 S2(costheta_o_squared,klk2)*
263 Fmu(k2,theta_i,theta_o,phi_o,b2,w2,AngCut)*std::sin(theta_o);
286 G4double a_max_thetas_o, max_thetas_o = theta_i;
287 G4double a_max_phis_o, max_phis_o = 0.;
293 G4double costheta_i=std::cos(theta_i);
294 G4double costheta_i_squared=costheta_i*costheta_i;
299 if (E*costheta_i_squared < fermipot)
return wkeit;
302 G4double kl4d4=neutron_mass_c2/hbarc_squared*neutron_mass_c2/
303 hbarc_squared*fermipot*fermipot;
308 G4double klks2=fermipot/(E-fermipot);
311 G4double ksdk=std::sqrt((E-fermipot)/E);
316 G4double k=std::sqrt(2*neutron_mass_c2*E/hbarc_squared);
321 for (thetas_o=0.*
degree; thetas_o<=90.*
degree+1e-6; thetas_o+=ang_steptheta)
323 costhetas_o_squared=std::cos(thetas_o)*std::cos(thetas_o);
325 for (phis_o=-180.*
degree; phis_o<=180.*
degree+1e-6; phis_o+=ang_stepphi)
328 if (costhetas_o_squared>=-klks2) {
333 if (std::fabs(std::sin(theta_i)/ksdk) <= 1.)
334 thetarefract = std::asin(std::sin(theta_i)/ksdk);
336 IntensS = kl4d4/costheta_i*ksdk*
S2(costheta_i_squared, klk2)*
337 SS2(costhetas_o_squared,klks2)*
338 FmuS(k,kS,theta_i,thetas_o,phis_o,b2,w2,AngCut,thetarefract)*
349 wkeit+=IntensS*ang_steptheta*ang_stepphi;
360 while (ang_stepphi>=AngCut*AngCut || ang_steptheta>=AngCut*AngCut)
362 a_max_thetas_o=max_thetas_o;
363 a_max_phis_o=max_phis_o;
368 for (thetas_o=a_max_thetas_o-ang_steptheta;
369 thetas_o<=a_max_thetas_o-ang_steptheta+1e-6;
370 thetas_o+=ang_steptheta)
372 costhetas_o_squared=std::cos(thetas_o)*std::cos(thetas_o);
373 for (phis_o=a_max_phis_o-ang_stepphi;
374 phis_o<=a_max_phis_o+ang_stepphi+1e-6;
378 if (std::fabs(std::sin(theta_i)/ksdk) <= 1.)
379 thetarefract = std::asin(std::sin(theta_i)/ksdk);
381 IntensS=kl4d4/costheta_i*ksdk*
S2(costheta_i_squared, klk2)*
382 SS2(costhetas_o_squared,klks2)*
383 FmuS(k,kS,theta_i,thetas_o,phis_o,b2,w2,AngCut,thetarefract)*
388 max_thetas_o=thetas_o;
408 G4double kl4d4=neutron_mass_c2/hbarc_squared*neutron_mass_c2/
409 hbarc_squared*fermipot*fermipot;
414 G4double costheta_i=std::cos(theta_i);
415 G4double costheta_o=std::cos(theta_o);
419 return kl4d4/costheta_i*
S2(costheta_i*costheta_i, klk2)*
420 S2(costheta_o*costheta_o,klk2)*
421 Fmu(2*neutron_mass_c2*E/hbarc_squared,theta_i,theta_o,phi_o,b*b,w*w,AngCut)*
434 G4double kl4d4=neutron_mass_c2/hbarc_squared*neutron_mass_c2/
435 hbarc_squared*fermipot*fermipot;
440 G4double klks2=fermipot/(E-fermipot);
448 G4double ksdk=std::sqrt((E-fermipot)/E);
451 G4double k=std::sqrt(2*neutron_mass_c2*E/hbarc_squared);
456 G4double costheta_i=std::cos(theta_i);
457 G4double costhetas_o=std::cos(thetas_o);
462 if(std::fabs(std::sin(theta_i)/ksdk) <= 1.)thetarefract = std::asin(std::sin(theta_i)/ksdk);
464 return kl4d4/costheta_i*ksdk*
S2(costheta_i*costheta_i, klk2)*
465 SS2(costhetas_o*costhetas_o,klks2)*
466 FmuS(k,kS,theta_i,thetas_o,phis_o,b*b,w*w,AngCut,thetarefract)*
G4double Fmu(G4double, G4double, G4double, G4double, G4double, G4double, G4double)
G4UCNMicroRoughnessHelper()
G4double IntIplus(G4double, G4double, G4double, G4int, G4int, G4double, G4double, G4double *, G4double)
const G4double w[NPOINTSGL]
G4double ProbIplus(G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double)
G4GLOB_DLL std::ostream G4cout
G4double S2(G4double, G4double)
static const double twopi
G4double FmuS(G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double)
static G4UCNMicroRoughnessHelper * GetInstance()
~G4UCNMicroRoughnessHelper()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static const double degree
static G4UCNMicroRoughnessHelper * fpInstance
G4double SS2(G4double, G4double)
G4double ProbIminus(G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double)
G4double IntIminus(G4double, G4double, G4double, G4int, G4int, G4double, G4double, G4double *, G4double)