Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | Static Public Member Functions | Protected Member Functions | List of all members
G4UCNMicroRoughnessHelper Class Reference

#include <G4UCNMicroRoughnessHelper.hh>

Public Member Functions

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

Static Public Member Functions

static G4UCNMicroRoughnessHelper * GetInstance ()
 

Protected Member Functions

 G4UCNMicroRoughnessHelper ()
 
 ~G4UCNMicroRoughnessHelper ()
 

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 fpInstance = nullptr;
60}

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 108 of file G4UCNMicroRoughnessHelper.cc.

112{
113 G4double mu_squared;
114
115 // Checks if the distribution is peaked around the specular direction
116
117 if((std::fabs(thetai - thetao) < AngCut) && (std::fabs(phio) < AngCut))
118 {
119 mu_squared=0.;
120 }
121 else
122 {
123 // cf. p. 177 of the Steyerl paper
124
125 G4double sinthetai = std::sin(thetai);
126 G4double sinthetao = std::sin(thetao);
127 mu_squared = k2 * (sinthetai * sinthetai + sinthetao * sinthetao -
128 2. * sinthetai * sinthetao * std::cos(phio));
129 }
130
131 // cf. p. 177 of the Steyerl paper
132
133 return b2*w2/twopi*std::exp(-mu_squared*w2/2);
134}
G4double
double G4double
Definition: G4Types.hh:83

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 139 of file G4UCNMicroRoughnessHelper.cc.

144{
145 G4double mu_squared;
146
147 // Checks if the distribution is peaked around the direction of
148 // unperturbed refraction
149
150 if((std::fabs(thetarefract - thetaSo) < AngCut) &&
151 (std::fabs(phiSo) < AngCut))
152 {
153 mu_squared=0.;
154 }
155 else
156 {
157 G4double sinthetai = std::sin(thetai);
158 G4double sinthetaSo = std::sin(thetaSo);
159
160 // cf. p. 177 of the Steyerl paper
161 mu_squared = k * k * sinthetai * sinthetai +
162 kS * kS * sinthetaSo * sinthetaSo -
163 2. * k * kS * sinthetai * sinthetaSo * std::cos(phiSo);
164 }
165
166 // cf. p. 177 of the Steyerl paper
167
168 return b2*w2/twopi*std::exp(-mu_squared*w2/2);
169}

Referenced by IntIminus(), and ProbIminus().

◆ GetInstance()

G4UCNMicroRoughnessHelper * G4UCNMicroRoughnessHelper::GetInstance ( )
static

Definition at line 64 of file G4UCNMicroRoughnessHelper.cc.

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

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 298 of file G4UCNMicroRoughnessHelper.cc.

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

◆ 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 174 of file G4UCNMicroRoughnessHelper.cc.

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

◆ 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 453 of file G4UCNMicroRoughnessHelper.cc.

458{
459 //k_l^4/4
460 G4double kl4d4=neutron_mass_c2/hbarc_squared*neutron_mass_c2/
461 hbarc_squared*fermipot*fermipot;
462 // (k_l/k)^2
463 G4double klk2=fermipot/E;
464
465 // (k_l/k')^2
466 G4double klks2=fermipot/(E-fermipot);
467
468 if (E < fermipot) {
469 G4cout << " ProbIminus E < fermipot " << G4endl;
470 return 0.;
471 }
472
473 // k'/k
474 G4double ksdk=std::sqrt((E-fermipot)/E);
475
476 // k
477 G4double k=std::sqrt(2*neutron_mass_c2*E/hbarc_squared);
478
479 // k'/k
480 G4double kS=ksdk*k;
481
482 G4double costheta_i=std::cos(theta_i);
483 G4double costhetas_o=std::cos(thetas_o);
484
485 // eq. 20 on page 176 in the steyerl paper
486
487 G4double thetarefract = thetas_o;
488 if(std::fabs(std::sin(theta_i) / ksdk) <= 1.)
489 {
490 thetarefract = std::asin(std::sin(theta_i) / ksdk);
491 }
492
493 return kl4d4/costheta_i*ksdk*S2(costheta_i*costheta_i, klk2)*
494 SS2(costhetas_o*costhetas_o,klks2)*
495 FmuS(k,kS,theta_i,thetas_o,phis_o,b*b,w*w,AngCut,thetarefract)*
496 std::sin(thetas_o);
497}
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout

◆ 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 426 of file G4UCNMicroRoughnessHelper.cc.

432{
433 //k_l^4/4
434 G4double kl4d4=neutron_mass_c2/hbarc_squared*neutron_mass_c2/
435 hbarc_squared*fermipot*fermipot;
436
437 // (k_l/k)^2
438 G4double klk2=fermipot/E;
439
440 G4double costheta_i=std::cos(theta_i);
441 G4double costheta_o=std::cos(theta_o);
442
443 // eq. 20 on page 176 in the steyerl paper
444
445 return kl4d4/costheta_i*S2(costheta_i*costheta_i, klk2)*
446 S2(costheta_o*costheta_o,klk2)*
447 Fmu(2*neutron_mass_c2*E/hbarc_squared,theta_i,theta_o,phi_o,b*b,w*w,AngCut)*
448 std::sin(theta_o);
449}

◆ S2()

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

Definition at line 76 of file G4UCNMicroRoughnessHelper.cc.

77{
78 // case 1: radicand is positive,
79 // case 2: radicand is negative, cf. p. 174 of the Steyerl paper
80
81 if(costheta2 >= klk2)
82 {
83 return 4 * costheta2 /
84 (2 * costheta2 - klk2 +
85 2 * std::sqrt(costheta2 * (costheta2 - klk2)));
86 }
87 else
88 {
89 return std::norm(2 * std::sqrt(costheta2) /
90 (std::sqrt(costheta2) +
91 std::sqrt(std::complex<G4double>(costheta2 - klk2))));
92 }
93}

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

◆ SS2()

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

Definition at line 98 of file G4UCNMicroRoughnessHelper.cc.

99{
100 // cf. p. 175 of the Steyerl paper
101
102 return 4*costhetas2/
103 (2*costhetas2+klks2+2*std::sqrt(costhetas2*(costhetas2+klks2)));
104}

Referenced by IntIminus(), and ProbIminus().

The documentation for this class was generated from the following files: