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