| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 << 25 //
>> 26 //
>> 27 //
>> 28 //
>> 29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 30 //
26 // This file contains the source code of vario 31 // This file contains the source code of various functions all related to the
27 // calculation of microroughness. 32 // calculation of microroughness.
28 // 33 //
29 // see A. Steyerl, Z. Physik 254 (1972) 169. 34 // see A. Steyerl, Z. Physik 254 (1972) 169.
30 // 35 //
31 // A description of the functions can be found 36 // A description of the functions can be found in the corresponding header file
>> 37 //
>> 38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 39
33 #include "G4UCNMicroRoughnessHelper.hh" 40 #include "G4UCNMicroRoughnessHelper.hh"
34 41
35 #include "G4PhysicalConstants.hh" <<
36 #include "G4SystemOfUnits.hh" <<
37 #include "globals.hh" 42 #include "globals.hh"
38 43
>> 44 #include "G4SystemOfUnits.hh"
>> 45 #include "G4PhysicalConstants.hh"
>> 46
39 G4UCNMicroRoughnessHelper* G4UCNMicroRoughness 47 G4UCNMicroRoughnessHelper* G4UCNMicroRoughnessHelper::fpInstance = nullptr;
40 48
41 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
42 50
43 G4UCNMicroRoughnessHelper::~G4UCNMicroRoughnes << 51 // Constructor
>> 52
>> 53 G4UCNMicroRoughnessHelper::G4UCNMicroRoughnessHelper() {;}
>> 54
>> 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 56
>> 57 G4UCNMicroRoughnessHelper::~G4UCNMicroRoughnessHelper()
>> 58 {
>> 59 fpInstance = nullptr;
>> 60 }
44 61
45 //....oooOO0OOooo........oooOO0OOooo........oo 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46 63
47 G4UCNMicroRoughnessHelper* G4UCNMicroRoughness 64 G4UCNMicroRoughnessHelper* G4UCNMicroRoughnessHelper::GetInstance()
48 { 65 {
49 if (fpInstance == nullptr) { << 66 if(fpInstance == nullptr)
>> 67 {
50 fpInstance = new G4UCNMicroRoughnessHelper 68 fpInstance = new G4UCNMicroRoughnessHelper;
51 } 69 }
52 return fpInstance; 70 return fpInstance;
53 } 71 }
54 72
55 //....oooOO0OOooo........oooOO0OOooo........oo 73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 74
57 G4double G4UCNMicroRoughnessHelper::S2(G4doubl << 75 G4double
>> 76 G4UCNMicroRoughnessHelper::S2(G4double costheta2, G4double klk2) const
58 { 77 {
59 // case 1: radicand is positive, 78 // case 1: radicand is positive,
60 // case 2: radicand is negative, cf. p. 174 79 // case 2: radicand is negative, cf. p. 174 of the Steyerl paper
61 80
62 if (costheta2 >= klk2) { << 81 if(costheta2 >= klk2)
63 return 4 * costheta2 / (2 * costheta2 - kl << 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))));
64 } 92 }
65 <<
66 return std::norm(2 * std::sqrt(costheta2) / <<
67 (std::sqrt(costheta2) + std <<
68 } 93 }
69 94
70 //....oooOO0OOooo........oooOO0OOooo........oo 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 96
72 G4double G4UCNMicroRoughnessHelper::SS2(G4doub << 97 G4double
>> 98 G4UCNMicroRoughnessHelper::SS2(G4double costhetas2, G4double klks2) const
73 { 99 {
74 // cf. p. 175 of the Steyerl paper 100 // cf. p. 175 of the Steyerl paper
75 101
76 return 4 * costhetas2 / << 102 return 4*costhetas2/
77 (2 * costhetas2 + klks2 + 2 * std::sq << 103 (2*costhetas2+klks2+2*std::sqrt(costhetas2*(costhetas2+klks2)));
78 } 104 }
79 105
80 //....oooOO0OOooo........oooOO0OOooo........oo 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81 107
82 G4double G4UCNMicroRoughnessHelper::Fmu(G4doub << 108 G4double G4UCNMicroRoughnessHelper::Fmu(G4double k2, G4double thetai,
83 G4double phio, G4double b2, G4double w2, G4d << 109 G4double thetao, G4double phio,
>> 110 G4double b2, G4double w2,
>> 111 G4double AngCut) const
84 { 112 {
85 G4double mu_squared; 113 G4double mu_squared;
86 114
87 // Checks if the distribution is peaked arou 115 // Checks if the distribution is peaked around the specular direction
88 116
89 if ((std::fabs(thetai - thetao) < AngCut) && << 117 if((std::fabs(thetai - thetao) < AngCut) && (std::fabs(phio) < AngCut))
90 mu_squared = 0.; << 118 {
>> 119 mu_squared=0.;
91 } 120 }
92 else { << 121 else
>> 122 {
93 // cf. p. 177 of the Steyerl paper 123 // cf. p. 177 of the Steyerl paper
94 124
95 G4double sinthetai = std::sin(thetai); 125 G4double sinthetai = std::sin(thetai);
96 G4double sinthetao = std::sin(thetao); 126 G4double sinthetao = std::sin(thetao);
97 mu_squared = k2 * (sinthetai * sinthetai + << 127 mu_squared = k2 * (sinthetai * sinthetai + sinthetao * sinthetao -
98 2. * sinthetai * sinth << 128 2. * sinthetai * sinthetao * std::cos(phio));
99 } << 129 }
100 130
101 // cf. p. 177 of the Steyerl paper 131 // cf. p. 177 of the Steyerl paper
102 132
103 return b2 * w2 / twopi * std::exp(-mu_square << 133 return b2*w2/twopi*std::exp(-mu_squared*w2/2);
104 } 134 }
105 135
106 //....oooOO0OOooo........oooOO0OOooo........oo 136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
107 137
108 G4double G4UCNMicroRoughnessHelper::FmuS(G4dou << 138 G4double
109 G4double phiSo, G4double b2, G4double w2, G4 << 139 G4UCNMicroRoughnessHelper::FmuS(G4double k, G4double kS,
>> 140 G4double thetai, G4double thetaSo,
>> 141 G4double phiSo,
>> 142 G4double b2, G4double w2,
>> 143 G4double AngCut, G4double thetarefract) const
110 { 144 {
111 G4double mu_squared; 145 G4double mu_squared;
112 146
113 // Checks if the distribution is peaked arou 147 // Checks if the distribution is peaked around the direction of
114 // unperturbed refraction 148 // unperturbed refraction
115 149
116 if ((std::fabs(thetarefract - thetaSo) < Ang << 150 if((std::fabs(thetarefract - thetaSo) < AngCut) &&
117 mu_squared = 0.; << 151 (std::fabs(phiSo) < AngCut))
118 } << 152 {
119 else { << 153 mu_squared=0.;
120 G4double sinthetai = std::sin(thetai); << 154 }
>> 155 else
>> 156 {
>> 157 G4double sinthetai = std::sin(thetai);
121 G4double sinthetaSo = std::sin(thetaSo); 158 G4double sinthetaSo = std::sin(thetaSo);
122 159
123 // cf. p. 177 of the Steyerl paper 160 // cf. p. 177 of the Steyerl paper
124 mu_squared = k * k * sinthetai * sinthetai << 161 mu_squared = k * k * sinthetai * sinthetai +
>> 162 kS * kS * sinthetaSo * sinthetaSo -
125 2. * k * kS * sinthetai * sin 163 2. * k * kS * sinthetai * sinthetaSo * std::cos(phiSo);
126 } << 164 }
127 165
128 // cf. p. 177 of the Steyerl paper 166 // cf. p. 177 of the Steyerl paper
129 167
130 return b2 * w2 / twopi * std::exp(-mu_square << 168 return b2*w2/twopi*std::exp(-mu_squared*w2/2);
131 } 169 }
132 170
133 //....oooOO0OOooo........oooOO0OOooo........oo 171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
134 172
135 G4double G4UCNMicroRoughnessHelper::IntIplus(G << 173 G4double
136 G4int AngNoTheta, G4int AngNoPhi, G4double b << 174 G4UCNMicroRoughnessHelper::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
137 { 179 {
138 *max = 0.; << 180 *max=0.;
139 181
140 // max_theta_o saves the theta-position of t 182 // max_theta_o saves the theta-position of the max probability,
141 // the previous value is saved in a_max_thet 183 // the previous value is saved in a_max_theta_o
142 184
143 G4double a_max_theta_o, max_theta_o = theta_ << 185 G4double a_max_theta_o, max_theta_o=theta_i, a_max_phi_o, max_phi_o=0.;
144 186
145 // max_phi_o saves the phi-position of the m 187 // max_phi_o saves the phi-position of the max probability,
146 // the previous value is saved in a_max_phi_ 188 // the previous value is saved in a_max_phi_o
147 189
148 // Definition of the stepsizes in theta_o an 190 // Definition of the stepsizes in theta_o and phi_o
149 191
150 G4double theta_o; 192 G4double theta_o;
151 G4double phi_o; 193 G4double phi_o;
152 G4double Intens; 194 G4double Intens;
153 G4double ang_steptheta = 90. * degree / (Ang << 195 G4double ang_steptheta=90.*degree/(AngNoTheta-1);
154 G4double ang_stepphi = 360. * degree / (AngN << 196 G4double ang_stepphi=360.*degree/(AngNoPhi-1);
155 G4double costheta_i = std::cos(theta_i); << 197 G4double costheta_i=std::cos(theta_i);
156 G4double costheta_i_squared = costheta_i * c << 198 G4double costheta_i_squared=costheta_i*costheta_i;
157 199
158 // (k_l/k)^2 200 // (k_l/k)^2
159 G4double kl4d4 = << 201 G4double kl4d4=neutron_mass_c2/hbarc_squared*neutron_mass_c2/
160 neutron_mass_c2 / hbarc_squared * neutron_ << 202 hbarc_squared*fermipot*fermipot;
161 203
162 // (k_l/k)^2 204 // (k_l/k)^2
163 G4double klk2 = fermipot / E; << 205 G4double klk2=fermipot/E;
164 206
165 G4double costheta_o_squared; 207 G4double costheta_o_squared;
166 208
167 // k^2 209 // k^2
168 G4double k2 = 2 * neutron_mass_c2 * E / hbar << 210 G4double k2=2*neutron_mass_c2*E/hbarc_squared;
169 211
170 G4double wkeit = 0.; << 212 G4double wkeit=0.;
171 213
172 // Loop through theta_o 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 }
173 242
174 for (theta_o = 0. * degree; theta_o <= 90. * << 243 // Adds the newly found probability to the integral probability
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 244
191 // Adds the newly found probability to t << 245 wkeit+=Intens*ang_steptheta*ang_stepphi;
192 << 246 }
193 wkeit += Intens * ang_steptheta * ang_st <<
194 } 247 }
195 } <<
196 248
197 // Fine-Iteration to find maximum of distrib 249 // Fine-Iteration to find maximum of distribution
198 // only if the energy is not zero 250 // only if the energy is not zero
199 251
200 if (E > 1e-10 * eV) { << 252 if (E>1e-10*eV)
201 // Break-condition for refining << 253 {
>> 254
>> 255 // Break-condition for refining
202 256
203 // Loop checking, 07-Aug-2015, Vladimir Iv << 257 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
204 while ((ang_stepphi >= AngCut * AngCut) || << 258 while ((ang_stepphi>=AngCut*AngCut) || (ang_steptheta>=AngCut*AngCut))
205 a_max_theta_o = max_theta_o; << 259 {
206 a_max_phi_o = max_phi_o; << 260 a_max_theta_o=max_theta_o;
>> 261 a_max_phi_o=max_phi_o;
207 ang_stepphi /= 2.; 262 ang_stepphi /= 2.;
208 ang_steptheta /= 2.; 263 ang_steptheta /= 2.;
209 264
210 for (theta_o = a_max_theta_o - ang_stept << 265 //G4cout << ang_stepphi/degree << ", "
211 theta_o += ang_steptheta) << 266 // << ang_steptheta/degree << ","
212 { << 267 // << AngCut/degree << G4endl;
213 // G4cout << "theta_o: " << theta_o/de << 268
214 costheta_o_squared = std::cos(theta_o) << 269 for (theta_o=a_max_theta_o-ang_steptheta;
215 for (phi_o = a_max_phi_o - ang_stepphi << 270 theta_o<=a_max_theta_o-ang_steptheta+1e-6;
216 phi_o += ang_stepphi) << 271 theta_o+=ang_steptheta)
217 { << 272 {
218 // G4cout << "phi_o: " << phi_o/degr << 273 //G4cout << "theta_o: " << theta_o/degree << G4endl;
219 Intens = kl4d4 / costheta_i * S2(cos << 274 costheta_o_squared=std::cos(theta_o)*std::cos(theta_o);
220 S2(costheta_o_squared, klk2 << 275 for (phi_o=a_max_phi_o-ang_stepphi;
221 std::sin(theta_o); << 276 phi_o<=a_max_phi_o+ang_stepphi+1e-6;
222 if (Intens > *max) { << 277 phi_o+=ang_stepphi)
223 *max = Intens; << 278 {
224 max_theta_o = theta_o; << 279 //G4cout << "phi_o: " << phi_o/degree << G4endl;
225 max_phi_o = phi_o; << 280 Intens=kl4d4/costheta_i*S2(costheta_i_squared, klk2)*
226 } << 281 S2(costheta_o_squared,klk2)*
227 } << 282 Fmu(k2,theta_i,theta_o,phi_o,b2,w2,AngCut)*std::sin(theta_o);
228 } << 283 if (Intens>*max)
>> 284 {
>> 285 *max=Intens;
>> 286 max_theta_o=theta_o;
>> 287 max_phi_o=phi_o;
>> 288 }
>> 289 }
>> 290 }
>> 291 }
229 } 292 }
230 } <<
231 return wkeit; 293 return wkeit;
232 } 294 }
233 295
234 //....oooOO0OOooo........oooOO0OOooo........oo 296 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
235 297
236 G4double G4UCNMicroRoughnessHelper::IntIminus( << 298 G4double G4UCNMicroRoughnessHelper::IntIminus(G4double E, G4double fermipot,
237 G4int AngNoTheta, G4int AngNoPhi, G4double b << 299 G4double theta_i,
>> 300 G4int AngNoTheta,
>> 301 G4int AngNoPhi, G4double b2,
>> 302 G4double w2, G4double* max,
>> 303 G4double AngCut) const
238 { 304 {
239 G4double a_max_thetas_o, max_thetas_o = thet 305 G4double a_max_thetas_o, max_thetas_o = theta_i;
240 G4double a_max_phis_o, max_phis_o = 0.; 306 G4double a_max_phis_o, max_phis_o = 0.;
241 G4double thetas_o; 307 G4double thetas_o;
242 G4double phis_o; 308 G4double phis_o;
243 G4double IntensS; 309 G4double IntensS;
244 G4double ang_steptheta = 180. * degree / (An << 310 G4double ang_steptheta=180.*degree/(AngNoTheta-1);
245 G4double ang_stepphi = 180. * degree / (AngN << 311 G4double ang_stepphi=180.*degree/(AngNoPhi-1);
246 G4double costheta_i = std::cos(theta_i); << 312 G4double costheta_i=std::cos(theta_i);
247 G4double costheta_i_squared = costheta_i * c << 313 G4double costheta_i_squared=costheta_i*costheta_i;
248 314
249 *max = 0.; 315 *max = 0.;
250 G4double wkeit = 0.; << 316 G4double wkeit=0.;
251 317
252 if (E < fermipot) { << 318 if(E < fermipot)
>> 319 {
253 return wkeit; 320 return wkeit;
254 } 321 }
255 322
256 // k_l^4/4 << 323 //k_l^4/4
257 G4double kl4d4 = << 324 G4double kl4d4=neutron_mass_c2/hbarc_squared*neutron_mass_c2/
258 neutron_mass_c2 / hbarc_squared * neutron_ << 325 hbarc_squared*fermipot*fermipot;
259 // (k_l/k)^2 326 // (k_l/k)^2
260 G4double klk2 = fermipot / E; << 327 G4double klk2=fermipot/E;
261 328
262 // (k_l/k')^2 329 // (k_l/k')^2
263 G4double klks2 = fermipot / (E - fermipot); << 330 G4double klks2=fermipot/(E-fermipot);
264 331
265 // k'/k 332 // k'/k
266 G4double ksdk = std::sqrt((E - fermipot) / E << 333 G4double ksdk=std::sqrt((E-fermipot)/E);
267 334
268 G4double costhetas_o_squared; 335 G4double costhetas_o_squared;
269 336
270 // k 337 // k
271 G4double k = std::sqrt(2 * neutron_mass_c2 * << 338 G4double k=std::sqrt(2*neutron_mass_c2*E/hbarc_squared);
272 339
273 // k' 340 // k'
274 G4double kS = ksdk * k; << 341 G4double kS=ksdk*k;
275 <<
276 for (thetas_o = 0. * degree; thetas_o <= 90. <<
277 costhetas_o_squared = std::cos(thetas_o) * <<
278 342
279 for (phis_o = -180. * degree; phis_o <= 18 << 343 for (thetas_o=0.*degree; thetas_o<=90.*degree+1e-6; thetas_o+=ang_steptheta)
280 // cf. Steyerl-paper p. 176, eq. 21 << 344 {
281 if (costhetas_o_squared >= -klks2) { << 345 costhetas_o_squared=std::cos(thetas_o)*std::cos(thetas_o);
282 // calculates probability for a certai << 346
283 << 347 for (phis_o=-180.*degree; phis_o<=180.*degree+1e-6; phis_o+=ang_stepphi)
284 G4double thetarefract = thetas_o; << 348 {
285 if (std::fabs(std::sin(theta_i) / ksdk << 349 //cf. Steyerl-paper p. 176, eq. 21
286 thetarefract = std::asin(std::sin(th << 350 if (costhetas_o_squared>=-klks2) {
287 } << 351
288 << 352 //calculates probability for a certain theta'_o, phi'_o pair
289 IntensS = kl4d4 / costheta_i * ksdk * << 353
290 SS2(costhetas_o_squared, klk << 354 G4double thetarefract = thetas_o;
291 FmuS(k, kS, theta_i, thetas_ << 355 if(std::fabs(std::sin(theta_i) / ksdk) <= 1.)
292 std::sin(thetas_o); << 356 {
293 } << 357 thetarefract = std::asin(std::sin(theta_i) / ksdk);
294 else { << 358 }
295 IntensS = 0.; << 359
296 } << 360 IntensS = kl4d4/costheta_i*ksdk*S2(costheta_i_squared, klk2)*
297 // checks if the new probability is larg << 361 SS2(costhetas_o_squared,klks2)*
298 // the maximum found so far << 362 FmuS(k,kS,theta_i,thetas_o,phis_o,b2,w2,AngCut,thetarefract)*
299 if (IntensS > *max) { << 363 std::sin(thetas_o);
300 *max = IntensS; << 364 } else {
301 } << 365 IntensS=0.;
302 wkeit += IntensS * ang_steptheta * ang_s << 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 }
303 } 375 }
304 } <<
305 376
306 // Fine-Iteration to find maximum of distrib 377 // Fine-Iteration to find maximum of distribution
307 378
308 if (E > 1e-10 * eV) { << 379 if (E>1e-10*eV)
309 // Break-condition for refining << 380 {
>> 381
>> 382 // Break-condition for refining
310 383
311 while (ang_stepphi >= AngCut * AngCut || a << 384 while (ang_stepphi>=AngCut*AngCut || ang_steptheta>=AngCut*AngCut)
312 a_max_thetas_o = max_thetas_o; << 385 {
313 a_max_phis_o = max_phis_o; << 386 a_max_thetas_o=max_thetas_o;
>> 387 a_max_phis_o=max_phis_o;
314 ang_stepphi /= 2.; 388 ang_stepphi /= 2.;
315 ang_steptheta /= 2.; 389 ang_steptheta /= 2.;
316 << 390 //G4cout << ang_stepphi/degree << ", " << ang_steptheta/degree
317 for (thetas_o = a_max_thetas_o - ang_ste << 391 // << ", " << AngCut/degree << G4endl;
318 thetas_o <= a_max_thetas_o - ang_st << 392 for (thetas_o=a_max_thetas_o-ang_steptheta;
319 { << 393 thetas_o<=a_max_thetas_o-ang_steptheta+1e-6;
320 costhetas_o_squared = std::cos(thetas_ << 394 thetas_o+=ang_steptheta)
321 for (phis_o = a_max_phis_o - ang_stepp << 395 {
322 phis_o += ang_stepphi) << 396 costhetas_o_squared=std::cos(thetas_o)*std::cos(thetas_o);
323 { << 397 for (phis_o=a_max_phis_o-ang_stepphi;
324 G4double thetarefract = thetas_o; << 398 phis_o<=a_max_phis_o+ang_stepphi+1e-6;
325 if (std::fabs(std::sin(theta_i) / ks << 399 phis_o+=ang_stepphi)
326 thetarefract = std::asin(std::sin( << 400 {
327 } << 401 G4double thetarefract = thetas_o;
328 << 402 if(std::fabs(std::sin(theta_i) / ksdk) <= 1.)
329 IntensS = kl4d4 / costheta_i * ksdk << 403 {
330 SS2(costhetas_o_squared, k << 404 thetarefract = std::asin(std::sin(theta_i) / ksdk);
331 FmuS(k, kS, theta_i, theta << 405 }
332 std::sin(thetas_o); << 406
333 if (IntensS > *max) { << 407 IntensS=kl4d4/costheta_i*ksdk*S2(costheta_i_squared, klk2)*
334 *max = IntensS; << 408 SS2(costhetas_o_squared,klks2)*
335 max_thetas_o = thetas_o; << 409 FmuS(k,kS,theta_i,thetas_o,phis_o,b2,w2,AngCut,thetarefract)*
336 max_phis_o = phis_o; << 410 std::sin(thetas_o);
337 } << 411 if (IntensS>*max)
338 } << 412 {
339 } << 413 *max=IntensS;
>> 414 max_thetas_o=thetas_o;
>> 415 max_phis_o=phis_o;
>> 416 }
>> 417 }
>> 418 }
>> 419 }
340 } 420 }
341 } <<
342 return wkeit; 421 return wkeit;
343 } 422 }
344 423
345 //....oooOO0OOooo........oooOO0OOooo........oo 424 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
346 425
347 G4double G4UCNMicroRoughnessHelper::ProbIplus( << 426 G4double G4UCNMicroRoughnessHelper::ProbIplus(G4double E, G4double fermipot,
348 G4double theta_o, G4double phi_o, G4double b << 427 G4double theta_i,
>> 428 G4double theta_o,
>> 429 G4double phi_o,
>> 430 G4double b, G4double w,
>> 431 G4double AngCut) const
349 { 432 {
350 // k_l^4/4 << 433 //k_l^4/4
351 G4double kl4d4 = << 434 G4double kl4d4=neutron_mass_c2/hbarc_squared*neutron_mass_c2/
352 neutron_mass_c2 / hbarc_squared * neutron_ << 435 hbarc_squared*fermipot*fermipot;
353 436
354 // (k_l/k)^2 437 // (k_l/k)^2
355 G4double klk2 = fermipot / E; << 438 G4double klk2=fermipot/E;
356 439
357 G4double costheta_i = std::cos(theta_i); << 440 G4double costheta_i=std::cos(theta_i);
358 G4double costheta_o = std::cos(theta_o); << 441 G4double costheta_o=std::cos(theta_o);
359 442
360 // eq. 20 on page 176 in the steyerl paper 443 // eq. 20 on page 176 in the steyerl paper
361 444
362 return kl4d4 / costheta_i * S2(costheta_i * << 445 return kl4d4/costheta_i*S2(costheta_i*costheta_i, klk2)*
363 S2(costheta_o * costheta_o, klk2) * << 446 S2(costheta_o*costheta_o,klk2)*
364 Fmu( << 447 Fmu(2*neutron_mass_c2*E/hbarc_squared,theta_i,theta_o,phi_o,b*b,w*w,AngCut)*
365 2 * neutron_mass_c2 * E / hbarc_squ << 448 std::sin(theta_o);
366 std::sin(theta_o); <<
367 } 449 }
368 450
369 //....oooOO0OOooo........oooOO0OOooo........oo 451 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
370 452
371 G4double G4UCNMicroRoughnessHelper::ProbIminus << 453 G4double G4UCNMicroRoughnessHelper::ProbIminus(G4double E, G4double fermipot,
372 G4double thetas_o, G4double phis_o, G4double << 454 G4double theta_i,
>> 455 G4double thetas_o,
>> 456 G4double phis_o, G4double b,
>> 457 G4double w, G4double AngCut) const
373 { 458 {
374 // k_l^4/4 << 459 //k_l^4/4
375 G4double kl4d4 = << 460 G4double kl4d4=neutron_mass_c2/hbarc_squared*neutron_mass_c2/
376 neutron_mass_c2 / hbarc_squared * neutron_ << 461 hbarc_squared*fermipot*fermipot;
377 // (k_l/k)^2 462 // (k_l/k)^2
378 G4double klk2 = fermipot / E; << 463 G4double klk2=fermipot/E;
379 464
380 // (k_l/k')^2 465 // (k_l/k')^2
381 G4double klks2 = fermipot / (E - fermipot); << 466 G4double klks2=fermipot/(E-fermipot);
382 467
383 if (E < fermipot) { 468 if (E < fermipot) {
384 G4cout << " ProbIminus E < fermipot " << G << 469 G4cout << " ProbIminus E < fermipot " << G4endl;
385 return 0.; << 470 return 0.;
386 } 471 }
387 472
388 // k'/k 473 // k'/k
389 G4double ksdk = std::sqrt((E - fermipot) / E << 474 G4double ksdk=std::sqrt((E-fermipot)/E);
390 475
391 // k 476 // k
392 G4double k = std::sqrt(2 * neutron_mass_c2 * << 477 G4double k=std::sqrt(2*neutron_mass_c2*E/hbarc_squared);
393 478
394 // k'/k 479 // k'/k
395 G4double kS = ksdk * k; << 480 G4double kS=ksdk*k;
396 481
397 G4double costheta_i = std::cos(theta_i); << 482 G4double costheta_i=std::cos(theta_i);
398 G4double costhetas_o = std::cos(thetas_o); << 483 G4double costhetas_o=std::cos(thetas_o);
399 484
400 // eq. 20 on page 176 in the steyerl paper 485 // eq. 20 on page 176 in the steyerl paper
401 486
402 G4double thetarefract = thetas_o; 487 G4double thetarefract = thetas_o;
403 if (std::fabs(std::sin(theta_i) / ksdk) <= 1 << 488 if(std::fabs(std::sin(theta_i) / ksdk) <= 1.)
>> 489 {
404 thetarefract = std::asin(std::sin(theta_i) 490 thetarefract = std::asin(std::sin(theta_i) / ksdk);
405 } 491 }
406 492
407 return kl4d4 / costheta_i * ksdk * S2(costhe << 493 return kl4d4/costheta_i*ksdk*S2(costheta_i*costheta_i, klk2)*
408 SS2(costhetas_o * costhetas_o, klks2) << 494 SS2(costhetas_o*costhetas_o,klks2)*
409 FmuS(k, kS, theta_i, thetas_o, phis_o << 495 FmuS(k,kS,theta_i,thetas_o,phis_o,b*b,w*w,AngCut,thetarefract)*
410 std::sin(thetas_o); 496 std::sin(thetas_o);
411 } 497 }
>> 498
>> 499 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
412 500