ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4UCNMicroRoughnessHelper.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4UCNMicroRoughnessHelper.cc
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 
47 G4UCNMicroRoughnessHelper* G4UCNMicroRoughnessHelper::fpInstance = nullptr;
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 
51 // Constructor
52 
53 G4UCNMicroRoughnessHelper::G4UCNMicroRoughnessHelper() {;}
54 
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 
57 G4UCNMicroRoughnessHelper::~G4UCNMicroRoughnessHelper()
58 {
59  delete fpInstance;
60  fpInstance = nullptr;
61 }
62 
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
64 
65 G4UCNMicroRoughnessHelper* G4UCNMicroRoughnessHelper::GetInstance()
66 {
67  if (fpInstance == nullptr) fpInstance = new G4UCNMicroRoughnessHelper;
68  return fpInstance;
69 }
70 
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72 
73 G4double
74 G4UCNMicroRoughnessHelper::S2(G4double costheta2, G4double klk2) const
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
88 G4UCNMicroRoughnessHelper::SS2(G4double costhetas2, G4double klks2) const
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 
98 G4double G4UCNMicroRoughnessHelper::Fmu(G4double k2, G4double thetai,
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
128 G4UCNMicroRoughnessHelper::FmuS(G4double k, G4double kS,
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
159 G4UCNMicroRoughnessHelper::IntIplus(G4double E, G4double fermipot,
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
186  G4double kl4d4=neutron_mass_c2/hbarc_squared*neutron_mass_c2/
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
195  G4double k2=2*neutron_mass_c2*E/hbarc_squared;
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 
283 G4double G4UCNMicroRoughnessHelper::IntIminus(G4double E, G4double fermipot,
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
306  G4double kl4d4=neutron_mass_c2/hbarc_squared*neutron_mass_c2/
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 
404 G4double G4UCNMicroRoughnessHelper::ProbIplus(G4double E, G4double fermipot,
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
412  G4double kl4d4=neutron_mass_c2/hbarc_squared*neutron_mass_c2/
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 
431 G4double G4UCNMicroRoughnessHelper::ProbIminus(G4double E, G4double fermipot,
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
438  G4double kl4d4=neutron_mass_c2/hbarc_squared*neutron_mass_c2/
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......