Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
41 
42 #include "globals.hh"
43 
44 #include "G4SystemOfUnits.hh"
45 #include "G4PhysicalConstants.hh"
46 
47 G4UCNMicroRoughnessHelper* G4UCNMicroRoughnessHelper::fpInstance = nullptr;
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 
51 // Constructor
52 
54 
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 
58 {
59  delete fpInstance;
60  fpInstance = nullptr;
61 }
62 
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64 
66 {
67  if (fpInstance == nullptr) fpInstance = new G4UCNMicroRoughnessHelper;
68  return fpInstance;
69 }
70 
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72 
73 G4double
75 {
76  // case 1: radicand is positive,
77  // case 2: radicand is negative, cf. p. 174 of the Steyerl paper
78 
79  if (costheta2>=klk2)
80  return 4*costheta2/(2*costheta2-klk2+2*std::sqrt(costheta2*(costheta2-klk2)));
81  else
82  return std::norm(2*std::sqrt(costheta2)/(std::sqrt(costheta2) + std::sqrt(std::complex<G4double>(costheta2 - klk2))));
83 }
84 
85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
86 
87 G4double
89 {
90  // cf. p. 175 of the Steyerl paper
91 
92  return 4*costhetas2/
93  (2*costhetas2+klks2+2*std::sqrt(costhetas2*(costhetas2+klks2)));
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
97 
99  G4double thetao, G4double phio,
100  G4double b2, G4double w2,
101  G4double AngCut) const
102 {
103  G4double mu_squared;
104 
105  // Checks if the distribution is peaked around the specular direction
106 
107  if ((std::fabs(thetai-thetao)<AngCut) && (std::fabs(phio)<AngCut))
108  mu_squared=0.;
109  else
110  {
111  // cf. p. 177 of the Steyerl paper
112 
113  G4double sinthetai=std::sin(thetai);
114  G4double sinthetao=std::sin(thetao);
115  mu_squared=k2*
116  (sinthetai*sinthetai+sinthetao*sinthetao-
117  2.*sinthetai*sinthetao*std::cos(phio));
118  }
119 
120  // cf. p. 177 of the Steyerl paper
121 
122  return b2*w2/twopi*std::exp(-mu_squared*w2/2);
123 }
124 
125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
126 
127 G4double
129  G4double thetai, G4double thetaSo,
130  G4double phiSo,
131  G4double b2, G4double w2,
132  G4double AngCut, G4double thetarefract) const
133 {
134  G4double mu_squared;
135 
136  // Checks if the distribution is peaked around the direction of
137  // unperturbed refraction
138 
139  if ((std::fabs(thetarefract-thetaSo)<AngCut) && (std::fabs(phiSo)<AngCut))
140  mu_squared=0.;
141  else
142  {
143  G4double sinthetai=std::sin(thetai);
144  G4double sinthetaSo=std::sin(thetaSo);
145 
146  // cf. p. 177 of the Steyerl paper
147  mu_squared=k*k*sinthetai*sinthetai+kS*kS*sinthetaSo*sinthetaSo-
148  2.*k*kS*sinthetai*sinthetaSo*std::cos(phiSo);
149  }
150 
151  // cf. p. 177 of the Steyerl paper
152 
153  return b2*w2/twopi*std::exp(-mu_squared*w2/2);
154 }
155 
156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
157 
158 G4double
160  G4double theta_i, G4int AngNoTheta,
161  G4int AngNoPhi, G4double b2,
162  G4double w2, G4double* max,
163  G4double AngCut) const
164 {
165  *max=0.;
166 
167  // max_theta_o saves the theta-position of the max probability,
168  // the previous value is saved in a_max_theta_o
169 
170  G4double a_max_theta_o, max_theta_o=theta_i, a_max_phi_o, max_phi_o=0.;
171 
172  // max_phi_o saves the phi-position of the max probability,
173  // the previous value is saved in a_max_phi_o
174 
175  // Definition of the stepsizes in theta_o and phi_o
176 
177  G4double theta_o;
178  G4double phi_o;
179  G4double Intens;
180  G4double ang_steptheta=90.*degree/(AngNoTheta-1);
181  G4double ang_stepphi=360.*degree/(AngNoPhi-1);
182  G4double costheta_i=std::cos(theta_i);
183  G4double costheta_i_squared=costheta_i*costheta_i;
184 
185  // (k_l/k)^2
187  hbarc_squared*fermipot*fermipot;
188 
189  // (k_l/k)^2
190  G4double klk2=fermipot/E;
191 
192  G4double costheta_o_squared;
193 
194  // k^2
196 
197  G4double wkeit=0.;
198 
199  // Loop through theta_o
200 
201  for (theta_o=0.*degree; theta_o<=90.*degree+1e-6; theta_o+=ang_steptheta)
202  {
203  costheta_o_squared=std::cos(theta_o)*std::cos(theta_o);
204 
205  // Loop through phi_o
206 
207  for (phi_o=-180.*degree; phi_o<=180.*degree+1e-6; phi_o+=ang_stepphi)
208  {
209  //calculates probability for a certain theta_o,phi_o pair
210 
211  Intens=kl4d4/costheta_i*S2(costheta_i_squared,klk2)*
212  S2(costheta_o_squared,klk2)*
213  Fmu(k2,theta_i,theta_o,phi_o,b2,w2,AngCut)*std::sin(theta_o);
214 
215  //G4cout << "S2: " << S2(costheta_i_squared,klk2) << G4endl;
216  //G4cout << "S2: " << S2(costheta_o_squared,klk2) << G4endl;
217  //G4cout << "Fmu: " << Fmu(k2,theta_i,theta_o,phi_o,b2,w2,AngCut)*sin(theta_o) << G4endl;
218  // checks if the new probability is larger than the
219  // maximum found so far
220 
221  if (Intens>*max)
222  {
223  *max=Intens;
224  max_theta_o=theta_o;
225  max_phi_o=phi_o;
226  }
227 
228  // Adds the newly found probability to the integral probability
229 
230  wkeit+=Intens*ang_steptheta*ang_stepphi;
231  }
232  }
233 
234  // Fine-Iteration to find maximum of distribution
235  // only if the energy is not zero
236 
237  if (E>1e-10*eV)
238  {
239 
240  // Break-condition for refining
241 
242  // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
243  while ((ang_stepphi>=AngCut*AngCut) || (ang_steptheta>=AngCut*AngCut))
244  {
245  a_max_theta_o=max_theta_o;
246  a_max_phi_o=max_phi_o;
247  ang_stepphi /= 2.;
248  ang_steptheta /= 2.;
249 
250  //G4cout << ang_stepphi/degree << ", "
251  // << ang_steptheta/degree << ","
252  // << AngCut/degree << G4endl;
253 
254  for (theta_o=a_max_theta_o-ang_steptheta;
255  theta_o<=a_max_theta_o-ang_steptheta+1e-6;
256  theta_o+=ang_steptheta)
257  {
258  //G4cout << "theta_o: " << theta_o/degree << G4endl;
259  costheta_o_squared=std::cos(theta_o)*std::cos(theta_o);
260  for (phi_o=a_max_phi_o-ang_stepphi;
261  phi_o<=a_max_phi_o+ang_stepphi+1e-6;
262  phi_o+=ang_stepphi)
263  {
264  //G4cout << "phi_o: " << phi_o/degree << G4endl;
265  Intens=kl4d4/costheta_i*S2(costheta_i_squared, klk2)*
266  S2(costheta_o_squared,klk2)*
267  Fmu(k2,theta_i,theta_o,phi_o,b2,w2,AngCut)*std::sin(theta_o);
268  if (Intens>*max)
269  {
270  *max=Intens;
271  max_theta_o=theta_o;
272  max_phi_o=phi_o;
273  }
274  }
275  }
276  }
277  }
278  return wkeit;
279 }
280 
281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
282 
284  G4double theta_i,
285  G4int AngNoTheta,
286  G4int AngNoPhi, G4double b2,
287  G4double w2, G4double* max,
288  G4double AngCut) const
289 {
290  G4double a_max_thetas_o, max_thetas_o = theta_i;
291  G4double a_max_phis_o, max_phis_o = 0.;
292  G4double thetas_o;
293  G4double phis_o;
294  G4double IntensS;
295  G4double ang_steptheta=180.*degree/(AngNoTheta-1);
296  G4double ang_stepphi=180.*degree/(AngNoPhi-1);
297  G4double costheta_i=std::cos(theta_i);
298  G4double costheta_i_squared=costheta_i*costheta_i;
299 
300  *max = 0.;
301  G4double wkeit=0.;
302 
303  if (E < fermipot) return wkeit;
304 
305  //k_l^4/4
307  hbarc_squared*fermipot*fermipot;
308  // (k_l/k)^2
309  G4double klk2=fermipot/E;
310 
311  // (k_l/k')^2
312  G4double klks2=fermipot/(E-fermipot);
313 
314  // k'/k
315  G4double ksdk=std::sqrt((E-fermipot)/E);
316 
317  G4double costhetas_o_squared;
318 
319  // k
320  G4double k=std::sqrt(2*neutron_mass_c2*E/hbarc_squared);
321 
322  // k'
323  G4double kS=ksdk*k;
324 
325  for (thetas_o=0.*degree; thetas_o<=90.*degree+1e-6; thetas_o+=ang_steptheta)
326  {
327  costhetas_o_squared=std::cos(thetas_o)*std::cos(thetas_o);
328 
329  for (phis_o=-180.*degree; phis_o<=180.*degree+1e-6; phis_o+=ang_stepphi)
330  {
331  //cf. Steyerl-paper p. 176, eq. 21
332  if (costhetas_o_squared>=-klks2) {
333 
334  //calculates probability for a certain theta'_o, phi'_o pair
335 
336  G4double thetarefract = thetas_o;
337  if (std::fabs(std::sin(theta_i)/ksdk) <= 1.)
338  thetarefract = std::asin(std::sin(theta_i)/ksdk);
339 
340  IntensS = kl4d4/costheta_i*ksdk*S2(costheta_i_squared, klk2)*
341  SS2(costhetas_o_squared,klks2)*
342  FmuS(k,kS,theta_i,thetas_o,phis_o,b2,w2,AngCut,thetarefract)*
343  std::sin(thetas_o);
344  } else {
345  IntensS=0.;
346  }
347  // checks if the new probability is larger than
348  // the maximum found so far
349  if (IntensS>*max)
350  {
351  *max=IntensS;
352  }
353  wkeit+=IntensS*ang_steptheta*ang_stepphi;
354  }
355  }
356 
357  // Fine-Iteration to find maximum of distribution
358 
359  if (E>1e-10*eV)
360  {
361 
362  // Break-condition for refining
363 
364  while (ang_stepphi>=AngCut*AngCut || ang_steptheta>=AngCut*AngCut)
365  {
366  a_max_thetas_o=max_thetas_o;
367  a_max_phis_o=max_phis_o;
368  ang_stepphi /= 2.;
369  ang_steptheta /= 2.;
370  //G4cout << ang_stepphi/degree << ", " << ang_steptheta/degree
371  // << ", " << AngCut/degree << G4endl;
372  for (thetas_o=a_max_thetas_o-ang_steptheta;
373  thetas_o<=a_max_thetas_o-ang_steptheta+1e-6;
374  thetas_o+=ang_steptheta)
375  {
376  costhetas_o_squared=std::cos(thetas_o)*std::cos(thetas_o);
377  for (phis_o=a_max_phis_o-ang_stepphi;
378  phis_o<=a_max_phis_o+ang_stepphi+1e-6;
379  phis_o+=ang_stepphi)
380  {
381  G4double thetarefract = thetas_o;
382  if (std::fabs(std::sin(theta_i)/ksdk) <= 1.)
383  thetarefract = std::asin(std::sin(theta_i)/ksdk);
384 
385  IntensS=kl4d4/costheta_i*ksdk*S2(costheta_i_squared, klk2)*
386  SS2(costhetas_o_squared,klks2)*
387  FmuS(k,kS,theta_i,thetas_o,phis_o,b2,w2,AngCut,thetarefract)*
388  std::sin(thetas_o);
389  if (IntensS>*max)
390  {
391  *max=IntensS;
392  max_thetas_o=thetas_o;
393  max_phis_o=phis_o;
394  }
395  }
396  }
397  }
398  }
399  return wkeit;
400 }
401 
402 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
403 
405  G4double theta_i,
406  G4double theta_o,
407  G4double phi_o,
408  G4double b, G4double w,
409  G4double AngCut) const
410 {
411  //k_l^4/4
413  hbarc_squared*fermipot*fermipot;
414 
415  // (k_l/k)^2
416  G4double klk2=fermipot/E;
417 
418  G4double costheta_i=std::cos(theta_i);
419  G4double costheta_o=std::cos(theta_o);
420 
421  // eq. 20 on page 176 in the steyerl paper
422 
423  return kl4d4/costheta_i*S2(costheta_i*costheta_i, klk2)*
424  S2(costheta_o*costheta_o,klk2)*
425  Fmu(2*neutron_mass_c2*E/hbarc_squared,theta_i,theta_o,phi_o,b*b,w*w,AngCut)*
426  std::sin(theta_o);
427 }
428 
429 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
430 
432  G4double theta_i,
433  G4double thetas_o,
434  G4double phis_o, G4double b,
435  G4double w, G4double AngCut) const
436 {
437  //k_l^4/4
439  hbarc_squared*fermipot*fermipot;
440  // (k_l/k)^2
441  G4double klk2=fermipot/E;
442 
443  // (k_l/k')^2
444  G4double klks2=fermipot/(E-fermipot);
445 
446  if (E < fermipot) {
447  G4cout << " ProbIminus E < fermipot " << G4endl;
448  return 0.;
449  }
450 
451  // k'/k
452  G4double ksdk=std::sqrt((E-fermipot)/E);
453 
454  // k
455  G4double k=std::sqrt(2*neutron_mass_c2*E/hbarc_squared);
456 
457  // k'/k
458  G4double kS=ksdk*k;
459 
460  G4double costheta_i=std::cos(theta_i);
461  G4double costhetas_o=std::cos(thetas_o);
462 
463  // eq. 20 on page 176 in the steyerl paper
464 
465  G4double thetarefract = thetas_o;
466  if(std::fabs(std::sin(theta_i)/ksdk) <= 1.)thetarefract = std::asin(std::sin(theta_i)/ksdk);
467 
468  return kl4d4/costheta_i*ksdk*S2(costheta_i*costheta_i, klk2)*
469  SS2(costhetas_o*costhetas_o,klks2)*
470  FmuS(k,kS,theta_i,thetas_o,phis_o,b*b,w*w,AngCut,thetarefract)*
471  std::sin(thetas_o);
472 }
473 
474 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4double S2(G4double, G4double) const
G4double ProbIminus(G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double) const
int G4int
Definition: G4Types.hh:78
G4double Fmu(G4double, G4double, G4double, G4double, G4double, G4double, G4double) const
static constexpr double twopi
Definition: G4SIunits.hh:76
tuple b
Definition: test.py:12
G4GLOB_DLL std::ostream G4cout
static constexpr double degree
Definition: G4SIunits.hh:144
G4double ProbIplus(G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double) const
static constexpr double eV
Definition: G4SIunits.hh:215
static G4UCNMicroRoughnessHelper * GetInstance()
float neutron_mass_c2
Definition: hepunit.py:276
G4double IntIplus(G4double, G4double, G4double, G4int, G4int, G4double, G4double, G4double *, G4double) const
G4double IntIminus(G4double, G4double, G4double, G4int, G4int, G4double, G4double, G4double *, G4double) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double SS2(G4double, G4double) const
G4double FmuS(G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double, G4double) const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76