source: trunk/source/processes/electromagnetic/lowenergy/src/G4Generator2BN.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

geant4 tag 9.4
File size: 16.2 KB
Line 
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,v 1.9 2010/10/14 14:01:02 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name:     G4Generator2BN
35//
36// Author:        Andreia Trindade (andreia@lip.pt)
37//                Pedro Rodrigues  (psilva@lip.pt)
38//                Luis Peralta     (luis@lip.pt)
39//
40// Creation date: 21 June 2003
41//
42// Modifications:
43// 21 Jun 2003                                 First implementation acording with new design
44// 05 Nov 2003  MGP                            Fixed compilation warning
45// 14 Mar 2004                                 Code optimization
46//
47// Class Description:
48//
49// Concrete base class for Bremsstrahlung Angular Distribution Generation
50// 2BN Distribution
51//
52// Class Description: End
53//
54// -------------------------------------------------------------------
55//   
56
57#include "G4Generator2BN.hh"
58#include "Randomize.hh"
59//   
60
61G4double G4Generator2BN::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
104G4double G4Generator2BN::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
154G4Generator2BN::G4Generator2BN(const G4String&)
155 : G4VBremAngularDistribution("AngularGen2BN")
156{
157  b = 1.2;
158  index_min = -300;
159  index_max = 20;
160
161  // Set parameters minimum limits Ekelectron = 250 eV and kphoton = 100 eV
162  kmin = 100*eV;
163  Ekmin = 250*eV;
164  kcut = 100*eV;
165
166  // Increment Theta value for surface interpolation
167  dtheta = 0.1*rad;
168
169  // Construct Majorant Surface to 2BN cross-section
170  // (we dont need this. Pre-calculated values are used instead due to performance issues
171  // ConstructMajorantSurface();
172}
173
174//   
175
176G4Generator2BN::~G4Generator2BN() 
177{}
178
179//
180
181G4double G4Generator2BN::PolarAngle(const G4double initial_energy,
182                                    const G4double final_energy,
183                                    const G4int) // Z
184{
185
186  G4double theta = 0;
187
188  G4double k = initial_energy - final_energy;
189
190  theta = Generate2BN(initial_energy, k);
191
192  return theta;
193}
194//
195
196G4double G4Generator2BN::CalculateFkt(G4double k, G4double theta, G4double A, G4double c) const
197{
198  G4double Fkt_value = 0;
199  Fkt_value = A*std::pow(k,-b)*theta/(1+c*theta*theta);
200  return Fkt_value;
201}
202
203G4double G4Generator2BN::Calculatedsdkdt(G4double kout, G4double theta, G4double Eel) const
204{
205
206  G4double dsdkdt_value = 0.;
207  G4double Z = 1;
208  // classic radius (in cm)
209  G4double r0 = 2.82E-13;
210  //squared classic radius (in barn)
211  G4double r02 = r0*r0*1.0E+24;
212
213  // Photon energy cannot be greater than electron kinetic energy
214  if(kout > (Eel-electron_mass_c2)){
215    dsdkdt_value = 0;
216    return dsdkdt_value;
217  }
218
219     G4double E0 = Eel/electron_mass_c2;
220     G4double k =  kout/electron_mass_c2;
221     G4double E =  E0 - k;
222
223     if(E <= 1*MeV ){                                 // Kinematic limit at 1 MeV !!
224       dsdkdt_value = 0;
225       return dsdkdt_value;
226     }
227
228
229     G4double p0 = std::sqrt(E0*E0-1);
230     G4double p = std::sqrt(E*E-1);
231     G4double L = std::log((E*E0-1+p*p0)/(E*E0-1-p*p0));
232     G4double delta0 = E0 - p0*std::cos(theta);
233     G4double epsilon = std::log((E+p)/(E-p));
234     G4double Z2 = Z*Z;
235     G4double sintheta2 = std::sin(theta)*std::sin(theta);
236     G4double E02 = E0*E0;
237     G4double E2 = E*E;
238     G4double p02 = E0*E0-1;
239     G4double k2 = k*k;
240     G4double delta02 = delta0*delta0;
241     G4double delta04 =  delta02* delta02;
242     G4double Q = std::sqrt(p02+k2-2*k*p0*std::cos(theta));
243     G4double Q2 = Q*Q;
244     G4double epsilonQ = std::log((Q+p)/(Q-p));
245
246
247     dsdkdt_value = Z2 * (r02/(8*pi*137)) * (1/k) * (p/p0) *
248       ( (8 * (sintheta2*(2*E02+1))/(p02*delta04)) -
249         ((2*(5*E02+2*E*E0+3))/(p02 * delta02)) -
250         ((2*(p02-k2))/((Q2*delta02))) +
251         ((4*E)/(p02*delta0)) +
252         (L/(p*p0))*(
253                 ((4*E0*sintheta2*(3*k-p02*E))/(p02*delta04)) +
254                 ((4*E02*(E02+E2))/(p02*delta02)) +
255                 ((2-2*(7*E02-3*E*E0+E2))/(p02*delta02)) +
256                 (2*k*(E02+E*E0-1))/((p02*delta0))
257                 ) -
258         ((4*epsilon)/(p*delta0)) +
259         ((epsilonQ)/(p*Q))*
260         (4/delta02-(6*k/delta0)-(2*k*(p02-k2))/(Q2*delta0))
261         );
262
263
264
265     dsdkdt_value = dsdkdt_value*std::sin(theta);
266     return dsdkdt_value;
267}
268
269void G4Generator2BN::ConstructMajorantSurface()
270{
271
272  G4double Eel;
273  G4int vmax;
274  G4double Ek;
275  G4double k, theta, k0, theta0, thetamax;
276  G4double ds, df, dsmax, dk, dt;
277  G4double ratmin;
278  G4double ratio = 0;
279  G4double Vds, Vdf;
280  G4double A, c;
281
282  G4cout << "**** Constructing Majorant Surface for 2BN Distribution ****" << G4endl;
283
284  if(kcut > kmin) kmin = kcut;
285
286  G4int i = 0;
287  // Loop on energy index
288  for(G4int index = index_min; index < index_max; index++){
289
290  G4double fraction = index/100.;
291  Ek = std::pow(10.,fraction);
292  Eel = Ek + electron_mass_c2;
293
294  // find x-section maximum at k=kmin
295  dsmax = 0.;
296  thetamax = 0.;
297
298  for(theta = 0.; theta < pi; theta = theta + dtheta){
299
300    ds = Calculatedsdkdt(kmin, theta, Eel);
301
302    if(ds > dsmax){
303      dsmax = ds;
304      thetamax = theta;
305    }
306  }
307
308  // Compute surface paremeters at kmin
309  if(Ek < kmin || thetamax == 0){
310    c = 0;
311    A = 0;
312  }else{
313    c = 1/(thetamax*thetamax);
314    A = 2*std::sqrt(c)*dsmax/(std::pow(kmin,-b));
315  }
316
317  // look for correction factor to normalization at kmin
318  ratmin = 1.;
319
320  // Volume under surfaces
321  Vds = 0.;
322  Vdf = 0.;
323  k0 = 0.;
324  theta0 = 0.;
325
326  vmax = G4int(100.*std::log10(Ek/kmin));
327
328  for(G4int v = 0; v < vmax; v++){
329    G4double fraction = (v/100.);
330    k = std::pow(10.,fraction)*kmin;
331
332    for(theta = 0.; theta < pi; theta = theta + dtheta){
333      dk = k - k0;
334      dt = theta - theta0;
335      ds = Calculatedsdkdt(k,theta, Eel);
336      Vds = Vds + ds*dk*dt;
337      df = CalculateFkt(k,theta, A, c);
338      Vdf = Vdf + df*dk*dt;
339
340      if(ds != 0.){
341        if(df != 0.) ratio = df/ds;
342      }
343
344      if(ratio < ratmin && ratio != 0.){
345        ratmin = ratio;
346      }
347    }
348  }
349
350
351  // sampling efficiency
352  Atab[i] = A/ratmin * 1.04;
353  ctab[i] = c;
354
355//  G4cout << Ek << " " << i << " " << index << " " << Atab[i] << " " << ctab[i] << " " << G4endl;
356  i++;
357  }
358
359}
360
361G4double G4Generator2BN::Generate2BN(G4double Ek, G4double k) const 
362{
363
364  G4double Eel;
365  G4double kmin2; 
366  G4double kmax, t;
367  G4double cte2;
368  G4double y, u;
369  G4double fk, ft;
370  G4double ds;
371  G4double A2;
372  G4double A, c;
373
374  G4int trials = 0;
375  G4int index;
376
377  // find table index
378  index = G4int(std::log10(Ek)*100) - index_min;
379  Eel = Ek + electron_mass_c2;
380
381  kmax = Ek;
382  kmin2 = kcut;
383
384  c = ctab[index];
385  A = Atab[index];
386  if(index < index_max){
387    A2 = Atab[index+1];
388    if(A2 > A) A = A2;
389  }
390
391  do{
392  // generate k accordimg to std::pow(k,-b)
393  trials++;
394
395  // normalization constant
396//  cte1 = (1-b)/(std::pow(kmax,(1-b))-std::pow(kmin2,(1-b)));
397//  y = G4UniformRand();
398//  k = std::pow(((1-b)/cte1*y+std::pow(kmin2,(1-b))),(1/(1-b)));
399
400  // generate theta accordimg to theta/(1+c*std::pow(theta,2)
401  // Normalization constant
402  cte2 = 2*c/std::log(1+c*pi2);
403
404  y = G4UniformRand();
405  t = std::sqrt((std::exp(2*c*y/cte2)-1)/c);
406  u = G4UniformRand();
407
408  // point acceptance algorithm
409  fk = std::pow(k,-b);
410  ft = t/(1+c*t*t);
411  ds = Calculatedsdkdt(k,t,Eel);
412
413  // violation point
414  if(ds > (A*fk*ft)) G4cout << "WARNING IN G4Generator2BN !!!" << Ek << " " <<  (ds-A*fk*ft)/ds << G4endl;
415
416  }while(u*A*fk*ft > ds);
417
418  return t;
419
420}
421
422void G4Generator2BN::PrintGeneratorInformation() const
423{
424  G4cout << "\n" << G4endl;
425  G4cout << "Bremsstrahlung Angular Generator is 2BN Generator from 2BN Koch & Motz distribution (Rev Mod Phys 31(4), 920 (1959))" << G4endl;
426  G4cout << "\n" << G4endl;
427} 
428
Note: See TracBrowser for help on using the repository browser.