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