| [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 |
|
|---|
