| 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 // $Id: G4Generator2BN.cc 68038 2013-03-13 14:16:11Z gcosmo $
26 // 27 //
27 // ------------------------------------------- 28 // -------------------------------------------------------------------
28 // 29 //
29 // GEANT4 Class file 30 // GEANT4 Class file
30 // 31 //
31 // 32 //
32 // File name: G4Generator2BN 33 // File name: G4Generator2BN
33 // 34 //
34 // Author: Andreia Trindade (andreia@li 35 // Author: Andreia Trindade (andreia@lip.pt)
35 // Pedro Rodrigues (psilva@lip 36 // Pedro Rodrigues (psilva@lip.pt)
36 // Luis Peralta (luis@lip.p 37 // Luis Peralta (luis@lip.pt)
37 // 38 //
38 // Creation date: 21 June 2003 39 // Creation date: 21 June 2003
39 // 40 //
40 // Modifications: << 41 // Modifications:
41 // 21 Jun 2003 42 // 21 Jun 2003 First implementation acording with new design
42 // 05 Nov 2003 MGP 43 // 05 Nov 2003 MGP Fixed compilation warning
43 // 14 Mar 2004 44 // 14 Mar 2004 Code optimization
44 // 45 //
45 // Class Description: << 46 // Class Description:
46 // 47 //
47 // Concrete base class for Bremsstrahlung Angu << 48 // Concrete base class for Bremsstrahlung Angular Distribution Generation
48 // 2BN Distribution 49 // 2BN Distribution
49 // 50 //
50 // Class Description: End << 51 // Class Description: End
51 // 52 //
52 // ------------------------------------------- 53 // -------------------------------------------------------------------
53 // << 54 //
54 55
55 #include "G4Generator2BN.hh" 56 #include "G4Generator2BN.hh"
56 #include "Randomize.hh" 57 #include "Randomize.hh"
57 #include "G4PhysicalConstants.hh" 58 #include "G4PhysicalConstants.hh"
58 #include "G4SystemOfUnits.hh" 59 #include "G4SystemOfUnits.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 <<
156 155
157 G4Generator2BN::G4Generator2BN(const G4String& 156 G4Generator2BN::G4Generator2BN(const G4String&)
158 : G4VEmAngularDistribution("AngularGen2BN") 157 : G4VEmAngularDistribution("AngularGen2BN")
159 { 158 {
160 b = 1.2; 159 b = 1.2;
161 index_min = -300; 160 index_min = -300;
162 index_max = 319; 161 index_max = 319;
163 162
164 // Set parameters minimum limits Ekelectron 163 // Set parameters minimum limits Ekelectron = 250 eV and kphoton = 100 eV
165 kmin = 100*eV; 164 kmin = 100*eV;
166 Ekmin = 250*eV; 165 Ekmin = 250*eV;
167 kcut = 100*eV; 166 kcut = 100*eV;
168 167
169 // Increment Theta value for surface interpo 168 // Increment Theta value for surface interpolation
170 dtheta = 0.1*rad; 169 dtheta = 0.1*rad;
171 170
172 nwarn = 0; 171 nwarn = 0;
173 172
174 // Construct Majorant Surface to 2BN cross-s 173 // Construct Majorant Surface to 2BN cross-section
175 // (we dont need this. Pre-calculated values 174 // (we dont need this. Pre-calculated values are used instead due to performance issues
176 // ConstructMajorantSurface(); 175 // ConstructMajorantSurface();
177 } 176 }
178 177
179 //....oooOO0OOooo........oooOO0OOooo........oo << 178 //
>> 179
>> 180 G4Generator2BN::~G4Generator2BN()
>> 181 {}
180 182
181 G4ThreeVector& G4Generator2BN::SampleDirection 183 G4ThreeVector& G4Generator2BN::SampleDirection(const G4DynamicParticle* dp,
182 G4double out_energy, 184 G4double out_energy,
183 G4int, 185 G4int,
184 const G4Material*) 186 const G4Material*)
185 { 187 {
186 G4double Ek = dp->GetKineticEnergy(); 188 G4double Ek = dp->GetKineticEnergy();
187 G4double Eel = dp->GetTotalEnergy(); 189 G4double Eel = dp->GetTotalEnergy();
188 if(Eel > 2*MeV) { 190 if(Eel > 2*MeV) {
189 return fGenerator2BS.SampleDirection(dp, o 191 return fGenerator2BS.SampleDirection(dp, out_energy, 0);
190 } 192 }
191 193
192 G4double k = Eel - out_energy; 194 G4double k = Eel - out_energy;
193 195
194 G4double t; 196 G4double t;
195 G4double cte2; 197 G4double cte2;
196 G4double y, u; 198 G4double y, u;
197 G4double ds; 199 G4double ds;
198 G4double dmax; 200 G4double dmax;
199 201
>> 202 G4int trials = 0;
>> 203
200 // find table index 204 // find table index
201 G4int index = G4int(std::log10(Ek/MeV)*100) 205 G4int index = G4int(std::log10(Ek/MeV)*100) - index_min;
202 if(index > index_max) { index = index_max; } 206 if(index > index_max) { index = index_max; }
203 else if(index < 0) { index = 0; } 207 else if(index < 0) { index = 0; }
204 208
>> 209 //kmax = Ek;
>> 210 //kmin2 = kcut;
>> 211
205 G4double c = ctab[index]; 212 G4double c = ctab[index];
206 G4double A = Atab[index]; 213 G4double A = Atab[index];
207 if(index < index_max) { A = std::max(A, Atab 214 if(index < index_max) { A = std::max(A, Atab[index+1]); }
208 215
209 do{ 216 do{
210 // generate k accordimg to std::pow(k,-b) 217 // generate k accordimg to std::pow(k,-b)
>> 218 trials++;
211 219
212 // normalization constant << 220 // normalization constant
213 // cte1 = (1-b)/(std::pow(kmax,(1-b))-std 221 // cte1 = (1-b)/(std::pow(kmax,(1-b))-std::pow(kmin2,(1-b)));
214 // y = G4UniformRand(); 222 // y = G4UniformRand();
215 // k = std::pow(((1-b)/cte1*y+std::pow(km 223 // k = std::pow(((1-b)/cte1*y+std::pow(kmin2,(1-b))),(1/(1-b)));
216 224
217 // generate theta accordimg to theta/(1+c* 225 // generate theta accordimg to theta/(1+c*std::pow(theta,2)
218 // Normalization constant 226 // Normalization constant
219 cte2 = 2*c/std::log(1+c*pi2); 227 cte2 = 2*c/std::log(1+c*pi2);
220 228
221 y = G4UniformRand(); 229 y = G4UniformRand();
222 t = std::sqrt((G4Exp(2*c*y/cte2)-1)/c); << 230 t = std::sqrt((std::exp(2*c*y/cte2)-1)/c);
223 u = G4UniformRand(); 231 u = G4UniformRand();
224 232
225 // point acceptance algorithm 233 // point acceptance algorithm
226 //fk = std::pow(k,-b); 234 //fk = std::pow(k,-b);
227 //ft = t/(1+c*t*t); 235 //ft = t/(1+c*t*t);
228 dmax = A*std::pow(k,-b)*t/(1+c*t*t); 236 dmax = A*std::pow(k,-b)*t/(1+c*t*t);
229 ds = Calculatedsdkdt(k,t,Eel); 237 ds = Calculatedsdkdt(k,t,Eel);
230 238
231 // violation point 239 // violation point
232 if(ds > dmax && nwarn >= 20) { 240 if(ds > dmax && nwarn >= 20) {
233 ++nwarn; 241 ++nwarn;
234 G4cout << "### WARNING in G4Generator2BN << 242 G4cout << "### WARNING in G4Generator2BN: Ekin(MeV)= " << Ek/MeV
235 << " D(Ekin,k)/Dmax-1= " << (ds/dmax - 243 << " D(Ekin,k)/Dmax-1= " << (ds/dmax - 1)
236 << " results are not reliable!" << 244 << " results are not reliable!"
237 << G4endl; 245 << G4endl;
238 if(20 == nwarn) { << 246 if(20 == nwarn) {
239 G4cout << " WARNING in G4Generator2BN is c << 247 G4cout << " WARNING in G4Generator2BN is closed" << G4endl;
240 } 248 }
241 } 249 }
242 250
243 } while(u*dmax > ds || t > CLHEP::pi); 251 } while(u*dmax > ds || t > CLHEP::pi);
244 252
245 G4double sint = std::sin(t); 253 G4double sint = std::sin(t);
246 G4double phi = twopi*G4UniformRand(); << 254 G4double phi = twopi*G4UniformRand();
247 255
248 fLocalDirection.set(sint*std::cos(phi), sint 256 fLocalDirection.set(sint*std::cos(phi), sint*std::sin(phi),std::cos(t));
249 fLocalDirection.rotateUz(dp->GetMomentumDire 257 fLocalDirection.rotateUz(dp->GetMomentumDirection());
250 258
251 return fLocalDirection; 259 return fLocalDirection;
252 } 260 }
253 261
254 //....oooOO0OOooo........oooOO0OOooo........oo <<
255 <<
256 G4double G4Generator2BN::CalculateFkt(G4double 262 G4double G4Generator2BN::CalculateFkt(G4double k, G4double theta, G4double A, G4double c) const
257 { 263 {
258 G4double Fkt_value = 0; 264 G4double Fkt_value = 0;
259 Fkt_value = A*std::pow(k,-b)*theta/(1+c*thet 265 Fkt_value = A*std::pow(k,-b)*theta/(1+c*theta*theta);
260 return Fkt_value; 266 return Fkt_value;
261 } 267 }
262 268
263 //....oooOO0OOooo........oooOO0OOooo........oo <<
264 <<
265 G4double G4Generator2BN::Calculatedsdkdt(G4dou 269 G4double G4Generator2BN::Calculatedsdkdt(G4double kout, G4double theta, G4double Eel) const
266 { 270 {
>> 271
267 G4double dsdkdt_value = 0.; 272 G4double dsdkdt_value = 0.;
268 G4double Z = 1; 273 G4double Z = 1;
269 // classic radius (in cm) 274 // classic radius (in cm)
270 G4double r0 = 2.82E-13; 275 G4double r0 = 2.82E-13;
271 //squared classic radius (in barn) 276 //squared classic radius (in barn)
272 G4double r02 = r0*r0*1.0E+24; 277 G4double r02 = r0*r0*1.0E+24;
273 278
274 // Photon energy cannot be greater than elec 279 // Photon energy cannot be greater than electron kinetic energy
275 if(kout > (Eel-electron_mass_c2)){ 280 if(kout > (Eel-electron_mass_c2)){
276 dsdkdt_value = 0; 281 dsdkdt_value = 0;
277 return dsdkdt_value; 282 return dsdkdt_value;
278 } 283 }
279 284
280 G4double E0 = Eel/electron_mass_c2; << 285 G4double E0 = Eel/electron_mass_c2;
281 G4double k = kout/electron_mass_c2; << 286 G4double k = kout/electron_mass_c2;
282 G4double E = E0 - k; << 287 G4double E = E0 - k;
283 << 288
284 if(E <= 1*MeV ){ << 289 if(E <= 1*MeV ){ // Kinematic limit at 1 MeV !!
285 dsdkdt_value = 0; << 290 dsdkdt_value = 0;
286 return dsdkdt_value; << 291 return dsdkdt_value;
287 } << 292 }
288 << 293
289 G4double p0 = std::sqrt(E0*E0-1); << 294
290 G4double p = std::sqrt(E*E-1); << 295 G4double p0 = std::sqrt(E0*E0-1);
291 G4double LL = std::log((E*E0-1+p*p0)/(E*E0-1 << 296 G4double p = std::sqrt(E*E-1);
292 G4double delta0 = E0 - p0*std::cos(theta); << 297 G4double L = std::log((E*E0-1+p*p0)/(E*E0-1-p*p0));
293 G4double epsilon = std::log((E+p)/(E-p)); << 298 G4double delta0 = E0 - p0*std::cos(theta);
294 G4double Z2 = Z*Z; << 299 G4double epsilon = std::log((E+p)/(E-p));
295 G4double sintheta2 = std::sin(theta)*std::si << 300 G4double Z2 = Z*Z;
296 G4double E02 = E0*E0; << 301 G4double sintheta2 = std::sin(theta)*std::sin(theta);
297 G4double E2 = E*E; << 302 G4double E02 = E0*E0;
298 G4double p02 = E0*E0-1; << 303 G4double E2 = E*E;
299 G4double k2 = k*k; << 304 G4double p02 = E0*E0-1;
300 G4double delta02 = delta0*delta0; << 305 G4double k2 = k*k;
301 G4double delta04 = delta02* delta02; << 306 G4double delta02 = delta0*delta0;
302 G4double Q = std::sqrt(p02+k2-2*k*p0*std::co << 307 G4double delta04 = delta02* delta02;
303 G4double Q2 = Q*Q; << 308 G4double Q = std::sqrt(p02+k2-2*k*p0*std::cos(theta));
304 G4double epsilonQ = std::log((Q+p)/(Q-p)); << 309 G4double Q2 = Q*Q;
305 << 310 G4double epsilonQ = std::log((Q+p)/(Q-p));
306 << 311
307 dsdkdt_value = Z2 * (r02/(8*pi*137)) * (1/k) << 312
308 ( (8 * (sintheta2*(2*E02+1))/(p02*delta04) << 313 dsdkdt_value = Z2 * (r02/(8*pi*137)) * (1/k) * (p/p0) *
309 ((2*(5*E02+2*E*E0+3))/(p02 * delta02)) - << 314 ( (8 * (sintheta2*(2*E02+1))/(p02*delta04)) -
310 ((2*(p02-k2))/((Q2*delta02))) + << 315 ((2*(5*E02+2*E*E0+3))/(p02 * delta02)) -
311 ((4*E)/(p02*delta0)) + << 316 ((2*(p02-k2))/((Q2*delta02))) +
312 (LL/(p*p0))*( << 317 ((4*E)/(p02*delta0)) +
313 ((4*E0*sintheta2*(3*k-p02*E))/(p02*delt << 318 (L/(p*p0))*(
314 ((4*E02*(E02+E2))/(p02*delta02)) + << 319 ((4*E0*sintheta2*(3*k-p02*E))/(p02*delta04)) +
315 ((2-2*(7*E02-3*E*E0+E2))/(p02*delta02)) << 320 ((4*E02*(E02+E2))/(p02*delta02)) +
316 (2*k*(E02+E*E0-1))/((p02*delta0)) << 321 ((2-2*(7*E02-3*E*E0+E2))/(p02*delta02)) +
317 ) - << 322 (2*k*(E02+E*E0-1))/((p02*delta0))
318 ((4*epsilon)/(p*delta0)) + << 323 ) -
319 ((epsilonQ)/(p*Q))* << 324 ((4*epsilon)/(p*delta0)) +
320 (4/delta02-(6*k/delta0)-(2*k*(p02-k2))/( << 325 ((epsilonQ)/(p*Q))*
321 ); << 326 (4/delta02-(6*k/delta0)-(2*k*(p02-k2))/(Q2*delta0))
322 << 327 );
323 dsdkdt_value = dsdkdt_value*std::sin(theta); << 328
324 return dsdkdt_value; <<
325 } <<
326 329
327 //....oooOO0OOooo........oooOO0OOooo........oo << 330
>> 331 dsdkdt_value = dsdkdt_value*std::sin(theta);
>> 332 return dsdkdt_value;
>> 333 }
328 334
329 void G4Generator2BN::ConstructMajorantSurface( 335 void G4Generator2BN::ConstructMajorantSurface()
330 { 336 {
>> 337
331 G4double Eel; 338 G4double Eel;
332 G4int vmax; 339 G4int vmax;
333 G4double Ek; 340 G4double Ek;
334 G4double k, theta, k0, theta0, thetamax; 341 G4double k, theta, k0, theta0, thetamax;
335 G4double ds, df, dsmax, dk, dt; 342 G4double ds, df, dsmax, dk, dt;
336 G4double ratmin; 343 G4double ratmin;
337 G4double ratio = 0; 344 G4double ratio = 0;
338 G4double Vds, Vdf; 345 G4double Vds, Vdf;
339 G4double A, c; 346 G4double A, c;
340 347
341 G4cout << "**** Constructing Majorant Surfac 348 G4cout << "**** Constructing Majorant Surface for 2BN Distribution ****" << G4endl;
342 349
343 if(kcut > kmin) kmin = kcut; 350 if(kcut > kmin) kmin = kcut;
344 351
345 G4int i = 0; 352 G4int i = 0;
346 // Loop on energy index 353 // Loop on energy index
347 for(G4int index = index_min; index < index_m 354 for(G4int index = index_min; index < index_max; index++){
348 355
349 G4double fraction = index/100.; 356 G4double fraction = index/100.;
350 Ek = std::pow(10.,fraction); 357 Ek = std::pow(10.,fraction);
351 Eel = Ek + electron_mass_c2; 358 Eel = Ek + electron_mass_c2;
352 359
353 // find x-section maximum at k=kmin 360 // find x-section maximum at k=kmin
354 dsmax = 0.; 361 dsmax = 0.;
355 thetamax = 0.; 362 thetamax = 0.;
356 363
357 for(theta = 0.; theta < pi; theta = theta + 364 for(theta = 0.; theta < pi; theta = theta + dtheta){
358 365
359 ds = Calculatedsdkdt(kmin, theta, Eel); 366 ds = Calculatedsdkdt(kmin, theta, Eel);
360 367
361 if(ds > dsmax){ 368 if(ds > dsmax){
362 dsmax = ds; 369 dsmax = ds;
363 thetamax = theta; 370 thetamax = theta;
364 } 371 }
365 } 372 }
366 373
367 // Compute surface paremeters at kmin 374 // Compute surface paremeters at kmin
368 if(Ek < kmin || thetamax == 0){ 375 if(Ek < kmin || thetamax == 0){
369 c = 0; 376 c = 0;
370 A = 0; 377 A = 0;
371 }else{ 378 }else{
372 c = 1/(thetamax*thetamax); 379 c = 1/(thetamax*thetamax);
373 A = 2*std::sqrt(c)*dsmax/(std::pow(kmin,-b 380 A = 2*std::sqrt(c)*dsmax/(std::pow(kmin,-b));
374 } 381 }
375 382
376 // look for correction factor to normalizati << 383 // look for correction factor to normalization at kmin
377 ratmin = 1.; 384 ratmin = 1.;
378 385
379 // Volume under surfaces 386 // Volume under surfaces
380 Vds = 0.; 387 Vds = 0.;
381 Vdf = 0.; 388 Vdf = 0.;
382 k0 = 0.; 389 k0 = 0.;
383 theta0 = 0.; 390 theta0 = 0.;
384 391
385 vmax = G4int(100.*std::log10(Ek/kmin)); 392 vmax = G4int(100.*std::log10(Ek/kmin));
386 393
387 for(G4int v = 0; v < vmax; v++){ 394 for(G4int v = 0; v < vmax; v++){
388 G4double fractionLocal = (v/100.); 395 G4double fractionLocal = (v/100.);
389 k = std::pow(10.,fractionLocal)*kmin; 396 k = std::pow(10.,fractionLocal)*kmin;
390 397
391 for(theta = 0.; theta < pi; theta = theta 398 for(theta = 0.; theta < pi; theta = theta + dtheta){
392 dk = k - k0; 399 dk = k - k0;
393 dt = theta - theta0; 400 dt = theta - theta0;
394 ds = Calculatedsdkdt(k,theta, Eel); 401 ds = Calculatedsdkdt(k,theta, Eel);
395 Vds = Vds + ds*dk*dt; 402 Vds = Vds + ds*dk*dt;
396 df = CalculateFkt(k,theta, A, c); 403 df = CalculateFkt(k,theta, A, c);
397 Vdf = Vdf + df*dk*dt; 404 Vdf = Vdf + df*dk*dt;
398 405
399 if(ds != 0.){ 406 if(ds != 0.){
400 if(df != 0.) ratio = df/ds; 407 if(df != 0.) ratio = df/ds;
401 } 408 }
402 409
403 if(ratio < ratmin && ratio != 0.){ 410 if(ratio < ratmin && ratio != 0.){
404 ratmin = ratio; 411 ratmin = ratio;
405 } 412 }
406 } 413 }
407 } 414 }
408 415
>> 416
409 // sampling efficiency 417 // sampling efficiency
410 Atab[i] = A/ratmin * 1.04; 418 Atab[i] = A/ratmin * 1.04;
411 ctab[i] = c; 419 ctab[i] = c;
412 420
413 // G4cout << Ek << " " << i << " " << index < 421 // G4cout << Ek << " " << i << " " << index << " " << Atab[i] << " " << ctab[i] << " " << G4endl;
414 i++; 422 i++;
415 } 423 }
416 } <<
417 424
418 //....oooOO0OOooo........oooOO0OOooo........oo << 425 }
419 426
420 void G4Generator2BN::PrintGeneratorInformation 427 void G4Generator2BN::PrintGeneratorInformation() const
421 { 428 {
422 G4cout << "\n" << G4endl; 429 G4cout << "\n" << G4endl;
423 G4cout << "Bremsstrahlung Angular Generator 430 G4cout << "Bremsstrahlung Angular Generator is 2BN Generator from 2BN Koch & Motz distribution (Rev Mod Phys 31(4), 920 (1959))" << G4endl;
424 G4cout << "\n" << G4endl; 431 G4cout << "\n" << G4endl;
425 } << 432 }
>> 433
426 434