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

Last change on this file since 1196 was 819, checked in by garnier, 16 years ago

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