| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25
26 // This file contains the source code of various functions all related to the
27 // calculation of microroughness.
28 //
29 // see A. Steyerl, Z. Physik 254 (1972) 169.
30 //
31 // A description of the functions can be found in the corresponding header file
32
33 #include "G4UCNMicroRoughnessHelper.hh"
34
35 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "globals.hh"
38
39 G4UCNMicroRoughnessHelper* G4UCNMicroRoughnessHelper::fpInstance = nullptr;
40
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42
43 G4UCNMicroRoughnessHelper::~G4UCNMicroRoughnessHelper() { fpInstance = nullptr; }
44
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46
47 G4UCNMicroRoughnessHelper* G4UCNMicroRoughnessHelper::GetInstance()
48 {
49 if (fpInstance == nullptr) {
50 fpInstance = new G4UCNMicroRoughnessHelper;
51 }
52 return fpInstance;
53 }
54
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56
57 G4double G4UCNMicroRoughnessHelper::S2(G4double costheta2, G4double klk2) const
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 }
69
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71
72 G4double G4UCNMicroRoughnessHelper::SS2(G4double costhetas2, G4double klks2) const
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 }
79
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81
82 G4double G4UCNMicroRoughnessHelper::Fmu(G4double k2, G4double thetai, G4double thetao,
83 G4double phio, G4double b2, G4double w2, G4double AngCut) const
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 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107
108 G4double G4UCNMicroRoughnessHelper::FmuS(G4double k, G4double kS, G4double thetai, G4double thetaSo,
109 G4double phiSo, G4double b2, G4double w2, G4double AngCut, G4double thetarefract) const
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 }
132
133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
134
135 G4double G4UCNMicroRoughnessHelper::IntIplus(G4double E, G4double fermipot, G4double theta_i,
136 G4int AngNoTheta, G4int AngNoPhi, G4double b2, G4double w2, G4double* max, G4double AngCut) const
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 }
233
234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
235
236 G4double G4UCNMicroRoughnessHelper::IntIminus(G4double E, G4double fermipot, G4double theta_i,
237 G4int AngNoTheta, G4int AngNoPhi, G4double b2, G4double w2, G4double* max, G4double AngCut) const
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 }
344
345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
346
347 G4double G4UCNMicroRoughnessHelper::ProbIplus(G4double E, G4double fermipot, G4double theta_i,
348 G4double theta_o, G4double phi_o, G4double b, G4double w, G4double AngCut) const
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 }
368
369 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
370
371 G4double G4UCNMicroRoughnessHelper::ProbIminus(G4double E, G4double fermipot, G4double theta_i,
372 G4double thetas_o, G4double phis_o, G4double b, G4double w, G4double AngCut) const
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 }
412