Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4UCNMicroRoughnessHelper.cc
Go to the documentation of this file.
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//
27//
28//
29//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
30//
31// This file contains the source code of various functions all related to the
32// calculation of microroughness.
33//
34// see A. Steyerl, Z. Physik 254 (1972) 169.
35//
36// A description of the functions can be found in the corresponding header file
37//
38//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
39
40#include "G4UCNMicroRoughnessHelper.hh"
41
42#include "globals.hh"
43
44#include "G4SystemOfUnits.hh"
45#include "G4PhysicalConstants.hh"
46
47G4UCNMicroRoughnessHelper* G4UCNMicroRoughnessHelper::fpInstance = nullptr;
48
49//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50
51// Constructor
52
53G4UCNMicroRoughnessHelper::G4UCNMicroRoughnessHelper() {;}
54
55//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56
57G4UCNMicroRoughnessHelper::~G4UCNMicroRoughnessHelper()
58{
59 fpInstance = nullptr;
60}
61
62//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
63
64G4UCNMicroRoughnessHelper* G4UCNMicroRoughnessHelper::GetInstance()
65{
66 if(fpInstance == nullptr)
67 {
68 fpInstance = new G4UCNMicroRoughnessHelper;
69 }
70 return fpInstance;
71}
72
73//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74
75G4double
76G4UCNMicroRoughnessHelper::S2(G4double costheta2, G4double klk2) const
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}
94
95//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
96
97G4double
98G4UCNMicroRoughnessHelper::SS2(G4double costhetas2, G4double klks2) const
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}
105
106//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107
108G4double G4UCNMicroRoughnessHelper::Fmu(G4double k2, G4double thetai,
109 G4double thetao, G4double phio,
110 G4double b2, G4double w2,
111 G4double AngCut) const
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}
135
136//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
137
138G4double
139G4UCNMicroRoughnessHelper::FmuS(G4double k, G4double kS,
140 G4double thetai, G4double thetaSo,
141 G4double phiSo,
142 G4double b2, G4double w2,
143 G4double AngCut, G4double thetarefract) const
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}
170
171//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
172
173G4double
174G4UCNMicroRoughnessHelper::IntIplus(G4double E, G4double fermipot,
175 G4double theta_i, G4int AngNoTheta,
176 G4int AngNoPhi, G4double b2,
177 G4double w2, G4double* max,
178 G4double AngCut) const
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}
295
296//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
297
298G4double G4UCNMicroRoughnessHelper::IntIminus(G4double E, G4double fermipot,
299 G4double theta_i,
300 G4int AngNoTheta,
301 G4int AngNoPhi, G4double b2,
302 G4double w2, G4double* max,
303 G4double AngCut) const
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}
423
424//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
425
426G4double G4UCNMicroRoughnessHelper::ProbIplus(G4double E, G4double fermipot,
427 G4double theta_i,
428 G4double theta_o,
429 G4double phi_o,
430 G4double b, G4double w,
431 G4double AngCut) const
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}
450
451//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
452
453G4double G4UCNMicroRoughnessHelper::ProbIminus(G4double E, G4double fermipot,
454 G4double theta_i,
455 G4double thetas_o,
456 G4double phis_o, G4double b,
457 G4double w, G4double AngCut) const
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}
498
499//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4UCNMicroRoughnessHelper::ProbIplus
G4double ProbIplus(G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double) const
Definition: G4UCNMicroRoughnessHelper.cc:426
G4UCNMicroRoughnessHelper::FmuS
G4double FmuS(G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double) const
Definition: G4UCNMicroRoughnessHelper.cc:139
G4UCNMicroRoughnessHelper::IntIminus
G4double IntIminus(G4double, G4double, G4double, G4int, G4int, G4double, G4double, G4double *, G4double) const
Definition: G4UCNMicroRoughnessHelper.cc:298
G4UCNMicroRoughnessHelper::S2
G4double S2(G4double, G4double) const
Definition: G4UCNMicroRoughnessHelper.cc:76
G4UCNMicroRoughnessHelper::SS2
G4double SS2(G4double, G4double) const
Definition: G4UCNMicroRoughnessHelper.cc:98
G4UCNMicroRoughnessHelper::GetInstance
static G4UCNMicroRoughnessHelper * GetInstance()
Definition: G4UCNMicroRoughnessHelper.cc:64
G4UCNMicroRoughnessHelper::Fmu
G4double Fmu(G4double, G4double, G4double, G4double, G4double, G4double, G4double) const
Definition: G4UCNMicroRoughnessHelper.cc:108
G4UCNMicroRoughnessHelper::ProbIminus
G4double ProbIminus(G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double) const
Definition: G4UCNMicroRoughnessHelper.cc:453
G4UCNMicroRoughnessHelper::IntIplus
G4double IntIplus(G4double, G4double, G4double, G4int, G4int, G4double, G4double, G4double *, G4double) const
Definition: G4UCNMicroRoughnessHelper.cc:174