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