| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25
26 // This file contains the source code of vario
27 // calculation of microroughness.
28 //
29 // see A. Steyerl, Z. Physik 254 (1972) 169.
30 //
31 // A description of the functions can be found
32
33 #include "G4UCNMicroRoughnessHelper.hh"
34
35 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "globals.hh"
38
39 G4UCNMicroRoughnessHelper* G4UCNMicroRoughness
40
41 //....oooOO0OOooo........oooOO0OOooo........oo
42
43 G4UCNMicroRoughnessHelper::~G4UCNMicroRoughnes
44
45 //....oooOO0OOooo........oooOO0OOooo........oo
46
47 G4UCNMicroRoughnessHelper* G4UCNMicroRoughness
48 {
49 if (fpInstance == nullptr) {
50 fpInstance = new G4UCNMicroRoughnessHelper
51 }
52 return fpInstance;
53 }
54
55 //....oooOO0OOooo........oooOO0OOooo........oo
56
57 G4double G4UCNMicroRoughnessHelper::S2(G4doubl
58 {
59 // case 1: radicand is positive,
60 // case 2: radicand is negative, cf. p. 174
61
62 if (costheta2 >= klk2) {
63 return 4 * costheta2 / (2 * costheta2 - kl
64 }
65
66 return std::norm(2 * std::sqrt(costheta2) /
67 (std::sqrt(costheta2) + std
68 }
69
70 //....oooOO0OOooo........oooOO0OOooo........oo
71
72 G4double G4UCNMicroRoughnessHelper::SS2(G4doub
73 {
74 // cf. p. 175 of the Steyerl paper
75
76 return 4 * costhetas2 /
77 (2 * costhetas2 + klks2 + 2 * std::sq
78 }
79
80 //....oooOO0OOooo........oooOO0OOooo........oo
81
82 G4double G4UCNMicroRoughnessHelper::Fmu(G4doub
83 G4double phio, G4double b2, G4double w2, G4d
84 {
85 G4double mu_squared;
86
87 // Checks if the distribution is peaked arou
88
89 if ((std::fabs(thetai - thetao) < 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 +
98 2. * sinthetai * sinth
99 }
100
101 // cf. p. 177 of the Steyerl paper
102
103 return b2 * w2 / twopi * std::exp(-mu_square
104 }
105
106 //....oooOO0OOooo........oooOO0OOooo........oo
107
108 G4double G4UCNMicroRoughnessHelper::FmuS(G4dou
109 G4double phiSo, G4double b2, G4double w2, G4
110 {
111 G4double mu_squared;
112
113 // Checks if the distribution is peaked arou
114 // unperturbed refraction
115
116 if ((std::fabs(thetarefract - thetaSo) < Ang
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
125 2. * k * kS * sinthetai * sin
126 }
127
128 // cf. p. 177 of the Steyerl paper
129
130 return b2 * w2 / twopi * std::exp(-mu_square
131 }
132
133 //....oooOO0OOooo........oooOO0OOooo........oo
134
135 G4double G4UCNMicroRoughnessHelper::IntIplus(G
136 G4int AngNoTheta, G4int AngNoPhi, G4double b
137 {
138 *max = 0.;
139
140 // max_theta_o saves the theta-position of t
141 // the previous value is saved in a_max_thet
142
143 G4double a_max_theta_o, max_theta_o = theta_
144
145 // max_phi_o saves the phi-position of the m
146 // the previous value is saved in a_max_phi_
147
148 // Definition of the stepsizes in theta_o an
149
150 G4double theta_o;
151 G4double phi_o;
152 G4double Intens;
153 G4double ang_steptheta = 90. * degree / (Ang
154 G4double ang_stepphi = 360. * degree / (AngN
155 G4double costheta_i = std::cos(theta_i);
156 G4double costheta_i_squared = costheta_i * c
157
158 // (k_l/k)^2
159 G4double kl4d4 =
160 neutron_mass_c2 / hbarc_squared * neutron_
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 / hbar
169
170 G4double wkeit = 0.;
171
172 // Loop through theta_o
173
174 for (theta_o = 0. * degree; theta_o <= 90. *
175 costheta_o_squared = std::cos(theta_o) * s
176
177 // Loop through phi_o
178
179 for (phi_o = -180. * degree; phi_o <= 180.
180 // calculates probability for a certain
181
182 Intens = kl4d4 / costheta_i * S2(costhet
183 Fmu(k2, theta_i, theta_o, phi_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 t
192
193 wkeit += Intens * ang_steptheta * ang_st
194 }
195 }
196
197 // Fine-Iteration to find maximum of distrib
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 Iv
204 while ((ang_stepphi >= 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_stept
211 theta_o += ang_steptheta)
212 {
213 // G4cout << "theta_o: " << theta_o/de
214 costheta_o_squared = std::cos(theta_o)
215 for (phi_o = a_max_phi_o - ang_stepphi
216 phi_o += ang_stepphi)
217 {
218 // G4cout << "phi_o: " << phi_o/degr
219 Intens = kl4d4 / costheta_i * S2(cos
220 S2(costheta_o_squared, klk2
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........oo
235
236 G4double G4UCNMicroRoughnessHelper::IntIminus(
237 G4int AngNoTheta, G4int AngNoPhi, G4double b
238 {
239 G4double a_max_thetas_o, max_thetas_o = thet
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 / (An
245 G4double ang_stepphi = 180. * degree / (AngN
246 G4double costheta_i = std::cos(theta_i);
247 G4double costheta_i_squared = costheta_i * c
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_
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 *
272
273 // k'
274 G4double kS = ksdk * k;
275
276 for (thetas_o = 0. * degree; thetas_o <= 90.
277 costhetas_o_squared = std::cos(thetas_o) *
278
279 for (phis_o = -180. * degree; phis_o <= 18
280 // cf. Steyerl-paper p. 176, eq. 21
281 if (costhetas_o_squared >= -klks2) {
282 // calculates probability for a certai
283
284 G4double thetarefract = thetas_o;
285 if (std::fabs(std::sin(theta_i) / ksdk
286 thetarefract = std::asin(std::sin(th
287 }
288
289 IntensS = kl4d4 / costheta_i * ksdk *
290 SS2(costhetas_o_squared, klk
291 FmuS(k, kS, theta_i, thetas_
292 std::sin(thetas_o);
293 }
294 else {
295 IntensS = 0.;
296 }
297 // checks if the new probability is larg
298 // the maximum found so far
299 if (IntensS > *max) {
300 *max = IntensS;
301 }
302 wkeit += IntensS * ang_steptheta * ang_s
303 }
304 }
305
306 // Fine-Iteration to find maximum of distrib
307
308 if (E > 1e-10 * eV) {
309 // Break-condition for refining
310
311 while (ang_stepphi >= AngCut * AngCut || a
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_ste
318 thetas_o <= a_max_thetas_o - ang_st
319 {
320 costhetas_o_squared = std::cos(thetas_
321 for (phis_o = a_max_phis_o - ang_stepp
322 phis_o += ang_stepphi)
323 {
324 G4double thetarefract = thetas_o;
325 if (std::fabs(std::sin(theta_i) / ks
326 thetarefract = std::asin(std::sin(
327 }
328
329 IntensS = kl4d4 / costheta_i * ksdk
330 SS2(costhetas_o_squared, k
331 FmuS(k, kS, theta_i, theta
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........oo
346
347 G4double G4UCNMicroRoughnessHelper::ProbIplus(
348 G4double theta_o, G4double phi_o, G4double b
349 {
350 // k_l^4/4
351 G4double kl4d4 =
352 neutron_mass_c2 / hbarc_squared * neutron_
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 *
363 S2(costheta_o * costheta_o, klk2) *
364 Fmu(
365 2 * neutron_mass_c2 * E / hbarc_squ
366 std::sin(theta_o);
367 }
368
369 //....oooOO0OOooo........oooOO0OOooo........oo
370
371 G4double G4UCNMicroRoughnessHelper::ProbIminus
372 G4double thetas_o, G4double phis_o, G4double
373 {
374 // k_l^4/4
375 G4double kl4d4 =
376 neutron_mass_c2 / hbarc_squared * neutron_
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 " << G
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 *
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)
405 }
406
407 return kl4d4 / costheta_i * ksdk * S2(costhe
408 SS2(costhetas_o * costhetas_o, klks2)
409 FmuS(k, kS, theta_i, thetas_o, phis_o
410 std::sin(thetas_o);
411 }
412