ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RDGenerator2BN.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RDGenerator2BN.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: 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