| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4RDGenerator2BN
33 //
34 // Author: Andreia Trindade (andreia@li
35 // Pedro Rodrigues (psilva@lip
36 // Luis Peralta (luis@lip.p
37 //
38 // Creation date: 21 June 2003
39 //
40 // Modifications:
41 // 21 Jun 2003
42 // 05 Nov 2003 MGP
43 // 14 Mar 2004
44 //
45 // Class Description:
46 //
47 // Concrete base class for Bremsstrahlung Angu
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.
63 0.023209, 0.0227892, 0.0225362, 0.
64 0.0204935, 0.0201227, 0.0197588, 0.
65 0.0179763, 0.0177866, 0.0174649, 0.
66 0.0160283, 0.0157384, 0.0155801, 0.
67 0.01432, 0.0140607, 0.0139267, 0.
68 0.0128205, 0.0125881, 0.012475, 0.
69 0.0115032, 0.0114067, 0.0111995, 0.
70 0.0104547, 0.0102646, 0.0101865, 0.
71 0.00943171,0.00936643,0.00930328,0.
72 0.00863611,0.00858428,0.00842757,0.
73 0.00784058,0.00779958,0.00776046,0.
74 0.00724201,0.00710969,0.00708004,0.
75 0.00662484,0.00650396,0.00648286,0.
76 0.00608412,0.00597331,0.00595991,0.
77 0.0056119, 0.00551005,0.00550377,0.
78 0.0051091, 0.00510777,0.00501582,0.
79 0.00475622,0.00476127,0.00467625,0.
80 0.00445522,0.00437664,0.00438768,0.
81 0.004203, 0.00413, 0.00405869,0.
82 0.00390795,0.00392975,0.00386339,0.
83 0.00374678,0.00368538,0.00371363,0.
84 0.00354139,0.00357481,0.00351921,0.
85 0.0034712, 0.00342026,0.00346165,0.
86 0.00335885,0.00340692,0.00336191,0.
87 0.00339463,0.00335506,0.00341401,0.
88 0.0033969, 0.00346557,0.00343302,0.
89 0.00357485,0.00354807,0.00352395,0.
90 0.00373078,0.00371234,0.00381532,0.
91 0.00398425,0.00411183,0.00410042,0.
92 0.00436082,0.00452075,0.00451934,0.
93 0.00490102,0.00510782,0.00511801,0.
94 0.0056677, 0.00594482,0.00596999,0.
95 0.00676236,0.00680914,0.00719407,0.
96 0.00835577,0.0084276, 0.00850759,0.
97 0.0101059, 0.0102249, 0.0109763, 0.
98 0.0125898, 0.0136311, 0.0138081, 0.
99 0.016168, 0.0176664, 0.0179339, 0.
100 0.0214374, 0.0237377, 0.0241275, 0.
101 0.0294621, 0.0300526, 0.0338619, 0.
102 };
103
104 G4double G4RDGenerator2BN::ctab[320] =
105 { 0.482253, 0.482253, 0.489021, 0.489021,
106 0.495933, 0.495933, 0.495933, 0.502993,
107 0.510204, 0.510204, 0.517572, 0.517572,
108 0.5251, 0.532793, 0.532793, 0.532793,
109 0.548697, 0.548697, 0.556917, 0.556917,
110 0.573921, 0.582717, 0.582717, 0.591716,
111 0.610352, 0.610352, 0.620001, 0.620001,
112 0.64, 0.650364, 0.650364, 0.660982,
113 0.683013, 0.694444, 0.694444, 0.706165,
114 0.730514, 0.743163, 0.743163, 0.756144,
115 0.783147, 0.797194, 0.797194, 0.811622,
116 0.84168, 0.857339, 0.857339, 0.873439,
117 0.907029, 0.924556, 0.924556, 0.942596,
118 0.980296, 1, 1, 1.0203,
119 1.06281, 1.06281, 1.08507, 1.08507,
120 1.13173, 1.1562, 1.1562, 1.18147,
121 1.23457, 1.23457, 1.26247, 1.26247,
122 1.32118, 1.35208, 1.35208, 1.38408,
123 1.45159, 1.45159, 1.45159, 1.48721,
124 1.5625, 1.5625, 1.60231, 1.60231,
125 1.68663, 1.68663, 1.7313, 1.7313,
126 1.82615, 1.87652, 1.87652, 1.92901,
127 1.98373, 2.04082, 2.04082, 2.1004,
128 2.22767, 2.22767, 2.22767, 2.29568,
129 2.44141, 2.44141, 2.51953, 2.51953,
130 2.68745, 2.68745, 2.77778, 2.77778,
131 2.97265, 2.97265, 3.07787, 3.07787,
132 3.30579, 3.30579, 3.42936, 3.42936,
133 3.69822, 3.84468, 3.84468, 3.84468,
134 4.16493, 4.34028, 4.34028, 4.34028,
135 4.7259, 4.93827, 4.93827, 4.93827,
136 5.40833, 5.40833, 5.66893, 5.66893,
137 6.25, 6.25, 6.57462, 6.57462,
138 7.3046, 7.3046, 7.71605, 7.71605,
139 8.65052, 8.65052, 8.65052, 9.18274,
140 9.76562, 10.4058, 10.4058, 10.4058,
141 11.8906, 11.8906, 12.7551, 12.7551,
142 13.7174, 14.7929, 14.7929, 14.7929,
143 17.3611, 17.3611, 17.3611, 18.9036,
144 20.6612, 20.6612, 22.6757, 22.6757,
145 25, 25, 27.7008, 27.7008,
146 30.8642, 30.8642, 34.6021, 34.6021,
147 39.0625, 39.0625, 39.0625, 44.4444,
148 51.0204, 51.0204, 51.0204, 51.0204,
149 59.1716, 59.1716, 69.4444, 69.4444,
150 82.6446, 82.6446, 82.6446, 82.6446,
151 };
152
153
154 G4RDGenerator2BN::G4RDGenerator2BN(const G4Str
155 {
156 b = 1.2;
157 index_min = -300;
158 index_max = 20;
159
160 // Set parameters minimum limits Ekelectron
161 kmin = 100*eV;
162 Ekmin = 250*eV;
163 kcut = 100*eV;
164
165 // Increment Theta value for surface interpo
166 dtheta = 0.1*rad;
167
168 // Construct Majorant Surface to 2BN cross-s
169 // (we dont need this. Pre-calculated values
170 // ConstructMajorantSurface();
171 }
172
173 //
174
175 G4RDGenerator2BN::~G4RDGenerator2BN()
176 {;}
177
178 //
179
180 G4double G4RDGenerator2BN::PolarAngle(const G4
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(G4doub
196 {
197 G4double Fkt_value = 0;
198 Fkt_value = A*std::pow(k,-b)*theta/(1+c*thet
199 return Fkt_value;
200 }
201
202 G4double G4RDGenerator2BN::Calculatedsdkdt(G4d
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 elec
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 ){
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
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:
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:
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
247 ( (8 * (sintheta2*(2*E02+1))/(p02*delta
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))
253 ((4*E02*(E02+E2))/(p02*delta0
254 ((2-2*(7*E02-3*E*E0+E2))/(p02
255 (2*k*(E02+E*E0-1))/((p02*delt
256 ) -
257 ((4*epsilon)/(p*delta0)) +
258 ((epsilonQ)/(p*Q))*
259 (4/delta02-(6*k/delta0)-(2*k*(p02-k2)
260 );
261
262
263
264 dsdkdt_value = dsdkdt_value*std::sin(thet
265 return dsdkdt_value;
266 }
267
268 void G4RDGenerator2BN::ConstructMajorantSurfac
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 Surfac
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_m
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 +
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 normalizati
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
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 <
355 i++;
356 }
357
358 }
359
360 G4double G4RDGenerator2BN::Generate2BN(G4doubl
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_mi
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::po
392 // y = G4UniformRand();
393 // k = std::pow(((1-b)/cte1*y+std::pow(kmin2,
394
395 // generate theta accordimg to theta/(1+c*st
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 G4R
410
411 }while(u*A*fk*ft > ds);
412
413 return t;
414
415 }
416
417 void G4RDGenerator2BN::PrintGeneratorInformati
418 {
419 G4cout << "\n" << G4endl;
420 G4cout << "Bremsstrahlung Angular Generator
421 G4cout << "\n" << G4endl;
422 }
423
424