Geant4 11.2.2
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 ()=default
 
 ~G4UCNMicroRoughnessHelper ()
 

Detailed Description

Definition at line 54 of file G4UCNMicroRoughnessHelper.hh.

Constructor & Destructor Documentation

◆ G4UCNMicroRoughnessHelper()

G4UCNMicroRoughnessHelper::G4UCNMicroRoughnessHelper ( )
protecteddefault

Referenced by GetInstance().

◆ ~G4UCNMicroRoughnessHelper()

G4UCNMicroRoughnessHelper::~G4UCNMicroRoughnessHelper ( )
protected

Definition at line 43 of file G4UCNMicroRoughnessHelper.cc.

43{ fpInstance = nullptr; }

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

84{
85 G4double mu_squared;
86
87 // Checks if the distribution is peaked around the specular direction
88
89 if ((std::fabs(thetai - thetao) < AngCut) && (std::fabs(phio) < AngCut)) {
90 mu_squared = 0.;
91 }
92 else {
93 // cf. p. 177 of the Steyerl paper
94
95 G4double sinthetai = std::sin(thetai);
96 G4double sinthetao = std::sin(thetao);
97 mu_squared = k2 * (sinthetai * sinthetai + sinthetao * sinthetao -
98 2. * sinthetai * sinthetao * std::cos(phio));
99 }
100
101 // cf. p. 177 of the Steyerl paper
102
103 return b2 * w2 / twopi * std::exp(-mu_squared * w2 / 2);
104}
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 108 of file G4UCNMicroRoughnessHelper.cc.

110{
111 G4double mu_squared;
112
113 // Checks if the distribution is peaked around the direction of
114 // unperturbed refraction
115
116 if ((std::fabs(thetarefract - thetaSo) < AngCut) && (std::fabs(phiSo) < AngCut)) {
117 mu_squared = 0.;
118 }
119 else {
120 G4double sinthetai = std::sin(thetai);
121 G4double sinthetaSo = std::sin(thetaSo);
122
123 // cf. p. 177 of the Steyerl paper
124 mu_squared = k * k * sinthetai * sinthetai + kS * kS * sinthetaSo * sinthetaSo -
125 2. * k * kS * sinthetai * sinthetaSo * std::cos(phiSo);
126 }
127
128 // cf. p. 177 of the Steyerl paper
129
130 return b2 * w2 / twopi * std::exp(-mu_squared * w2 / 2);
131}

Referenced by IntIminus(), and ProbIminus().

◆ GetInstance()

G4UCNMicroRoughnessHelper * G4UCNMicroRoughnessHelper::GetInstance ( )
static

Definition at line 47 of file G4UCNMicroRoughnessHelper.cc.

48{
49 if (fpInstance == nullptr) {
50 fpInstance = new G4UCNMicroRoughnessHelper;
51 }
52 return fpInstance;
53}
G4UCNMicroRoughnessHelper::G4UCNMicroRoughnessHelper
G4UCNMicroRoughnessHelper()=default

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

238{
239 G4double a_max_thetas_o, max_thetas_o = theta_i;
240 G4double a_max_phis_o, max_phis_o = 0.;
241 G4double thetas_o;
242 G4double phis_o;
243 G4double IntensS;
244 G4double ang_steptheta = 180. * degree / (AngNoTheta - 1);
245 G4double ang_stepphi = 180. * degree / (AngNoPhi - 1);
246 G4double costheta_i = std::cos(theta_i);
247 G4double costheta_i_squared = costheta_i * costheta_i;
248
249 *max = 0.;
250 G4double wkeit = 0.;
251
252 if (E < fermipot) {
253 return wkeit;
254 }
255
256 // k_l^4/4
257 G4double kl4d4 =
258 neutron_mass_c2 / hbarc_squared * neutron_mass_c2 / hbarc_squared * fermipot * fermipot;
259 // (k_l/k)^2
260 G4double klk2 = fermipot / E;
261
262 // (k_l/k')^2
263 G4double klks2 = fermipot / (E - fermipot);
264
265 // k'/k
266 G4double ksdk = std::sqrt((E - fermipot) / E);
267
268 G4double costhetas_o_squared;
269
270 // k
271 G4double k = std::sqrt(2 * neutron_mass_c2 * E / hbarc_squared);
272
273 // k'
274 G4double kS = ksdk * k;
275
276 for (thetas_o = 0. * degree; thetas_o <= 90. * degree + 1e-6; thetas_o += ang_steptheta) {
277 costhetas_o_squared = std::cos(thetas_o) * std::cos(thetas_o);
278
279 for (phis_o = -180. * degree; phis_o <= 180. * degree + 1e-6; phis_o += ang_stepphi) {
280 // cf. Steyerl-paper p. 176, eq. 21
281 if (costhetas_o_squared >= -klks2) {
282 // calculates probability for a certain theta'_o, phi'_o pair
283
284 G4double thetarefract = thetas_o;
285 if (std::fabs(std::sin(theta_i) / ksdk) <= 1.) {
286 thetarefract = std::asin(std::sin(theta_i) / ksdk);
287 }
288
289 IntensS = kl4d4 / costheta_i * ksdk * S2(costheta_i_squared, klk2) *
290 SS2(costhetas_o_squared, klks2) *
291 FmuS(k, kS, theta_i, thetas_o, phis_o, b2, w2, AngCut, thetarefract) *
292 std::sin(thetas_o);
293 }
294 else {
295 IntensS = 0.;
296 }
297 // checks if the new probability is larger than
298 // the maximum found so far
299 if (IntensS > *max) {
300 *max = IntensS;
301 }
302 wkeit += IntensS * ang_steptheta * ang_stepphi;
303 }
304 }
305
306 // Fine-Iteration to find maximum of distribution
307
308 if (E > 1e-10 * eV) {
309 // Break-condition for refining
310
311 while (ang_stepphi >= AngCut * AngCut || ang_steptheta >= AngCut * AngCut) {
312 a_max_thetas_o = max_thetas_o;
313 a_max_phis_o = max_phis_o;
314 ang_stepphi /= 2.;
315 ang_steptheta /= 2.;
316
317 for (thetas_o = a_max_thetas_o - ang_steptheta;
318 thetas_o <= a_max_thetas_o - ang_steptheta + 1e-6; thetas_o += ang_steptheta)
319 {
320 costhetas_o_squared = std::cos(thetas_o) * std::cos(thetas_o);
321 for (phis_o = a_max_phis_o - ang_stepphi; phis_o <= a_max_phis_o + ang_stepphi + 1e-6;
322 phis_o += ang_stepphi)
323 {
324 G4double thetarefract = thetas_o;
325 if (std::fabs(std::sin(theta_i) / ksdk) <= 1.) {
326 thetarefract = std::asin(std::sin(theta_i) / ksdk);
327 }
328
329 IntensS = kl4d4 / costheta_i * ksdk * S2(costheta_i_squared, klk2) *
330 SS2(costhetas_o_squared, klks2) *
331 FmuS(k, kS, theta_i, thetas_o, phis_o, b2, w2, AngCut, thetarefract) *
332 std::sin(thetas_o);
333 if (IntensS > *max) {
334 *max = IntensS;
335 max_thetas_o = thetas_o;
336 max_phis_o = phis_o;
337 }
338 }
339 }
340 }
341 }
342 return wkeit;
343}
G4UCNMicroRoughnessHelper::FmuS
G4double FmuS(G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double) const
Definition G4UCNMicroRoughnessHelper.cc:108
G4UCNMicroRoughnessHelper::S2
G4double S2(G4double, G4double) const
Definition G4UCNMicroRoughnessHelper.cc:57
G4UCNMicroRoughnessHelper::SS2
G4double SS2(G4double, G4double) const
Definition G4UCNMicroRoughnessHelper.cc:72
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition G4INCLGlobals.hh:112

Referenced by G4UCNMaterialPropertiesTable::ComputeMicroRoughnessTables().

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

137{
138 *max = 0.;
139
140 // max_theta_o saves the theta-position of the max probability,
141 // the previous value is saved in a_max_theta_o
142
143 G4double a_max_theta_o, max_theta_o = theta_i, a_max_phi_o, max_phi_o = 0.;
144
145 // max_phi_o saves the phi-position of the max probability,
146 // the previous value is saved in a_max_phi_o
147
148 // Definition of the stepsizes in theta_o and phi_o
149
150 G4double theta_o;
151 G4double phi_o;
152 G4double Intens;
153 G4double ang_steptheta = 90. * degree / (AngNoTheta - 1);
154 G4double ang_stepphi = 360. * degree / (AngNoPhi - 1);
155 G4double costheta_i = std::cos(theta_i);
156 G4double costheta_i_squared = costheta_i * costheta_i;
157
158 // (k_l/k)^2
159 G4double kl4d4 =
160 neutron_mass_c2 / hbarc_squared * neutron_mass_c2 / hbarc_squared * fermipot * fermipot;
161
162 // (k_l/k)^2
163 G4double klk2 = fermipot / E;
164
165 G4double costheta_o_squared;
166
167 // k^2
168 G4double k2 = 2 * neutron_mass_c2 * E / hbarc_squared;
169
170 G4double wkeit = 0.;
171
172 // Loop through theta_o
173
174 for (theta_o = 0. * degree; theta_o <= 90. * degree + 1e-6; theta_o += ang_steptheta) {
175 costheta_o_squared = std::cos(theta_o) * std::cos(theta_o);
176
177 // Loop through phi_o
178
179 for (phi_o = -180. * degree; phi_o <= 180. * degree + 1e-6; phi_o += ang_stepphi) {
180 // calculates probability for a certain theta_o,phi_o pair
181
182 Intens = kl4d4 / costheta_i * S2(costheta_i_squared, klk2) * S2(costheta_o_squared, klk2) *
183 Fmu(k2, theta_i, theta_o, phi_o, b2, w2, AngCut) * std::sin(theta_o);
184
185 if (Intens > *max) {
186 *max = Intens;
187 max_theta_o = theta_o;
188 max_phi_o = phi_o;
189 }
190
191 // Adds the newly found probability to the integral probability
192
193 wkeit += Intens * ang_steptheta * ang_stepphi;
194 }
195 }
196
197 // Fine-Iteration to find maximum of distribution
198 // only if the energy is not zero
199
200 if (E > 1e-10 * eV) {
201 // Break-condition for refining
202
203 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
204 while ((ang_stepphi >= AngCut * AngCut) || (ang_steptheta >= AngCut * AngCut)) {
205 a_max_theta_o = max_theta_o;
206 a_max_phi_o = max_phi_o;
207 ang_stepphi /= 2.;
208 ang_steptheta /= 2.;
209
210 for (theta_o = a_max_theta_o - ang_steptheta; theta_o <= a_max_theta_o - ang_steptheta + 1e-6;
211 theta_o += ang_steptheta)
212 {
213 // G4cout << "theta_o: " << theta_o/degree << G4endl;
214 costheta_o_squared = std::cos(theta_o) * std::cos(theta_o);
215 for (phi_o = a_max_phi_o - ang_stepphi; phi_o <= a_max_phi_o + ang_stepphi + 1e-6;
216 phi_o += ang_stepphi)
217 {
218 // G4cout << "phi_o: " << phi_o/degree << G4endl;
219 Intens = kl4d4 / costheta_i * S2(costheta_i_squared, klk2) *
220 S2(costheta_o_squared, klk2) * Fmu(k2, theta_i, theta_o, phi_o, b2, w2, AngCut) *
221 std::sin(theta_o);
222 if (Intens > *max) {
223 *max = Intens;
224 max_theta_o = theta_o;
225 max_phi_o = phi_o;
226 }
227 }
228 }
229 }
230 }
231 return wkeit;
232}
G4UCNMicroRoughnessHelper::Fmu
G4double Fmu(G4double, G4double, G4double, G4double, G4double, G4double, G4double) const
Definition G4UCNMicroRoughnessHelper.cc:82

Referenced by G4UCNMaterialPropertiesTable::ComputeMicroRoughnessTables().

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

373{
374 // k_l^4/4
375 G4double kl4d4 =
376 neutron_mass_c2 / hbarc_squared * neutron_mass_c2 / hbarc_squared * fermipot * fermipot;
377 // (k_l/k)^2
378 G4double klk2 = fermipot / E;
379
380 // (k_l/k')^2
381 G4double klks2 = fermipot / (E - fermipot);
382
383 if (E < fermipot) {
384 G4cout << " ProbIminus E < fermipot " << G4endl;
385 return 0.;
386 }
387
388 // k'/k
389 G4double ksdk = std::sqrt((E - fermipot) / E);
390
391 // k
392 G4double k = std::sqrt(2 * neutron_mass_c2 * E / hbarc_squared);
393
394 // k'/k
395 G4double kS = ksdk * k;
396
397 G4double costheta_i = std::cos(theta_i);
398 G4double costhetas_o = std::cos(thetas_o);
399
400 // eq. 20 on page 176 in the steyerl paper
401
402 G4double thetarefract = thetas_o;
403 if (std::fabs(std::sin(theta_i) / ksdk) <= 1.) {
404 thetarefract = std::asin(std::sin(theta_i) / ksdk);
405 }
406
407 return kl4d4 / costheta_i * ksdk * S2(costheta_i * costheta_i, klk2) *
408 SS2(costhetas_o * costhetas_o, klks2) *
409 FmuS(k, kS, theta_i, thetas_o, phis_o, b * b, w * w, AngCut, thetarefract) *
410 std::sin(thetas_o);
411}
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout

Referenced by G4UCNMaterialPropertiesTable::GetMRTransProbability().

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

349{
350 // k_l^4/4
351 G4double kl4d4 =
352 neutron_mass_c2 / hbarc_squared * neutron_mass_c2 / hbarc_squared * fermipot * fermipot;
353
354 // (k_l/k)^2
355 G4double klk2 = fermipot / E;
356
357 G4double costheta_i = std::cos(theta_i);
358 G4double costheta_o = std::cos(theta_o);
359
360 // eq. 20 on page 176 in the steyerl paper
361
362 return kl4d4 / costheta_i * S2(costheta_i * costheta_i, klk2) *
363 S2(costheta_o * costheta_o, klk2) *
364 Fmu(
365 2 * neutron_mass_c2 * E / hbarc_squared, theta_i, theta_o, phi_o, b * b, w * w, AngCut) *
366 std::sin(theta_o);
367}

Referenced by G4UCNMaterialPropertiesTable::GetMRProbability().

◆ S2()

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

Definition at line 57 of file G4UCNMicroRoughnessHelper.cc.

58{
59 // case 1: radicand is positive,
60 // case 2: radicand is negative, cf. p. 174 of the Steyerl paper
61
62 if (costheta2 >= klk2) {
63 return 4 * costheta2 / (2 * costheta2 - klk2 + 2 * std::sqrt(costheta2 * (costheta2 - klk2)));
64 }
65
66 return std::norm(2 * std::sqrt(costheta2) /
67 (std::sqrt(costheta2) + std::sqrt(std::complex<G4double>(costheta2 - klk2))));
68}

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

◆ SS2()

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

Definition at line 72 of file G4UCNMicroRoughnessHelper.cc.

73{
74 // cf. p. 175 of the Steyerl paper
75
76 return 4 * costhetas2 /
77 (2 * costhetas2 + klks2 + 2 * std::sqrt(costhetas2 * (costhetas2 + klks2)));
78}

Referenced by IntIminus(), and ProbIminus().

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