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