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