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