[819] | 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 | // GEANT4 Class file |
---|
| 30 | // |
---|
| 31 | // |
---|
| 32 | // File name: G4Generator2BN |
---|
| 33 | // |
---|
| 34 | // Author: Andreia Trindade (andreia@lip.pt) |
---|
| 35 | // Pedro Rodrigues (psilva@lip.pt) |
---|
| 36 | // Luis Peralta (luis@lip.pt) |
---|
| 37 | // |
---|
| 38 | // Creation date: 21 June 2003 |
---|
| 39 | // |
---|
| 40 | // Modifications: |
---|
| 41 | // 21 Jun 2003 First implementation acording with new design |
---|
| 42 | // 05 Nov 2003 MGP Fixed compilation warning |
---|
| 43 | // 14 Mar 2004 Code optimization |
---|
| 44 | // |
---|
| 45 | // Class Description: |
---|
| 46 | // |
---|
| 47 | // Concrete base class for Bremsstrahlung Angular Distribution Generation - 2BN Distribution |
---|
| 48 | // |
---|
| 49 | // Class Description: End |
---|
| 50 | // |
---|
| 51 | // ------------------------------------------------------------------- |
---|
| 52 | // |
---|
| 53 | // |
---|
| 54 | |
---|
| 55 | #include "G4Generator2BN.hh" |
---|
| 56 | #include "Randomize.hh" |
---|
| 57 | // |
---|
| 58 | |
---|
| 59 | G4double G4Generator2BN::Atab[320] = |
---|
| 60 | { 0.0264767, 0.0260006, 0.0257112, 0.0252475, 0.024792, 0.0243443, 0.0240726, 0.0236367, |
---|
| 61 | 0.023209, 0.0227892, 0.0225362, 0.0221284, 0.0217282,0.0214894, 0.0211005, 0.0207189, |
---|
| 62 | 0.0204935, 0.0201227, 0.0197588, 0.019546, 0.0191923,0.0188454, 0.0186445, 0.0183072, |
---|
| 63 | 0.0179763, 0.0177866, 0.0174649, 0.0172828, 0.0169702,0.0166634, 0.0164915, 0.0161933, |
---|
| 64 | 0.0160283, 0.0157384, 0.0155801, 0.0152981, 0.0151463,0.0148721, 0.0147263, 0.0144598, |
---|
| 65 | 0.01432, 0.0140607, 0.0139267, 0.0136744, 0.0135459,0.0133005, 0.0131773, 0.0129385, |
---|
| 66 | 0.0128205, 0.0125881, 0.012475, 0.0122488, 0.0121406,0.0119204, 0.0118167, 0.0117158, |
---|
| 67 | 0.0115032, 0.0114067, 0.0111995, 0.0111072, 0.0110175,0.0108173, 0.0107316, 0.0105365, |
---|
| 68 | 0.0104547, 0.0102646, 0.0101865, 0.010111, 0.00992684,0.0098548,0.00967532,0.00960671, |
---|
| 69 | 0.00943171,0.00936643,0.00930328,0.0091337, 0.00907372,0.00890831,0.00885141,0.00869003, |
---|
| 70 | 0.00863611,0.00858428,0.00842757,0.00837854,0.0082256,0.00817931,0.00803,0.00798639, |
---|
| 71 | 0.00784058,0.00779958,0.00776046,0.00761866,0.00758201,0.00744346,0.00740928,0.00727384, |
---|
| 72 | 0.00724201,0.00710969,0.00708004,0.00695074,0.00692333,0.00679688,0.00677166,0.00664801, |
---|
| 73 | 0.00662484,0.00650396,0.00648286,0.00636458,0.00634545,0.00622977,0.00621258,0.00609936, |
---|
| 74 | 0.00608412,0.00597331,0.00595991,0.00585143,0.00583988,0.0057337,0.0057239,0.00561991, |
---|
| 75 | 0.0056119, 0.00551005,0.00550377,0.00540399,0.00539938,0.00530162,0.00529872,0.00520292, |
---|
| 76 | 0.0051091, 0.00510777,0.00501582,0.00501608,0.00492594,0.00492781,0.00483942,0.0048429, |
---|
| 77 | 0.00475622,0.00476127,0.00467625,0.00468287,0.00459947,0.00451785,0.00452581,0.00444573, |
---|
| 78 | 0.00445522,0.00437664,0.00438768,0.00431057,0.00432316,0.00424745,0.0042616,0.00418726, |
---|
| 79 | 0.004203, 0.00413, 0.00405869,0.00407563,0.00400561,0.00402414,0.00395536,0.00397553, |
---|
| 80 | 0.00390795,0.00392975,0.00386339,0.00379862,0.00382167,0.00375805,0.00378276,0.00372031, |
---|
| 81 | 0.00374678,0.00368538,0.00371363,0.00365335,0.00359463,0.0036242, 0.00356653,0.003598, |
---|
| 82 | 0.00354139,0.00357481,0.00351921,0.00355464,0.00350005,0.0034471,0.00348403,0.00343208, |
---|
| 83 | 0.0034712, 0.00342026,0.00346165,0.00341172,0.00345548,0.00340657,0.00335944,0.00340491, |
---|
| 84 | 0.00335885,0.00340692,0.00336191,0.00341273,0.00336879,0.00342249,0.00337962,0.00333889, |
---|
| 85 | 0.00339463,0.00335506,0.00341401,0.00337558,0.00343797,0.00340067,0.00336584,0.00343059, |
---|
| 86 | 0.0033969, 0.00346557,0.00343302,0.00350594,0.00347448,0.00344563,0.00352163,0.00349383, |
---|
| 87 | 0.00357485,0.00354807,0.00352395,0.00360885,0.00358571,0.00367661,0.00365446,0.00375194, |
---|
| 88 | 0.00373078,0.00371234,0.00381532,0.00379787,0.00390882,0.00389241,0.00387881,0.00399675, |
---|
| 89 | 0.00398425,0.00411183,0.00410042,0.00409197,0.00422843,0.00422123,0.00436974,0.0043637, |
---|
| 90 | 0.00436082,0.00452075,0.00451934,0.00452125,0.00469406,0.00469756,0.00488741,0.00489221, |
---|
| 91 | 0.00490102,0.00510782,0.00511801,0.00513271,0.0053589, 0.00537524,0.00562643,0.00564452, |
---|
| 92 | 0.0056677, 0.00594482,0.00596999,0.0059999, 0.00630758,0.00634014,0.00637849,0.00672136, |
---|
| 93 | 0.00676236,0.00680914,0.00719407,0.0072439, 0.00730063,0.0077349, 0.00779494,0.00786293, |
---|
| 94 | 0.00835577,0.0084276, 0.00850759,0.00907162,0.00915592,0.00924925,0.00935226,0.00999779, |
---|
| 95 | 0.0101059, 0.0102249, 0.0109763, 0.0111003, 0.0112367, 0.0113862, 0.0122637, 0.0124196, |
---|
| 96 | 0.0125898, 0.0136311, 0.0138081, 0.0140011, 0.0142112, 0.0154536, 0.0156723, 0.0159099, |
---|
| 97 | 0.016168, 0.0176664, 0.0179339, 0.0182246, 0.0185396, 0.020381, 0.0207026, 0.0210558, |
---|
| 98 | 0.0214374, 0.0237377, 0.0241275, 0.0245528, 0.0250106, 0.0255038, 0.0284158, 0.0289213, |
---|
| 99 | 0.0294621, 0.0300526, 0.0338619, 0.0344537, 0.0351108, 0.0358099, 0.036554, 0.0416399 |
---|
| 100 | }; |
---|
| 101 | |
---|
| 102 | G4double G4Generator2BN::ctab[320] = |
---|
| 103 | { 0.482253, 0.482253, 0.489021, 0.489021, 0.489021, 0.489021, 0.495933, |
---|
| 104 | 0.495933, 0.495933, 0.495933, 0.502993, 0.502993, 0.502993, 0.510204, |
---|
| 105 | 0.510204, 0.510204, 0.517572, 0.517572, 0.517572, 0.5251, 0.5251, |
---|
| 106 | 0.5251, 0.532793, 0.532793, 0.532793, 0.540657, 0.540657, 0.548697, |
---|
| 107 | 0.548697, 0.548697, 0.556917, 0.556917, 0.565323, 0.565323, 0.573921, |
---|
| 108 | 0.573921, 0.582717, 0.582717, 0.591716, 0.591716, 0.600925, 0.600925, |
---|
| 109 | 0.610352, 0.610352, 0.620001, 0.620001, 0.629882, 0.629882, 0.64, |
---|
| 110 | 0.64, 0.650364, 0.650364, 0.660982, 0.660982, 0.671862, 0.683013, |
---|
| 111 | 0.683013, 0.694444, 0.694444, 0.706165, 0.718184, 0.718184, 0.730514, |
---|
| 112 | 0.730514, 0.743163, 0.743163, 0.756144, 0.769468, 0.769468, 0.783147, |
---|
| 113 | 0.783147, 0.797194, 0.797194, 0.811622, 0.826446, 0.826446, 0.84168, |
---|
| 114 | 0.84168, 0.857339, 0.857339, 0.873439, 0.889996, 0.889996, 0.907029, |
---|
| 115 | 0.907029, 0.924556, 0.924556, 0.942596, 0.942596, 0.961169, 0.980296, |
---|
| 116 | 0.980296, 1, 1, 1.0203, 1.0203, 1.04123, 1.04123, |
---|
| 117 | 1.06281, 1.06281, 1.08507, 1.08507, 1.10803, 1.10803, 1.13173, |
---|
| 118 | 1.13173, 1.1562, 1.1562, 1.18147, 1.18147, 1.20758, 1.20758, |
---|
| 119 | 1.23457, 1.23457, 1.26247, 1.26247, 1.29132, 1.29132, 1.32118, |
---|
| 120 | 1.32118, 1.35208, 1.35208, 1.38408, 1.38408, 1.41723, 1.41723, |
---|
| 121 | 1.45159, 1.45159, 1.45159, 1.48721, 1.48721, 1.52416, 1.52416, |
---|
| 122 | 1.5625, 1.5625, 1.60231, 1.60231, 1.64366, 1.64366, 1.68663, |
---|
| 123 | 1.68663, 1.68663, 1.7313, 1.7313, 1.77778, 1.77778, 1.82615, |
---|
| 124 | 1.82615, 1.87652, 1.87652, 1.92901, 1.92901, 1.98373, 1.98373, |
---|
| 125 | 1.98373, 2.04082, 2.04082, 2.1004, 2.1004, 2.16263, 2.16263, |
---|
| 126 | 2.22767, 2.22767, 2.22767, 2.29568, 2.29568, 2.36686, 2.36686, |
---|
| 127 | 2.44141, 2.44141, 2.51953, 2.51953, 2.51953, 2.60146, 2.60146, |
---|
| 128 | 2.68745, 2.68745, 2.77778, 2.77778, 2.87274, 2.87274, 2.87274, |
---|
| 129 | 2.97265, 2.97265, 3.07787, 3.07787, 3.18878, 3.18878, 3.30579, |
---|
| 130 | 3.30579, 3.30579, 3.42936, 3.42936, 3.55999, 3.55999, 3.69822, |
---|
| 131 | 3.69822, 3.84468, 3.84468, 3.84468, 4, 4, 4.16493, |
---|
| 132 | 4.16493, 4.34028, 4.34028, 4.34028, 4.52694, 4.52694, 4.7259, |
---|
| 133 | 4.7259, 4.93827, 4.93827, 4.93827, 5.16529, 5.16529, 5.40833, |
---|
| 134 | 5.40833, 5.40833, 5.66893, 5.66893, 5.94884, 5.94884, 6.25, |
---|
| 135 | 6.25, 6.25, 6.57462, 6.57462, 6.92521, 6.92521, 6.92521, |
---|
| 136 | 7.3046, 7.3046, 7.71605, 7.71605, 7.71605, 8.16327, 8.16327, |
---|
| 137 | 8.65052, 8.65052, 8.65052, 9.18274, 9.18274, 9.18274, 9.76562, |
---|
| 138 | 9.76562, 10.4058, 10.4058, 10.4058, 11.1111, 11.1111, 11.1111, |
---|
| 139 | 11.8906, 11.8906, 12.7551, 12.7551, 12.7551, 13.7174, 13.7174, |
---|
| 140 | 13.7174, 14.7929, 14.7929, 14.7929, 16, 16, 16, |
---|
| 141 | 17.3611, 17.3611, 17.3611, 18.9036, 18.9036, 18.9036, 20.6612, |
---|
| 142 | 20.6612, 20.6612, 22.6757, 22.6757, 22.6757, 22.6757, 25, |
---|
| 143 | 25, 25, 27.7008, 27.7008, 27.7008, 27.7008, 30.8642, |
---|
| 144 | 30.8642, 30.8642, 34.6021, 34.6021, 34.6021, 34.6021, 39.0625, |
---|
| 145 | 39.0625, 39.0625, 39.0625, 44.4444, 44.4444, 44.4444, 44.4444, |
---|
| 146 | 51.0204, 51.0204, 51.0204, 51.0204, 59.1716, 59.1716, 59.1716, |
---|
| 147 | 59.1716, 59.1716, 69.4444, 69.4444, 69.4444, 69.4444, 82.6446, |
---|
| 148 | 82.6446, 82.6446, 82.6446, 82.6446, 100 |
---|
| 149 | }; |
---|
| 150 | |
---|
| 151 | |
---|
| 152 | G4Generator2BN::G4Generator2BN(const G4String& name):G4VBremAngularDistribution(name) |
---|
| 153 | { |
---|
| 154 | b = 1.2; |
---|
| 155 | index_min = -300; |
---|
| 156 | index_max = 20; |
---|
| 157 | |
---|
| 158 | // Set parameters minimum limits Ekelectron = 250 eV and kphoton = 100 eV |
---|
| 159 | kmin = 100*eV; |
---|
| 160 | Ekmin = 250*eV; |
---|
| 161 | kcut = 100*eV; |
---|
| 162 | |
---|
| 163 | // Increment Theta value for surface interpolation |
---|
| 164 | dtheta = 0.1*rad; |
---|
| 165 | |
---|
| 166 | // Construct Majorant Surface to 2BN cross-section |
---|
| 167 | // (we dont need this. Pre-calculated values are used instead due to performance issues |
---|
| 168 | // ConstructMajorantSurface(); |
---|
| 169 | } |
---|
| 170 | |
---|
| 171 | // |
---|
| 172 | |
---|
| 173 | G4Generator2BN::~G4Generator2BN() |
---|
| 174 | {;} |
---|
| 175 | |
---|
| 176 | // |
---|
| 177 | |
---|
| 178 | G4double G4Generator2BN::PolarAngle(const G4double initial_energy, |
---|
| 179 | const G4double final_energy, |
---|
| 180 | const G4int) // Z |
---|
| 181 | { |
---|
| 182 | |
---|
| 183 | G4double theta = 0; |
---|
| 184 | |
---|
| 185 | G4double k = initial_energy - final_energy; |
---|
| 186 | |
---|
| 187 | theta = Generate2BN(initial_energy, k); |
---|
| 188 | |
---|
| 189 | return theta; |
---|
| 190 | } |
---|
| 191 | // |
---|
| 192 | |
---|
| 193 | G4double G4Generator2BN::CalculateFkt(G4double k, G4double theta, G4double A, G4double c) const |
---|
| 194 | { |
---|
| 195 | G4double Fkt_value = 0; |
---|
| 196 | Fkt_value = A*std::pow(k,-b)*theta/(1+c*theta*theta); |
---|
| 197 | return Fkt_value; |
---|
| 198 | } |
---|
| 199 | |
---|
| 200 | G4double G4Generator2BN::Calculatedsdkdt(G4double kout, G4double theta, G4double Eel) const |
---|
| 201 | { |
---|
| 202 | |
---|
| 203 | G4double dsdkdt_value = 0.; |
---|
| 204 | G4double Z = 1; |
---|
| 205 | // classic radius (in cm) |
---|
| 206 | G4double r0 = 2.82E-13; |
---|
| 207 | //squared classic radius (in barn) |
---|
| 208 | G4double r02 = r0*r0*1.0E+24; |
---|
| 209 | |
---|
| 210 | // Photon energy cannot be greater than electron kinetic energy |
---|
| 211 | if(kout > (Eel-electron_mass_c2)){ |
---|
| 212 | dsdkdt_value = 0; |
---|
| 213 | return dsdkdt_value; |
---|
| 214 | } |
---|
| 215 | |
---|
| 216 | G4double E0 = Eel/electron_mass_c2; |
---|
| 217 | G4double k = kout/electron_mass_c2; |
---|
| 218 | G4double E = E0 - k; |
---|
| 219 | |
---|
| 220 | if(E <= 1*MeV ){ // Kinematic limit at 1 MeV !! |
---|
| 221 | dsdkdt_value = 0; |
---|
| 222 | return dsdkdt_value; |
---|
| 223 | } |
---|
| 224 | |
---|
| 225 | |
---|
| 226 | G4double p0 = std::sqrt(E0*E0-1); |
---|
| 227 | G4double p = std::sqrt(E*E-1); |
---|
| 228 | G4double L = std::log((E*E0-1+p*p0)/(E*E0-1-p*p0)); |
---|
| 229 | G4double delta0 = E0 - p0*std::cos(theta); |
---|
| 230 | G4double epsilon = std::log((E+p)/(E-p)); |
---|
| 231 | G4double Z2 = Z*Z; |
---|
| 232 | G4double sintheta2 = std::sin(theta)*std::sin(theta); |
---|
| 233 | G4double E02 = E0*E0; |
---|
| 234 | G4double E2 = E*E; |
---|
| 235 | G4double p02 = E0*E0-1; |
---|
| 236 | G4double k2 = k*k; |
---|
| 237 | G4double delta02 = delta0*delta0; |
---|
| 238 | G4double delta04 = delta02* delta02; |
---|
| 239 | G4double Q = std::sqrt(p02+k2-2*k*p0*std::cos(theta)); |
---|
| 240 | G4double Q2 = Q*Q; |
---|
| 241 | G4double epsilonQ = std::log((Q+p)/(Q-p)); |
---|
| 242 | |
---|
| 243 | |
---|
| 244 | dsdkdt_value = Z2 * (r02/(8*pi*137)) * (1/k) * (p/p0) * |
---|
| 245 | ( (8 * (sintheta2*(2*E02+1))/(p02*delta04)) - |
---|
| 246 | ((2*(5*E02+2*E*E0+3))/(p02 * delta02)) - |
---|
| 247 | ((2*(p02-k2))/((Q2*delta02))) + |
---|
| 248 | ((4*E)/(p02*delta0)) + |
---|
| 249 | (L/(p*p0))*( |
---|
| 250 | ((4*E0*sintheta2*(3*k-p02*E))/(p02*delta04)) + |
---|
| 251 | ((4*E02*(E02+E2))/(p02*delta02)) + |
---|
| 252 | ((2-2*(7*E02-3*E*E0+E2))/(p02*delta02)) + |
---|
| 253 | (2*k*(E02+E*E0-1))/((p02*delta0)) |
---|
| 254 | ) - |
---|
| 255 | ((4*epsilon)/(p*delta0)) + |
---|
| 256 | ((epsilonQ)/(p*Q))* |
---|
| 257 | (4/delta02-(6*k/delta0)-(2*k*(p02-k2))/(Q2*delta0)) |
---|
| 258 | ); |
---|
| 259 | |
---|
| 260 | |
---|
| 261 | |
---|
| 262 | dsdkdt_value = dsdkdt_value*std::sin(theta); |
---|
| 263 | return dsdkdt_value; |
---|
| 264 | } |
---|
| 265 | |
---|
| 266 | void G4Generator2BN::ConstructMajorantSurface() |
---|
| 267 | { |
---|
| 268 | |
---|
| 269 | G4double Eel; |
---|
| 270 | G4int vmax; |
---|
| 271 | G4double Ek; |
---|
| 272 | G4double k, theta, k0, theta0, thetamax; |
---|
| 273 | G4double ds, df, dsmax, dk, dt; |
---|
| 274 | G4double ratmin; |
---|
| 275 | G4double ratio = 0; |
---|
| 276 | G4double Vds, Vdf; |
---|
| 277 | G4double A, c; |
---|
| 278 | |
---|
| 279 | G4cout << "**** Constructing Majorant Surface for 2BN Distribution ****" << G4endl; |
---|
| 280 | |
---|
| 281 | if(kcut > kmin) kmin = kcut; |
---|
| 282 | |
---|
| 283 | G4int i = 0; |
---|
| 284 | // Loop on energy index |
---|
| 285 | for(G4int index = index_min; index < index_max; index++){ |
---|
| 286 | |
---|
| 287 | G4double fraction = index/100.; |
---|
| 288 | Ek = std::pow(10.,fraction); |
---|
| 289 | Eel = Ek + electron_mass_c2; |
---|
| 290 | |
---|
| 291 | // find x-section maximum at k=kmin |
---|
| 292 | dsmax = 0.; |
---|
| 293 | thetamax = 0.; |
---|
| 294 | |
---|
| 295 | for(theta = 0.; theta < pi; theta = theta + dtheta){ |
---|
| 296 | |
---|
| 297 | ds = Calculatedsdkdt(kmin, theta, Eel); |
---|
| 298 | |
---|
| 299 | if(ds > dsmax){ |
---|
| 300 | dsmax = ds; |
---|
| 301 | thetamax = theta; |
---|
| 302 | } |
---|
| 303 | } |
---|
| 304 | |
---|
| 305 | // Compute surface paremeters at kmin |
---|
| 306 | if(Ek < kmin || thetamax == 0){ |
---|
| 307 | c = 0; |
---|
| 308 | A = 0; |
---|
| 309 | }else{ |
---|
| 310 | c = 1/(thetamax*thetamax); |
---|
| 311 | A = 2*std::sqrt(c)*dsmax/(std::pow(kmin,-b)); |
---|
| 312 | } |
---|
| 313 | |
---|
| 314 | // look for correction factor to normalization at kmin |
---|
| 315 | ratmin = 1.; |
---|
| 316 | |
---|
| 317 | // Volume under surfaces |
---|
| 318 | Vds = 0.; |
---|
| 319 | Vdf = 0.; |
---|
| 320 | k0 = 0.; |
---|
| 321 | theta0 = 0.; |
---|
| 322 | |
---|
| 323 | vmax = G4int(100.*std::log10(Ek/kmin)); |
---|
| 324 | |
---|
| 325 | for(G4int v = 0; v < vmax; v++){ |
---|
| 326 | G4double fraction = (v/100.); |
---|
| 327 | k = std::pow(10.,fraction)*kmin; |
---|
| 328 | |
---|
| 329 | for(theta = 0.; theta < pi; theta = theta + dtheta){ |
---|
| 330 | dk = k - k0; |
---|
| 331 | dt = theta - theta0; |
---|
| 332 | ds = Calculatedsdkdt(k,theta, Eel); |
---|
| 333 | Vds = Vds + ds*dk*dt; |
---|
| 334 | df = CalculateFkt(k,theta, A, c); |
---|
| 335 | Vdf = Vdf + df*dk*dt; |
---|
| 336 | |
---|
| 337 | if(ds != 0.){ |
---|
| 338 | if(df != 0.) ratio = df/ds; |
---|
| 339 | } |
---|
| 340 | |
---|
| 341 | if(ratio < ratmin && ratio != 0.){ |
---|
| 342 | ratmin = ratio; |
---|
| 343 | } |
---|
| 344 | } |
---|
| 345 | } |
---|
| 346 | |
---|
| 347 | |
---|
| 348 | // sampling efficiency |
---|
| 349 | Atab[i] = A/ratmin * 1.04; |
---|
| 350 | ctab[i] = c; |
---|
| 351 | |
---|
| 352 | // G4cout << Ek << " " << i << " " << index << " " << Atab[i] << " " << ctab[i] << " " << G4endl; |
---|
| 353 | i++; |
---|
| 354 | } |
---|
| 355 | |
---|
| 356 | } |
---|
| 357 | |
---|
| 358 | G4double G4Generator2BN::Generate2BN(G4double Ek, G4double k) const |
---|
| 359 | { |
---|
| 360 | |
---|
| 361 | G4double Eel; |
---|
| 362 | G4double kmin2; |
---|
| 363 | G4double kmax, t; |
---|
| 364 | G4double cte2; |
---|
| 365 | G4double y, u; |
---|
| 366 | G4double fk, ft; |
---|
| 367 | G4double ds; |
---|
| 368 | G4double A2; |
---|
| 369 | G4double A, c; |
---|
| 370 | |
---|
| 371 | G4int trials = 0; |
---|
| 372 | G4int index; |
---|
| 373 | |
---|
| 374 | // find table index |
---|
| 375 | index = G4int(std::log10(Ek)*100) - index_min; |
---|
| 376 | Eel = Ek + electron_mass_c2; |
---|
| 377 | |
---|
| 378 | kmax = Ek; |
---|
| 379 | kmin2 = kcut; |
---|
| 380 | |
---|
| 381 | c = ctab[index]; |
---|
| 382 | A = Atab[index]; |
---|
| 383 | if(index < index_max){ |
---|
| 384 | A2 = Atab[index+1]; |
---|
| 385 | if(A2 > A) A = A2; |
---|
| 386 | } |
---|
| 387 | |
---|
| 388 | do{ |
---|
| 389 | // generate k accordimg to std::pow(k,-b) |
---|
| 390 | trials++; |
---|
| 391 | |
---|
| 392 | // normalization constant |
---|
| 393 | // cte1 = (1-b)/(std::pow(kmax,(1-b))-std::pow(kmin2,(1-b))); |
---|
| 394 | // y = G4UniformRand(); |
---|
| 395 | // k = std::pow(((1-b)/cte1*y+std::pow(kmin2,(1-b))),(1/(1-b))); |
---|
| 396 | |
---|
| 397 | // generate theta accordimg to theta/(1+c*std::pow(theta,2) |
---|
| 398 | // Normalization constant |
---|
| 399 | cte2 = 2*c/std::log(1+c*pi2); |
---|
| 400 | |
---|
| 401 | y = G4UniformRand(); |
---|
| 402 | t = std::sqrt((std::exp(2*c*y/cte2)-1)/c); |
---|
| 403 | u = G4UniformRand(); |
---|
| 404 | |
---|
| 405 | // point acceptance algorithm |
---|
| 406 | fk = std::pow(k,-b); |
---|
| 407 | ft = t/(1+c*t*t); |
---|
| 408 | ds = Calculatedsdkdt(k,t,Eel); |
---|
| 409 | |
---|
| 410 | // violation point |
---|
| 411 | if(ds > (A*fk*ft)) G4cout << "WARNING IN G4Generator2BN !!!" << Ek << " " << (ds-A*fk*ft)/ds << G4endl; |
---|
| 412 | |
---|
| 413 | }while(u*A*fk*ft > ds); |
---|
| 414 | |
---|
| 415 | return t; |
---|
| 416 | |
---|
| 417 | } |
---|
| 418 | |
---|
| 419 | void G4Generator2BN::PrintGeneratorInformation() const |
---|
| 420 | { |
---|
| 421 | G4cout << "\n" << G4endl; |
---|
| 422 | G4cout << "Bremsstrahlung Angular Generator is 2BN Generator from 2BN Koch & Motz distribution (Rev Mod Phys 31(4), 920 (1959))" << G4endl; |
---|
| 423 | G4cout << "\n" << G4endl; |
---|
| 424 | } |
---|
| 425 | |
---|