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