ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Generator2BN.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Generator2BN.cc
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: G4Generator2BN
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
48 // 2BN Distribution
49 //
50 // Class Description: End
51 //
52 // -------------------------------------------------------------------
53 //
54 
55 #include "G4Generator2BN.hh"
56 #include "Randomize.hh"
57 #include "G4PhysicalConstants.hh"
58 #include "G4SystemOfUnits.hh"
59 #include "G4Exp.hh"
60 
61 //
62 
63 G4double G4Generator2BN::Atab[320] =
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.0221284, 0.0217282,0.0214894, 0.0211005, 0.0207189,
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.0172828, 0.0169702,0.0166634, 0.0164915, 0.0161933,
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.0136744, 0.0135459,0.0133005, 0.0131773, 0.0129385,
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.0111072, 0.0110175,0.0108173, 0.0107316, 0.0105365,
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.0091337, 0.00907372,0.00890831,0.00885141,0.00869003,
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.00761866,0.00758201,0.00744346,0.00740928,0.00727384,
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.00636458,0.00634545,0.00622977,0.00621258,0.00609936,
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.00540399,0.00539938,0.00530162,0.00529872,0.00520292,
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.00468287,0.00459947,0.00451785,0.00452581,0.00444573,
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.00407563,0.00400561,0.00402414,0.00395536,0.00397553,
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.00365335,0.00359463,0.0036242, 0.00356653,0.003598,
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.00341172,0.00345548,0.00340657,0.00335944,0.00340491,
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.00337558,0.00343797,0.00340067,0.00336584,0.00343059,
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.00360885,0.00358571,0.00367661,0.00365446,0.00375194,
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.00409197,0.00422843,0.00422123,0.00436974,0.0043637,
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.00513271,0.0053589, 0.00537524,0.00562643,0.00564452,
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.0072439, 0.00730063,0.0077349, 0.00779494,0.00786293,
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.0111003, 0.0112367, 0.0113862, 0.0122637, 0.0124196,
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.0182246, 0.0185396, 0.020381, 0.0207026, 0.0210558,
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.0344537, 0.0351108, 0.0358099, 0.036554, 0.0416399
104  };
105 
106 G4double G4Generator2BN::ctab[320] =
107  { 0.482253, 0.482253, 0.489021, 0.489021, 0.489021, 0.489021, 0.495933,
108  0.495933, 0.495933, 0.495933, 0.502993, 0.502993, 0.502993, 0.510204,
109  0.510204, 0.510204, 0.517572, 0.517572, 0.517572, 0.5251, 0.5251,
110  0.5251, 0.532793, 0.532793, 0.532793, 0.540657, 0.540657, 0.548697,
111  0.548697, 0.548697, 0.556917, 0.556917, 0.565323, 0.565323, 0.573921,
112  0.573921, 0.582717, 0.582717, 0.591716, 0.591716, 0.600925, 0.600925,
113  0.610352, 0.610352, 0.620001, 0.620001, 0.629882, 0.629882, 0.64,
114  0.64, 0.650364, 0.650364, 0.660982, 0.660982, 0.671862, 0.683013,
115  0.683013, 0.694444, 0.694444, 0.706165, 0.718184, 0.718184, 0.730514,
116  0.730514, 0.743163, 0.743163, 0.756144, 0.769468, 0.769468, 0.783147,
117  0.783147, 0.797194, 0.797194, 0.811622, 0.826446, 0.826446, 0.84168,
118  0.84168, 0.857339, 0.857339, 0.873439, 0.889996, 0.889996, 0.907029,
119  0.907029, 0.924556, 0.924556, 0.942596, 0.942596, 0.961169, 0.980296,
120  0.980296, 1, 1, 1.0203, 1.0203, 1.04123, 1.04123,
121  1.06281, 1.06281, 1.08507, 1.08507, 1.10803, 1.10803, 1.13173,
122  1.13173, 1.1562, 1.1562, 1.18147, 1.18147, 1.20758, 1.20758,
123  1.23457, 1.23457, 1.26247, 1.26247, 1.29132, 1.29132, 1.32118,
124  1.32118, 1.35208, 1.35208, 1.38408, 1.38408, 1.41723, 1.41723,
125  1.45159, 1.45159, 1.45159, 1.48721, 1.48721, 1.52416, 1.52416,
126  1.5625, 1.5625, 1.60231, 1.60231, 1.64366, 1.64366, 1.68663,
127  1.68663, 1.68663, 1.7313, 1.7313, 1.77778, 1.77778, 1.82615,
128  1.82615, 1.87652, 1.87652, 1.92901, 1.92901, 1.98373, 1.98373,
129  1.98373, 2.04082, 2.04082, 2.1004, 2.1004, 2.16263, 2.16263,
130  2.22767, 2.22767, 2.22767, 2.29568, 2.29568, 2.36686, 2.36686,
131  2.44141, 2.44141, 2.51953, 2.51953, 2.51953, 2.60146, 2.60146,
132  2.68745, 2.68745, 2.77778, 2.77778, 2.87274, 2.87274, 2.87274,
133  2.97265, 2.97265, 3.07787, 3.07787, 3.18878, 3.18878, 3.30579,
134  3.30579, 3.30579, 3.42936, 3.42936, 3.55999, 3.55999, 3.69822,
135  3.69822, 3.84468, 3.84468, 3.84468, 4, 4, 4.16493,
136  4.16493, 4.34028, 4.34028, 4.34028, 4.52694, 4.52694, 4.7259,
137  4.7259, 4.93827, 4.93827, 4.93827, 5.16529, 5.16529, 5.40833,
138  5.40833, 5.40833, 5.66893, 5.66893, 5.94884, 5.94884, 6.25,
139  6.25, 6.25, 6.57462, 6.57462, 6.92521, 6.92521, 6.92521,
140  7.3046, 7.3046, 7.71605, 7.71605, 7.71605, 8.16327, 8.16327,
141  8.65052, 8.65052, 8.65052, 9.18274, 9.18274, 9.18274, 9.76562,
142  9.76562, 10.4058, 10.4058, 10.4058, 11.1111, 11.1111, 11.1111,
143  11.8906, 11.8906, 12.7551, 12.7551, 12.7551, 13.7174, 13.7174,
144  13.7174, 14.7929, 14.7929, 14.7929, 16, 16, 16,
145  17.3611, 17.3611, 17.3611, 18.9036, 18.9036, 18.9036, 20.6612,
146  20.6612, 20.6612, 22.6757, 22.6757, 22.6757, 22.6757, 25,
147  25, 25, 27.7008, 27.7008, 27.7008, 27.7008, 30.8642,
148  30.8642, 30.8642, 34.6021, 34.6021, 34.6021, 34.6021, 39.0625,
149  39.0625, 39.0625, 39.0625, 44.4444, 44.4444, 44.4444, 44.4444,
150  51.0204, 51.0204, 51.0204, 51.0204, 59.1716, 59.1716, 59.1716,
151  59.1716, 59.1716, 69.4444, 69.4444, 69.4444, 69.4444, 82.6446,
152  82.6446, 82.6446, 82.6446, 82.6446, 100
153  };
154 
155 
156 G4Generator2BN::G4Generator2BN(const G4String&)
157  : G4VEmAngularDistribution("AngularGen2BN")
158 {
159  b = 1.2;
160  index_min = -300;
161  index_max = 319;
162 
163  // Set parameters minimum limits Ekelectron = 250 eV and kphoton = 100 eV
164  kmin = 100*eV;
165  Ekmin = 250*eV;
166  kcut = 100*eV;
167 
168  // Increment Theta value for surface interpolation
169  dtheta = 0.1*rad;
170 
171  nwarn = 0;
172 
173  // Construct Majorant Surface to 2BN cross-section
174  // (we dont need this. Pre-calculated values are used instead due to performance issues
175  // ConstructMajorantSurface();
176 }
177 
178 //
179 
180 G4Generator2BN::~G4Generator2BN()
181 {}
182 
183 G4ThreeVector& G4Generator2BN::SampleDirection(const G4DynamicParticle* dp,
184  G4double out_energy,
185  G4int,
186  const G4Material*)
187 {
188  G4double Ek = dp->GetKineticEnergy();
189  G4double Eel = dp->GetTotalEnergy();
190  if(Eel > 2*MeV) {
191  return fGenerator2BS.SampleDirection(dp, out_energy, 0);
192  }
193 
194  G4double k = Eel - out_energy;
195 
196  G4double t;
197  G4double cte2;
198  G4double y, u;
199  G4double ds;
200  G4double dmax;
201 
202  G4int trials = 0;
203 
204  // find table index
205  G4int index = G4int(std::log10(Ek/MeV)*100) - index_min;
206  if(index > index_max) { index = index_max; }
207  else if(index < 0) { index = 0; }
208 
209  //kmax = Ek;
210  //kmin2 = kcut;
211 
212  G4double c = ctab[index];
213  G4double A = Atab[index];
214  if(index < index_max) { A = std::max(A, Atab[index+1]); }
215 
216  do{
217  // generate k accordimg to std::pow(k,-b)
218  trials++;
219 
220  // normalization constant
221  // cte1 = (1-b)/(std::pow(kmax,(1-b))-std::pow(kmin2,(1-b)));
222  // y = G4UniformRand();
223  // k = std::pow(((1-b)/cte1*y+std::pow(kmin2,(1-b))),(1/(1-b)));
224 
225  // generate theta accordimg to theta/(1+c*std::pow(theta,2)
226  // Normalization constant
227  cte2 = 2*c/std::log(1+c*pi2);
228 
229  y = G4UniformRand();
230  t = std::sqrt((G4Exp(2*c*y/cte2)-1)/c);
231  u = G4UniformRand();
232 
233  // point acceptance algorithm
234  //fk = std::pow(k,-b);
235  //ft = t/(1+c*t*t);
236  dmax = A*std::pow(k,-b)*t/(1+c*t*t);
237  ds = Calculatedsdkdt(k,t,Eel);
238 
239  // violation point
240  if(ds > dmax && nwarn >= 20) {
241  ++nwarn;
242  G4cout << "### WARNING in G4Generator2BN: Ekin(MeV)= " << Ek/MeV
243  << " D(Ekin,k)/Dmax-1= " << (ds/dmax - 1)
244  << " results are not reliable!"
245  << G4endl;
246  if(20 == nwarn) {
247  G4cout << " WARNING in G4Generator2BN is closed" << G4endl;
248  }
249  }
250 
251  } while(u*dmax > ds || t > CLHEP::pi);
252 
253  G4double sint = std::sin(t);
254  G4double phi = twopi*G4UniformRand();
255 
256  fLocalDirection.set(sint*std::cos(phi), sint*std::sin(phi),std::cos(t));
257  fLocalDirection.rotateUz(dp->GetMomentumDirection());
258 
259  return fLocalDirection;
260 }
261 
262 G4double G4Generator2BN::CalculateFkt(G4double k, G4double theta, G4double A, G4double c) const
263 {
264  G4double Fkt_value = 0;
265  Fkt_value = A*std::pow(k,-b)*theta/(1+c*theta*theta);
266  return Fkt_value;
267 }
268 
269 G4double G4Generator2BN::Calculatedsdkdt(G4double kout, G4double theta, G4double Eel) const
270 {
271 
272  G4double dsdkdt_value = 0.;
273  G4double Z = 1;
274  // classic radius (in cm)
275  G4double r0 = 2.82E-13;
276  //squared classic radius (in barn)
277  G4double r02 = r0*r0*1.0E+24;
278 
279  // Photon energy cannot be greater than electron kinetic energy
280  if(kout > (Eel-electron_mass_c2)){
281  dsdkdt_value = 0;
282  return dsdkdt_value;
283  }
284 
285  G4double E0 = Eel/electron_mass_c2;
286  G4double k = kout/electron_mass_c2;
287  G4double E = E0 - k;
288 
289  if(E <= 1*MeV ){ // Kinematic limit at 1 MeV !!
290  dsdkdt_value = 0;
291  return dsdkdt_value;
292  }
293 
294 
295  G4double p0 = std::sqrt(E0*E0-1);
296  G4double p = std::sqrt(E*E-1);
297  G4double LL = std::log((E*E0-1+p*p0)/(E*E0-1-p*p0));
298  G4double delta0 = E0 - p0*std::cos(theta);
299  G4double epsilon = std::log((E+p)/(E-p));
300  G4double Z2 = Z*Z;
301  G4double sintheta2 = std::sin(theta)*std::sin(theta);
302  G4double E02 = E0*E0;
303  G4double E2 = E*E;
304  G4double p02 = E0*E0-1;
305  G4double k2 = k*k;
306  G4double delta02 = delta0*delta0;
307  G4double delta04 = delta02* delta02;
308  G4double Q = std::sqrt(p02+k2-2*k*p0*std::cos(theta));
309  G4double Q2 = Q*Q;
310  G4double epsilonQ = std::log((Q+p)/(Q-p));
311 
312 
313  dsdkdt_value = Z2 * (r02/(8*pi*137)) * (1/k) * (p/p0) *
314  ( (8 * (sintheta2*(2*E02+1))/(p02*delta04)) -
315  ((2*(5*E02+2*E*E0+3))/(p02 * delta02)) -
316  ((2*(p02-k2))/((Q2*delta02))) +
317  ((4*E)/(p02*delta0)) +
318  (LL/(p*p0))*(
319  ((4*E0*sintheta2*(3*k-p02*E))/(p02*delta04)) +
320  ((4*E02*(E02+E2))/(p02*delta02)) +
321  ((2-2*(7*E02-3*E*E0+E2))/(p02*delta02)) +
322  (2*k*(E02+E*E0-1))/((p02*delta0))
323  ) -
324  ((4*epsilon)/(p*delta0)) +
325  ((epsilonQ)/(p*Q))*
326  (4/delta02-(6*k/delta0)-(2*k*(p02-k2))/(Q2*delta0))
327  );
328 
329 
330 
331  dsdkdt_value = dsdkdt_value*std::sin(theta);
332  return dsdkdt_value;
333 }
334 
335 void G4Generator2BN::ConstructMajorantSurface()
336 {
337 
338  G4double Eel;
339  G4int vmax;
340  G4double Ek;
341  G4double k, theta, k0, theta0, thetamax;
342  G4double ds, df, dsmax, dk, dt;
343  G4double ratmin;
344  G4double ratio = 0;
345  G4double Vds, Vdf;
346  G4double A, c;
347 
348  G4cout << "**** Constructing Majorant Surface for 2BN Distribution ****" << G4endl;
349 
350  if(kcut > kmin) kmin = kcut;
351 
352  G4int i = 0;
353  // Loop on energy index
354  for(G4int index = index_min; index < index_max; index++){
355 
356  G4double fraction = index/100.;
357  Ek = std::pow(10.,fraction);
358  Eel = Ek + electron_mass_c2;
359 
360  // find x-section maximum at k=kmin
361  dsmax = 0.;
362  thetamax = 0.;
363 
364  for(theta = 0.; theta < pi; theta = theta + dtheta){
365 
366  ds = Calculatedsdkdt(kmin, theta, Eel);
367 
368  if(ds > dsmax){
369  dsmax = ds;
370  thetamax = theta;
371  }
372  }
373 
374  // Compute surface paremeters at kmin
375  if(Ek < kmin || thetamax == 0){
376  c = 0;
377  A = 0;
378  }else{
379  c = 1/(thetamax*thetamax);
380  A = 2*std::sqrt(c)*dsmax/(std::pow(kmin,-b));
381  }
382 
383  // look for correction factor to normalization at kmin
384  ratmin = 1.;
385 
386  // Volume under surfaces
387  Vds = 0.;
388  Vdf = 0.;
389  k0 = 0.;
390  theta0 = 0.;
391 
392  vmax = G4int(100.*std::log10(Ek/kmin));
393 
394  for(G4int v = 0; v < vmax; v++){
395  G4double fractionLocal = (v/100.);
396  k = std::pow(10.,fractionLocal)*kmin;
397 
398  for(theta = 0.; theta < pi; theta = theta + dtheta){
399  dk = k - k0;
400  dt = theta - theta0;
401  ds = Calculatedsdkdt(k,theta, Eel);
402  Vds = Vds + ds*dk*dt;
403  df = CalculateFkt(k,theta, A, c);
404  Vdf = Vdf + df*dk*dt;
405 
406  if(ds != 0.){
407  if(df != 0.) ratio = df/ds;
408  }
409 
410  if(ratio < ratmin && ratio != 0.){
411  ratmin = ratio;
412  }
413  }
414  }
415 
416 
417  // sampling efficiency
418  Atab[i] = A/ratmin * 1.04;
419  ctab[i] = c;
420 
421 // G4cout << Ek << " " << i << " " << index << " " << Atab[i] << " " << ctab[i] << " " << G4endl;
422  i++;
423  }
424 
425 }
426 
427 void G4Generator2BN::PrintGeneratorInformation() const
428 {
429  G4cout << "\n" << G4endl;
430  G4cout << "Bremsstrahlung Angular Generator is 2BN Generator from 2BN Koch & Motz distribution (Rev Mod Phys 31(4), 920 (1959))" << G4endl;
431  G4cout << "\n" << G4endl;
432 }