Geant4-11
Public Member Functions | Static Public Member Functions | Protected Member Functions | Static Private Attributes
G4UCNMicroRoughnessHelper Class Reference

#include <G4UCNMicroRoughnessHelper.hh>

Public Member Functions

G4double Fmu (G4double, G4double, G4double, G4double, G4double, G4double, G4double) const
 
G4double FmuS (G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double) const
 
G4double IntIminus (G4double, G4double, G4double, G4int, G4int, G4double, G4double, G4double *, G4double) const
 
G4double IntIplus (G4double, G4double, G4double, G4int, G4int, G4double, G4double, G4double *, G4double) const
 
G4double ProbIminus (G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double) const
 
G4double ProbIplus (G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double) const
 
G4double S2 (G4double, G4double) const
 
G4double SS2 (G4double, G4double) const
 

Static Public Member Functions

static G4UCNMicroRoughnessHelper * GetInstance ()
 

Protected Member Functions

 G4UCNMicroRoughnessHelper ()
 
 ~G4UCNMicroRoughnessHelper ()
 

Static Private Attributes

static G4UCNMicroRoughnessHelper * fpInstance = nullptr
 

Detailed Description

Definition at line 60 of file G4UCNMicroRoughnessHelper.hh.

Constructor & Destructor Documentation

◆ G4UCNMicroRoughnessHelper()

G4UCNMicroRoughnessHelper::G4UCNMicroRoughnessHelper ( )
protected

Definition at line 53 of file G4UCNMicroRoughnessHelper.cc.

53{;}

Referenced by GetInstance().

◆ ~G4UCNMicroRoughnessHelper()

G4UCNMicroRoughnessHelper::~G4UCNMicroRoughnessHelper ( )
protected

Definition at line 57 of file G4UCNMicroRoughnessHelper.cc.

58{
59 delete fpInstance;
60 fpInstance = nullptr;
61}
G4UCNMicroRoughnessHelper::fpInstance
static G4UCNMicroRoughnessHelper * fpInstance
Definition: G4UCNMicroRoughnessHelper.hh:184

References fpInstance.

Member Function Documentation

◆ Fmu()

G4double G4UCNMicroRoughnessHelper::Fmu ( G4double  k2,
G4double  thetai,
G4double  thetao,
G4double  phio,
G4double  b2,
G4double  w2,
G4double  AngCut 
) const

Definition at line 98 of file G4UCNMicroRoughnessHelper.cc.

102{
103 G4double mu_squared;
104
105 // Checks if the distribution is peaked around the specular direction
106
107 if ((std::fabs(thetai-thetao)<AngCut) && (std::fabs(phio)<AngCut))
108 mu_squared=0.;
109 else
110 {
111 // cf. p. 177 of the Steyerl paper
112
113 G4double sinthetai=std::sin(thetai);
114 G4double sinthetao=std::sin(thetao);
115 mu_squared=k2*
116 (sinthetai*sinthetai+sinthetao*sinthetao-
117 2.*sinthetai*sinthetao*std::cos(phio));
118 }
119
120 // cf. p. 177 of the Steyerl paper
121
122 return b2*w2/twopi*std::exp(-mu_squared*w2/2);
123}
twopi
static constexpr double twopi
Definition: G4SIunits.hh:56
G4double
double G4double
Definition: G4Types.hh:83

References twopi.

Referenced by IntIplus(), and ProbIplus().

◆ FmuS()

G4double G4UCNMicroRoughnessHelper::FmuS ( G4double  k,
G4double  kS,
G4double  thetai,
G4double  thetaSo,
G4double  phiSo,
G4double  b2,
G4double  w2,
G4double  AngCut,
G4double  thetarefract 
) const

Definition at line 128 of file G4UCNMicroRoughnessHelper.cc.

133{
134 G4double mu_squared;
135
136 // Checks if the distribution is peaked around the direction of
137 // unperturbed refraction
138
139 if ((std::fabs(thetarefract-thetaSo)<AngCut) && (std::fabs(phiSo)<AngCut))
140 mu_squared=0.;
141 else
142 {
143 G4double sinthetai=std::sin(thetai);
144 G4double sinthetaSo=std::sin(thetaSo);
145
146 // cf. p. 177 of the Steyerl paper
147 mu_squared=k*k*sinthetai*sinthetai+kS*kS*sinthetaSo*sinthetaSo-
148 2.*k*kS*sinthetai*sinthetaSo*std::cos(phiSo);
149 }
150
151 // cf. p. 177 of the Steyerl paper
152
153 return b2*w2/twopi*std::exp(-mu_squared*w2/2);
154}

References twopi.

Referenced by IntIminus(), and ProbIminus().

◆ GetInstance()

G4UCNMicroRoughnessHelper * G4UCNMicroRoughnessHelper::GetInstance ( )
static

Definition at line 65 of file G4UCNMicroRoughnessHelper.cc.

66{
67 if (fpInstance == nullptr) fpInstance = new G4UCNMicroRoughnessHelper;
68 return fpInstance;
69}

References fpInstance, and G4UCNMicroRoughnessHelper().

Referenced by G4UCNMaterialPropertiesTable::ComputeMicroRoughnessTables(), G4UCNMaterialPropertiesTable::GetMRProbability(), and G4UCNMaterialPropertiesTable::GetMRTransProbability().

◆ IntIminus()

G4double G4UCNMicroRoughnessHelper::IntIminus ( G4double  E,
G4double  fermipot,
G4double  theta_i,
G4int  AngNoTheta,
G4int  AngNoPhi,
G4double  b2,
G4double  w2,
G4double *  max,
G4double  AngCut 
) const

Definition at line 283 of file G4UCNMicroRoughnessHelper.cc.

289{
290 G4double a_max_thetas_o, max_thetas_o = theta_i;
291 G4double a_max_phis_o, max_phis_o = 0.;
292 G4double thetas_o;
293 G4double phis_o;
294 G4double IntensS;
295 G4double ang_steptheta=180.*degree/(AngNoTheta-1);
296 G4double ang_stepphi=180.*degree/(AngNoPhi-1);
297 G4double costheta_i=std::cos(theta_i);
298 G4double costheta_i_squared=costheta_i*costheta_i;
299
300 *max = 0.;
301 G4double wkeit=0.;
302
303 if (E < fermipot) return wkeit;
304
305 //k_l^4/4
306 G4double kl4d4=neutron_mass_c2/hbarc_squared*neutron_mass_c2/
307 hbarc_squared*fermipot*fermipot;
308 // (k_l/k)^2
309 G4double klk2=fermipot/E;
310
311 // (k_l/k')^2
312 G4double klks2=fermipot/(E-fermipot);
313
314 // k'/k
315 G4double ksdk=std::sqrt((E-fermipot)/E);
316
317 G4double costhetas_o_squared;
318
319 // k
320 G4double k=std::sqrt(2*neutron_mass_c2*E/hbarc_squared);
321
322 // k'
323 G4double kS=ksdk*k;
324
325 for (thetas_o=0.*degree; thetas_o<=90.*degree+1e-6; thetas_o+=ang_steptheta)
326 {
327 costhetas_o_squared=std::cos(thetas_o)*std::cos(thetas_o);
328
329 for (phis_o=-180.*degree; phis_o<=180.*degree+1e-6; phis_o+=ang_stepphi)
330 {
331 //cf. Steyerl-paper p. 176, eq. 21
332 if (costhetas_o_squared>=-klks2) {
333
334 //calculates probability for a certain theta'_o, phi'_o pair
335
336 G4double thetarefract = thetas_o;
337 if (std::fabs(std::sin(theta_i)/ksdk) <= 1.)
338 thetarefract = std::asin(std::sin(theta_i)/ksdk);
339
340 IntensS = kl4d4/costheta_i*ksdk*S2(costheta_i_squared, klk2)*
341 SS2(costhetas_o_squared,klks2)*
342 FmuS(k,kS,theta_i,thetas_o,phis_o,b2,w2,AngCut,thetarefract)*
343 std::sin(thetas_o);
344 } else {
345 IntensS=0.;
346 }
347 // checks if the new probability is larger than
348 // the maximum found so far
349 if (IntensS>*max)
350 {
351 *max=IntensS;
352 }
353 wkeit+=IntensS*ang_steptheta*ang_stepphi;
354 }
355 }
356
357 // Fine-Iteration to find maximum of distribution
358
359 if (E>1e-10*eV)
360 {
361
362 // Break-condition for refining
363
364 while (ang_stepphi>=AngCut*AngCut || ang_steptheta>=AngCut*AngCut)
365 {
366 a_max_thetas_o=max_thetas_o;
367 a_max_phis_o=max_phis_o;
368 ang_stepphi /= 2.;
369 ang_steptheta /= 2.;
370 //G4cout << ang_stepphi/degree << ", " << ang_steptheta/degree
371 // << ", " << AngCut/degree << G4endl;
372 for (thetas_o=a_max_thetas_o-ang_steptheta;
373 thetas_o<=a_max_thetas_o-ang_steptheta+1e-6;
374 thetas_o+=ang_steptheta)
375 {
376 costhetas_o_squared=std::cos(thetas_o)*std::cos(thetas_o);
377 for (phis_o=a_max_phis_o-ang_stepphi;
378 phis_o<=a_max_phis_o+ang_stepphi+1e-6;
379 phis_o+=ang_stepphi)
380 {
381 G4double thetarefract = thetas_o;
382 if (std::fabs(std::sin(theta_i)/ksdk) <= 1.)
383 thetarefract = std::asin(std::sin(theta_i)/ksdk);
384
385 IntensS=kl4d4/costheta_i*ksdk*S2(costheta_i_squared, klk2)*
386 SS2(costhetas_o_squared,klks2)*
387 FmuS(k,kS,theta_i,thetas_o,phis_o,b2,w2,AngCut,thetarefract)*
388 std::sin(thetas_o);
389 if (IntensS>*max)
390 {
391 *max=IntensS;
392 max_thetas_o=thetas_o;
393 max_phis_o=phis_o;
394 }
395 }
396 }
397 }
398 }
399 return wkeit;
400}
eV
static constexpr double eV
Definition: G4SIunits.hh:201
degree
static constexpr double degree
Definition: G4SIunits.hh:124
G4UCNMicroRoughnessHelper::FmuS
G4double FmuS(G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double) const
Definition: G4UCNMicroRoughnessHelper.cc:128
G4UCNMicroRoughnessHelper::S2
G4double S2(G4double, G4double) const
Definition: G4UCNMicroRoughnessHelper.cc:74
G4UCNMicroRoughnessHelper::SS2
G4double SS2(G4double, G4double) const
Definition: G4UCNMicroRoughnessHelper.cc:88
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
source.hepunit.hbarc_squared
float hbarc_squared
Definition: hepunit.py:265
source.hepunit.neutron_mass_c2
float neutron_mass_c2
Definition: hepunit.py:275

References degree, eV, FmuS(), source.hepunit::hbarc_squared, G4INCL::Math::max(), source.hepunit::neutron_mass_c2, S2(), and SS2().

◆ IntIplus()

G4double G4UCNMicroRoughnessHelper::IntIplus ( G4double  E,
G4double  fermipot,
G4double  theta_i,
G4int  AngNoTheta,
G4int  AngNoPhi,
G4double  b2,
G4double  w2,
G4double *  max,
G4double  AngCut 
) const

Definition at line 159 of file G4UCNMicroRoughnessHelper.cc.

164{
165 *max=0.;
166
167 // max_theta_o saves the theta-position of the max probability,
168 // the previous value is saved in a_max_theta_o
169
170 G4double a_max_theta_o, max_theta_o=theta_i, a_max_phi_o, max_phi_o=0.;
171
172 // max_phi_o saves the phi-position of the max probability,
173 // the previous value is saved in a_max_phi_o
174
175 // Definition of the stepsizes in theta_o and phi_o
176
177 G4double theta_o;
178 G4double phi_o;
179 G4double Intens;
180 G4double ang_steptheta=90.*degree/(AngNoTheta-1);
181 G4double ang_stepphi=360.*degree/(AngNoPhi-1);
182 G4double costheta_i=std::cos(theta_i);
183 G4double costheta_i_squared=costheta_i*costheta_i;
184
185 // (k_l/k)^2
186 G4double kl4d4=neutron_mass_c2/hbarc_squared*neutron_mass_c2/
187 hbarc_squared*fermipot*fermipot;
188
189 // (k_l/k)^2
190 G4double klk2=fermipot/E;
191
192 G4double costheta_o_squared;
193
194 // k^2
195 G4double k2=2*neutron_mass_c2*E/hbarc_squared;
196
197 G4double wkeit=0.;
198
199 // Loop through theta_o
200
201 for (theta_o=0.*degree; theta_o<=90.*degree+1e-6; theta_o+=ang_steptheta)
202 {
203 costheta_o_squared=std::cos(theta_o)*std::cos(theta_o);
204
205 // Loop through phi_o
206
207 for (phi_o=-180.*degree; phi_o<=180.*degree+1e-6; phi_o+=ang_stepphi)
208 {
209 //calculates probability for a certain theta_o,phi_o pair
210
211 Intens=kl4d4/costheta_i*S2(costheta_i_squared,klk2)*
212 S2(costheta_o_squared,klk2)*
213 Fmu(k2,theta_i,theta_o,phi_o,b2,w2,AngCut)*std::sin(theta_o);
214
215 //G4cout << "S2: " << S2(costheta_i_squared,klk2) << G4endl;
216 //G4cout << "S2: " << S2(costheta_o_squared,klk2) << G4endl;
217 //G4cout << "Fmu: " << Fmu(k2,theta_i,theta_o,phi_o,b2,w2,AngCut)*sin(theta_o) << G4endl;
218 // checks if the new probability is larger than the
219 // maximum found so far
220
221 if (Intens>*max)
222 {
223 *max=Intens;
224 max_theta_o=theta_o;
225 max_phi_o=phi_o;
226 }
227
228 // Adds the newly found probability to the integral probability
229
230 wkeit+=Intens*ang_steptheta*ang_stepphi;
231 }
232 }
233
234 // Fine-Iteration to find maximum of distribution
235 // only if the energy is not zero
236
237 if (E>1e-10*eV)
238 {
239
240 // Break-condition for refining
241
242 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
243 while ((ang_stepphi>=AngCut*AngCut) || (ang_steptheta>=AngCut*AngCut))
244 {
245 a_max_theta_o=max_theta_o;
246 a_max_phi_o=max_phi_o;
247 ang_stepphi /= 2.;
248 ang_steptheta /= 2.;
249
250 //G4cout << ang_stepphi/degree << ", "
251 // << ang_steptheta/degree << ","
252 // << AngCut/degree << G4endl;
253
254 for (theta_o=a_max_theta_o-ang_steptheta;
255 theta_o<=a_max_theta_o-ang_steptheta+1e-6;
256 theta_o+=ang_steptheta)
257 {
258 //G4cout << "theta_o: " << theta_o/degree << G4endl;
259 costheta_o_squared=std::cos(theta_o)*std::cos(theta_o);
260 for (phi_o=a_max_phi_o-ang_stepphi;
261 phi_o<=a_max_phi_o+ang_stepphi+1e-6;
262 phi_o+=ang_stepphi)
263 {
264 //G4cout << "phi_o: " << phi_o/degree << G4endl;
265 Intens=kl4d4/costheta_i*S2(costheta_i_squared, klk2)*
266 S2(costheta_o_squared,klk2)*
267 Fmu(k2,theta_i,theta_o,phi_o,b2,w2,AngCut)*std::sin(theta_o);
268 if (Intens>*max)
269 {
270 *max=Intens;
271 max_theta_o=theta_o;
272 max_phi_o=phi_o;
273 }
274 }
275 }
276 }
277 }
278 return wkeit;
279}
G4UCNMicroRoughnessHelper::Fmu
G4double Fmu(G4double, G4double, G4double, G4double, G4double, G4double, G4double) const
Definition: G4UCNMicroRoughnessHelper.cc:98

References degree, eV, Fmu(), source.hepunit::hbarc_squared, G4INCL::Math::max(), source.hepunit::neutron_mass_c2, and S2().

◆ ProbIminus()

G4double G4UCNMicroRoughnessHelper::ProbIminus ( G4double  E,
G4double  fermipot,
G4double  theta_i,
G4double  thetas_o,
G4double  phis_o,
G4double  b,
G4double  w,
G4double  AngCut 
) const

Definition at line 431 of file G4UCNMicroRoughnessHelper.cc.

436{
437 //k_l^4/4
438 G4double kl4d4=neutron_mass_c2/hbarc_squared*neutron_mass_c2/
439 hbarc_squared*fermipot*fermipot;
440 // (k_l/k)^2
441 G4double klk2=fermipot/E;
442
443 // (k_l/k')^2
444 G4double klks2=fermipot/(E-fermipot);
445
446 if (E < fermipot) {
447 G4cout << " ProbIminus E < fermipot " << G4endl;
448 return 0.;
449 }
450
451 // k'/k
452 G4double ksdk=std::sqrt((E-fermipot)/E);
453
454 // k
455 G4double k=std::sqrt(2*neutron_mass_c2*E/hbarc_squared);
456
457 // k'/k
458 G4double kS=ksdk*k;
459
460 G4double costheta_i=std::cos(theta_i);
461 G4double costhetas_o=std::cos(thetas_o);
462
463 // eq. 20 on page 176 in the steyerl paper
464
465 G4double thetarefract = thetas_o;
466 if(std::fabs(std::sin(theta_i)/ksdk) <= 1.)thetarefract = std::asin(std::sin(theta_i)/ksdk);
467
468 return kl4d4/costheta_i*ksdk*S2(costheta_i*costheta_i, klk2)*
469 SS2(costhetas_o*costhetas_o,klks2)*
470 FmuS(k,kS,theta_i,thetas_o,phis_o,b*b,w*w,AngCut,thetarefract)*
471 std::sin(thetas_o);
472}
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout

References FmuS(), G4cout, G4endl, source.hepunit::hbarc_squared, source.hepunit::neutron_mass_c2, S2(), and SS2().

◆ ProbIplus()

G4double G4UCNMicroRoughnessHelper::ProbIplus ( G4double  E,
G4double  fermipot,
G4double  theta_i,
G4double  theta_o,
G4double  phi_o,
G4double  b,
G4double  w,
G4double  AngCut 
) const

Definition at line 404 of file G4UCNMicroRoughnessHelper.cc.

410{
411 //k_l^4/4
412 G4double kl4d4=neutron_mass_c2/hbarc_squared*neutron_mass_c2/
413 hbarc_squared*fermipot*fermipot;
414
415 // (k_l/k)^2
416 G4double klk2=fermipot/E;
417
418 G4double costheta_i=std::cos(theta_i);
419 G4double costheta_o=std::cos(theta_o);
420
421 // eq. 20 on page 176 in the steyerl paper
422
423 return kl4d4/costheta_i*S2(costheta_i*costheta_i, klk2)*
424 S2(costheta_o*costheta_o,klk2)*
425 Fmu(2*neutron_mass_c2*E/hbarc_squared,theta_i,theta_o,phi_o,b*b,w*w,AngCut)*
426 std::sin(theta_o);
427}

References Fmu(), source.hepunit::hbarc_squared, source.hepunit::neutron_mass_c2, and S2().

◆ S2()

G4double G4UCNMicroRoughnessHelper::S2 ( G4double  costheta2,
G4double  klk2 
) const

Definition at line 74 of file G4UCNMicroRoughnessHelper.cc.

75{
76 // case 1: radicand is positive,
77 // case 2: radicand is negative, cf. p. 174 of the Steyerl paper
78
79 if (costheta2>=klk2)
80 return 4*costheta2/(2*costheta2-klk2+2*std::sqrt(costheta2*(costheta2-klk2)));
81 else
82 return std::norm(2*std::sqrt(costheta2)/(std::sqrt(costheta2) + std::sqrt(std::complex<G4double>(costheta2 - klk2))));
83}

Referenced by IntIminus(), IntIplus(), ProbIminus(), and ProbIplus().

◆ SS2()

G4double G4UCNMicroRoughnessHelper::SS2 ( G4double  costhetas2,
G4double  klks2 
) const

Definition at line 88 of file G4UCNMicroRoughnessHelper.cc.

89{
90 // cf. p. 175 of the Steyerl paper
91
92 return 4*costhetas2/
93 (2*costhetas2+klks2+2*std::sqrt(costhetas2*(costhetas2+klks2)));
94}

Referenced by IntIminus(), and ProbIminus().

Field Documentation

◆ fpInstance

G4UCNMicroRoughnessHelper * G4UCNMicroRoughnessHelper::fpInstance = nullptr
staticprivate

Definition at line 184 of file G4UCNMicroRoughnessHelper.hh.

Referenced by GetInstance(), and ~G4UCNMicroRoughnessHelper().

The documentation for this class was generated from the following files:
Generated by doxygen 1.9.3