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