physics/src/G4RDGenerator2BN.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 // -------------------------------------------------------------------
0028 //
0029 // GEANT4 Class file
0030 //
0031 //
0032 // File name:     G4RDGenerator2BN
0033 //
0034 // Author:        Andreia Trindade (andreia@lip.pt)
0035 //                Pedro Rodrigues  (psilva@lip.pt)
0036 //                Luis Peralta     (luis@lip.pt)
0037 //
0038 // Creation date: 21 June 2003
0039 //
0040 // Modifications:
0041 // 21 Jun 2003                                 First implementation acording with new design
0042 // 05 Nov 2003  MGP                            Fixed compilation warning
0043 // 14 Mar 2004                                 Code optimization
0044 //
0045 // Class Description:
0046 //
0047 // Concrete base class for Bremsstrahlung Angular Distribution Generation - 2BN Distribution
0048 //
0049 // Class Description: End
0050 //
0051 // -------------------------------------------------------------------
0052 //
0053 //
0054
0055 #include "G4RDGenerator2BN.hh"
0056 #include "Randomize.hh"
0057 #include "G4PhysicalConstants.hh"
0058 #include "G4SystemOfUnits.hh"
0059 //
0060
0061 G4double G4RDGenerator2BN::Atab[320] =
0062          { 0.0264767, 0.0260006, 0.0257112, 0.0252475, 0.024792, 0.0243443, 0.0240726, 0.0236367,
0063            0.023209,  0.0227892, 0.0225362, 0.0221284, 0.0217282,0.0214894, 0.0211005, 0.0207189,
0064            0.0204935, 0.0201227, 0.0197588, 0.019546,  0.0191923,0.0188454, 0.0186445, 0.0183072,
0065            0.0179763, 0.0177866, 0.0174649, 0.0172828, 0.0169702,0.0166634, 0.0164915, 0.0161933,
0066            0.0160283, 0.0157384, 0.0155801, 0.0152981, 0.0151463,0.0148721, 0.0147263, 0.0144598,
0067            0.01432,   0.0140607, 0.0139267, 0.0136744, 0.0135459,0.0133005, 0.0131773, 0.0129385,
0068            0.0128205, 0.0125881, 0.012475,  0.0122488, 0.0121406,0.0119204, 0.0118167, 0.0117158,
0069            0.0115032, 0.0114067, 0.0111995, 0.0111072, 0.0110175,0.0108173, 0.0107316, 0.0105365,
0070            0.0104547, 0.0102646, 0.0101865, 0.010111,  0.00992684,0.0098548,0.00967532,0.00960671,
0071            0.00943171,0.00936643,0.00930328,0.0091337, 0.00907372,0.00890831,0.00885141,0.00869003,
0072            0.00863611,0.00858428,0.00842757,0.00837854,0.0082256,0.00817931,0.00803,0.00798639,
0073            0.00784058,0.00779958,0.00776046,0.00761866,0.00758201,0.00744346,0.00740928,0.00727384,
0074            0.00724201,0.00710969,0.00708004,0.00695074,0.00692333,0.00679688,0.00677166,0.00664801,
0075            0.00662484,0.00650396,0.00648286,0.00636458,0.00634545,0.00622977,0.00621258,0.00609936,
0076            0.00608412,0.00597331,0.00595991,0.00585143,0.00583988,0.0057337,0.0057239,0.00561991,
0077            0.0056119, 0.00551005,0.00550377,0.00540399,0.00539938,0.00530162,0.00529872,0.00520292,
0078            0.0051091, 0.00510777,0.00501582,0.00501608,0.00492594,0.00492781,0.00483942,0.0048429,
0079            0.00475622,0.00476127,0.00467625,0.00468287,0.00459947,0.00451785,0.00452581,0.00444573,
0080            0.00445522,0.00437664,0.00438768,0.00431057,0.00432316,0.00424745,0.0042616,0.00418726,
0081            0.004203,  0.00413,   0.00405869,0.00407563,0.00400561,0.00402414,0.00395536,0.00397553,
0082            0.00390795,0.00392975,0.00386339,0.00379862,0.00382167,0.00375805,0.00378276,0.00372031,
0083            0.00374678,0.00368538,0.00371363,0.00365335,0.00359463,0.0036242, 0.00356653,0.003598,
0084            0.00354139,0.00357481,0.00351921,0.00355464,0.00350005,0.0034471,0.00348403,0.00343208,
0085            0.0034712, 0.00342026,0.00346165,0.00341172,0.00345548,0.00340657,0.00335944,0.00340491,
0086            0.00335885,0.00340692,0.00336191,0.00341273,0.00336879,0.00342249,0.00337962,0.00333889,
0087            0.00339463,0.00335506,0.00341401,0.00337558,0.00343797,0.00340067,0.00336584,0.00343059,
0088            0.0033969, 0.00346557,0.00343302,0.00350594,0.00347448,0.00344563,0.00352163,0.00349383,
0089            0.00357485,0.00354807,0.00352395,0.00360885,0.00358571,0.00367661,0.00365446,0.00375194,
0090            0.00373078,0.00371234,0.00381532,0.00379787,0.00390882,0.00389241,0.00387881,0.00399675,
0091            0.00398425,0.00411183,0.00410042,0.00409197,0.00422843,0.00422123,0.00436974,0.0043637,
0092            0.00436082,0.00452075,0.00451934,0.00452125,0.00469406,0.00469756,0.00488741,0.00489221,
0093            0.00490102,0.00510782,0.00511801,0.00513271,0.0053589, 0.00537524,0.00562643,0.00564452,
0094            0.0056677, 0.00594482,0.00596999,0.0059999, 0.00630758,0.00634014,0.00637849,0.00672136,
0095            0.00676236,0.00680914,0.00719407,0.0072439, 0.00730063,0.0077349, 0.00779494,0.00786293,
0096            0.00835577,0.0084276, 0.00850759,0.00907162,0.00915592,0.00924925,0.00935226,0.00999779,
0097            0.0101059, 0.0102249, 0.0109763, 0.0111003, 0.0112367, 0.0113862, 0.0122637, 0.0124196,
0098            0.0125898, 0.0136311, 0.0138081, 0.0140011, 0.0142112, 0.0154536, 0.0156723, 0.0159099,
0099            0.016168,  0.0176664, 0.0179339, 0.0182246, 0.0185396, 0.020381,  0.0207026, 0.0210558,
0100            0.0214374, 0.0237377, 0.0241275, 0.0245528, 0.0250106, 0.0255038, 0.0284158, 0.0289213,
0101            0.0294621, 0.0300526, 0.0338619, 0.0344537, 0.0351108, 0.0358099, 0.036554,  0.0416399
0102          };
0103
0104 G4double G4RDGenerator2BN::ctab[320] =
0105     { 0.482253, 0.482253, 0.489021, 0.489021, 0.489021, 0.489021, 0.495933,
0106       0.495933, 0.495933, 0.495933, 0.502993, 0.502993, 0.502993, 0.510204,
0107       0.510204, 0.510204, 0.517572, 0.517572, 0.517572, 0.5251,   0.5251,
0108       0.5251,   0.532793, 0.532793, 0.532793, 0.540657, 0.540657, 0.548697,
0109       0.548697, 0.548697, 0.556917, 0.556917, 0.565323, 0.565323, 0.573921,
0110       0.573921, 0.582717, 0.582717, 0.591716, 0.591716, 0.600925, 0.600925,
0111       0.610352, 0.610352, 0.620001, 0.620001, 0.629882, 0.629882, 0.64,
0112       0.64,     0.650364, 0.650364, 0.660982, 0.660982, 0.671862, 0.683013,
0113       0.683013, 0.694444, 0.694444, 0.706165, 0.718184, 0.718184, 0.730514,
0114       0.730514, 0.743163, 0.743163, 0.756144, 0.769468, 0.769468, 0.783147,
0115       0.783147, 0.797194, 0.797194, 0.811622, 0.826446, 0.826446, 0.84168,
0116       0.84168,  0.857339, 0.857339, 0.873439, 0.889996, 0.889996, 0.907029,
0117       0.907029, 0.924556, 0.924556, 0.942596, 0.942596, 0.961169, 0.980296,
0118       0.980296, 1,        1,        1.0203,   1.0203,   1.04123,  1.04123,
0119       1.06281,  1.06281,  1.08507,  1.08507,  1.10803,  1.10803,  1.13173,
0120       1.13173,  1.1562,   1.1562,   1.18147,  1.18147,  1.20758,  1.20758,
0121       1.23457,  1.23457,  1.26247,  1.26247,  1.29132,  1.29132,  1.32118,
0122       1.32118,  1.35208,  1.35208,  1.38408,  1.38408,  1.41723,  1.41723,
0123       1.45159,  1.45159,  1.45159,  1.48721,  1.48721,  1.52416,  1.52416,
0124       1.5625,   1.5625,   1.60231,  1.60231,  1.64366,  1.64366,  1.68663,
0125       1.68663,  1.68663,  1.7313,   1.7313,   1.77778,  1.77778,  1.82615,
0126       1.82615,  1.87652,  1.87652,  1.92901,  1.92901,  1.98373,  1.98373,
0127       1.98373,  2.04082,  2.04082,  2.1004,   2.1004,   2.16263,  2.16263,
0128       2.22767,  2.22767,  2.22767,  2.29568,  2.29568,  2.36686,  2.36686,
0129       2.44141,  2.44141,  2.51953,  2.51953,  2.51953,  2.60146,  2.60146,
0130       2.68745,  2.68745,  2.77778,  2.77778,  2.87274,  2.87274,  2.87274,
0131       2.97265,  2.97265,  3.07787,  3.07787,  3.18878,  3.18878,  3.30579,
0132       3.30579,  3.30579,  3.42936,  3.42936,  3.55999,  3.55999,  3.69822,
0133       3.69822,  3.84468,  3.84468,  3.84468,  4,        4,        4.16493,
0134       4.16493,  4.34028,  4.34028,  4.34028,  4.52694,  4.52694,  4.7259,
0135       4.7259,   4.93827,  4.93827,  4.93827,  5.16529,  5.16529,  5.40833,
0136       5.40833,  5.40833,  5.66893,  5.66893,  5.94884,  5.94884,  6.25,
0137       6.25,     6.25,     6.57462,  6.57462,  6.92521,  6.92521,  6.92521,
0138       7.3046,   7.3046,   7.71605,  7.71605,  7.71605,  8.16327,  8.16327,
0139       8.65052,  8.65052,  8.65052,  9.18274,  9.18274,  9.18274,  9.76562,
0140       9.76562,  10.4058,  10.4058,  10.4058,  11.1111,  11.1111,  11.1111,
0141       11.8906,  11.8906,  12.7551,  12.7551,  12.7551,  13.7174,  13.7174,
0142       13.7174,  14.7929,  14.7929,  14.7929,  16,       16,       16,
0143       17.3611,  17.3611,  17.3611,  18.9036,  18.9036,  18.9036,  20.6612,
0144       20.6612,  20.6612,  22.6757,  22.6757,  22.6757,  22.6757,  25,
0145       25,       25,       27.7008,  27.7008,  27.7008,  27.7008,  30.8642,
0146       30.8642,  30.8642,  34.6021,  34.6021,  34.6021,  34.6021,  39.0625,
0147       39.0625,  39.0625,  39.0625,  44.4444,  44.4444,  44.4444,  44.4444,
0148       51.0204,  51.0204,  51.0204,  51.0204,  59.1716,  59.1716,  59.1716,
0149       59.1716,  59.1716,  69.4444,  69.4444,  69.4444,  69.4444,  82.6446,
0150       82.6446,  82.6446,  82.6446,  82.6446,  100
0151     };
0152
0153
0154 G4RDGenerator2BN::G4RDGenerator2BN(const G4String& name):G4RDVBremAngularDistribution(name)
0155 {
0156   b = 1.2;
0157   index_min = -300;
0158   index_max = 20;
0159
0160   // Set parameters minimum limits Ekelectron = 250 eV and kphoton = 100 eV
0161   kmin = 100*eV;
0162   Ekmin = 250*eV;
0163   kcut = 100*eV;
0164
0165   // Increment Theta value for surface interpolation
0166   dtheta = 0.1*rad;
0167
0168   // Construct Majorant Surface to 2BN cross-section
0169   // (we dont need this. Pre-calculated values are used instead due to performance issues
0170   // ConstructMajorantSurface();
0171 }
0172
0173 //
0174
0175 G4RDGenerator2BN::~G4RDGenerator2BN()
0176 {;}
0177
0178 //
0179
0180 G4double G4RDGenerator2BN::PolarAngle(const G4double initial_energy,
0181                     const G4double final_energy,
0182                     const G4int) // Z
0183 {
0184
0185   G4double theta = 0;
0186
0187   G4double k = initial_energy - final_energy;
0188
0189   theta = Generate2BN(initial_energy, k);
0190
0191   return theta;
0192 }
0193 //
0194
0195 G4double G4RDGenerator2BN::CalculateFkt(G4double k, G4double theta, G4double A, G4double c) const
0196 {
0197   G4double Fkt_value = 0;
0198   Fkt_value = A*std::pow(k,-b)*theta/(1+c*theta*theta);
0199   return Fkt_value;
0200 }
0201
0202 G4double G4RDGenerator2BN::Calculatedsdkdt(G4double kout, G4double theta, G4double Eel) const
0203 {
0204
0205   G4double dsdkdt_value = 0.;
0206   G4double Z = 1;
0207   // classic radius (in cm)
0208   G4double r0 = 2.82E-13;
0209   //squared classic radius (in barn)
0210   G4double r02 = r0*r0*1.0E+24;
0211
0212   // Photon energy cannot be greater than electron kinetic energy
0213   if(kout > (Eel-electron_mass_c2)){
0214     dsdkdt_value = 0;
0215     return dsdkdt_value;
0216   }
0217
0218      G4double E0 = Eel/electron_mass_c2;
0219      G4double k =  kout/electron_mass_c2;
0220      G4double E =  E0 - k;
0221
0222      if(E <= 1*MeV ){                                 // Kinematic limit at 1 MeV !!
0223        dsdkdt_value = 0;
0224        return dsdkdt_value;
0225      }
0226
0227
0228      G4double p0 = std::sqrt(E0*E0-1);
0229      G4double p = std::sqrt(E*E-1);
0230      G4double L = std::log((E*E0-1+p*p0)/(E*E0-1-p*p0));
0231      G4double delta0 = E0 - p0*std::cos(theta);
0232      G4double epsilon = std::log((E+p)/(E-p));
0233      G4double Z2 = Z*Z;
0234      G4double sintheta2 = std::sin(theta)*std::sin(theta);
0235      G4double E02 = E0*E0;
0236      G4double E2 = E*E;
0237      G4double p02 = E0*E0-1;
0238      G4double k2 = k*k;
0239      G4double delta02 = delta0*delta0;
0240      G4double delta04 =  delta02* delta02;
0241      G4double Q = std::sqrt(p02+k2-2*k*p0*std::cos(theta));
0242      G4double Q2 = Q*Q;
0243      G4double epsilonQ = std::log((Q+p)/(Q-p));
0244
0245
0246      dsdkdt_value = Z2 * (r02/(8*pi*137)) * (1/k) * (p/p0) *
0247        ( (8 * (sintheta2*(2*E02+1))/(p02*delta04)) -
0248          ((2*(5*E02+2*E*E0+3))/(p02 * delta02)) -
0249          ((2*(p02-k2))/((Q2*delta02))) +
0250          ((4*E)/(p02*delta0)) +
0251          (L/(p*p0))*(
0252                  ((4*E0*sintheta2*(3*k-p02*E))/(p02*delta04)) +
0253                  ((4*E02*(E02+E2))/(p02*delta02)) +
0254                  ((2-2*(7*E02-3*E*E0+E2))/(p02*delta02)) +
0255                  (2*k*(E02+E*E0-1))/((p02*delta0))
0256                  ) -
0257          ((4*epsilon)/(p*delta0)) +
0258          ((epsilonQ)/(p*Q))*
0259          (4/delta02-(6*k/delta0)-(2*k*(p02-k2))/(Q2*delta0))
0260          );
0261
0262
0263
0264      dsdkdt_value = dsdkdt_value*std::sin(theta);
0265      return dsdkdt_value;
0266 }
0267
0268 void G4RDGenerator2BN::ConstructMajorantSurface()
0269 {
0270
0271   G4double Eel;
0272   G4int vmax;
0273   G4double Ek;
0274   G4double k, theta, k0, theta0, thetamax;
0275   G4double ds, df, dsmax, dk, dt;
0276   G4double ratmin;
0277   G4double ratio = 0;
0278   G4double Vds, Vdf;
0279   G4double A, c;
0280
0281   G4cout << "**** Constructing Majorant Surface for 2BN Distribution ****" << G4endl;
0282
0283   if(kcut > kmin) kmin = kcut;
0284
0285   G4int i = 0;
0286   // Loop on energy index
0287   for(G4int index = index_min; index < index_max; index++){
0288
0289   G4double fraction = index/100.;
0290   Ek = std::pow(10.,fraction);
0291   Eel = Ek + electron_mass_c2;
0292
0293   // find x-section maximum at k=kmin
0294   dsmax = 0.;
0295   thetamax = 0.;
0296
0297   for(theta = 0.; theta < pi; theta = theta + dtheta){
0298
0299     ds = Calculatedsdkdt(kmin, theta, Eel);
0300
0301     if(ds > dsmax){
0302       dsmax = ds;
0303       thetamax = theta;
0304     }
0305   }
0306
0307   // Compute surface paremeters at kmin
0308   if(Ek < kmin || thetamax == 0){
0309     c = 0;
0310     A = 0;
0311   }else{
0312     c = 1/(thetamax*thetamax);
0313     A = 2*std::sqrt(c)*dsmax/(std::pow(kmin,-b));
0314   }
0315
0316   // look for correction factor to normalization at kmin
0317   ratmin = 1.;
0318
0319   // Volume under surfaces
0320   Vds = 0.;
0321   Vdf = 0.;
0322   k0 = 0.;
0323   theta0 = 0.;
0324
0325   vmax = G4int(100.*std::log10(Ek/kmin));
0326
0327   for(G4int v = 0; v < vmax; v++){
0328     G4double fraction = (v/100.);
0329     k = std::pow(10.,fraction)*kmin;
0330
0331     for(theta = 0.; theta < pi; theta = theta + dtheta){
0332       dk = k - k0;
0333       dt = theta - theta0;
0334       ds = Calculatedsdkdt(k,theta, Eel);
0335       Vds = Vds + ds*dk*dt;
0336       df = CalculateFkt(k,theta, A, c);
0337       Vdf = Vdf + df*dk*dt;
0338
0339       if(ds != 0.){
0340     if(df != 0.) ratio = df/ds;
0341       }
0342
0343       if(ratio < ratmin && ratio != 0.){
0344     ratmin = ratio;
0345       }
0346     }
0347   }
0348
0349
0350   // sampling efficiency
0351   Atab[i] = A/ratmin * 1.04;
0352   ctab[i] = c;
0353
0354 //  G4cout << Ek << " " << i << " " << index << " " << Atab[i] << " " << ctab[i] << " " << G4endl;
0355   i++;
0356   }
0357
0358 }
0359
0360 G4double G4RDGenerator2BN::Generate2BN(G4double Ek, G4double k) const
0361 {
0362
0363   G4double Eel;
0364   G4double t;
0365   G4double cte2;
0366   G4double y, u;
0367   G4double fk, ft;
0368   G4double ds;
0369   G4double A2;
0370   G4double A, c;
0371
0372   G4int trials = 0;
0373   G4int index;
0374
0375   // find table index
0376   index = G4int(std::log10(Ek)*100) - index_min;
0377   Eel = Ek + electron_mass_c2;
0378
0379   c = ctab[index];
0380   A = Atab[index];
0381   if(index < index_max){
0382     A2 = Atab[index+1];
0383     if(A2 > A) A = A2;
0384   }
0385
0386   do{
0387   // generate k accordimg to std::pow(k,-b)
0388   trials++;
0389
0390   // normalization constant
0391 //  cte1 = (1-b)/(std::pow(kmax,(1-b))-std::pow(kmin2,(1-b)));
0392 //  y = G4UniformRand();
0393 //  k = std::pow(((1-b)/cte1*y+std::pow(kmin2,(1-b))),(1/(1-b)));
0394
0395   // generate theta accordimg to theta/(1+c*std::pow(theta,2)
0396   // Normalization constant
0397   cte2 = 2*c/std::log(1+c*pi2);
0398
0399   y = G4UniformRand();
0400   t = std::sqrt((std::exp(2*c*y/cte2)-1)/c);
0401   u = G4UniformRand();
0402
0403   // point acceptance algorithm
0404   fk = std::pow(k,-b);
0405   ft = t/(1+c*t*t);
0406   ds = Calculatedsdkdt(k,t,Eel);
0407
0408   // violation point
0409   if(ds > (A*fk*ft)) G4cout << "WARNING IN G4RDGenerator2BN !!!" << Ek << " " <<  (ds-A*fk*ft)/ds << G4endl;
0410
0411   }while(u*A*fk*ft > ds);
0412
0413   return t;
0414
0415 }
0416
0417 void G4RDGenerator2BN::PrintGeneratorInformation() const
0418 {
0419   G4cout << "\n" << G4endl;
0420   G4cout << "Bremsstrahlung Angular Generator is 2BN Generator from 2BN Koch & Motz distribution (Rev Mod Phys 31(4), 920 (1959))" << G4endl;
0421   G4cout << "\n" << G4endl;
0422 }
0423