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