00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056 #include "G4Generator2BN.hh"
00057 #include "Randomize.hh"
00058 #include "G4PhysicalConstants.hh"
00059 #include "G4SystemOfUnits.hh"
00060
00061
00062
00063 G4double G4Generator2BN::Atab[320] =
00064 { 0.0264767, 0.0260006, 0.0257112, 0.0252475, 0.024792, 0.0243443, 0.0240726, 0.0236367,
00065 0.023209, 0.0227892, 0.0225362, 0.0221284, 0.0217282,0.0214894, 0.0211005, 0.0207189,
00066 0.0204935, 0.0201227, 0.0197588, 0.019546, 0.0191923,0.0188454, 0.0186445, 0.0183072,
00067 0.0179763, 0.0177866, 0.0174649, 0.0172828, 0.0169702,0.0166634, 0.0164915, 0.0161933,
00068 0.0160283, 0.0157384, 0.0155801, 0.0152981, 0.0151463,0.0148721, 0.0147263, 0.0144598,
00069 0.01432, 0.0140607, 0.0139267, 0.0136744, 0.0135459,0.0133005, 0.0131773, 0.0129385,
00070 0.0128205, 0.0125881, 0.012475, 0.0122488, 0.0121406,0.0119204, 0.0118167, 0.0117158,
00071 0.0115032, 0.0114067, 0.0111995, 0.0111072, 0.0110175,0.0108173, 0.0107316, 0.0105365,
00072 0.0104547, 0.0102646, 0.0101865, 0.010111, 0.00992684,0.0098548,0.00967532,0.00960671,
00073 0.00943171,0.00936643,0.00930328,0.0091337, 0.00907372,0.00890831,0.00885141,0.00869003,
00074 0.00863611,0.00858428,0.00842757,0.00837854,0.0082256,0.00817931,0.00803,0.00798639,
00075 0.00784058,0.00779958,0.00776046,0.00761866,0.00758201,0.00744346,0.00740928,0.00727384,
00076 0.00724201,0.00710969,0.00708004,0.00695074,0.00692333,0.00679688,0.00677166,0.00664801,
00077 0.00662484,0.00650396,0.00648286,0.00636458,0.00634545,0.00622977,0.00621258,0.00609936,
00078 0.00608412,0.00597331,0.00595991,0.00585143,0.00583988,0.0057337,0.0057239,0.00561991,
00079 0.0056119, 0.00551005,0.00550377,0.00540399,0.00539938,0.00530162,0.00529872,0.00520292,
00080 0.0051091, 0.00510777,0.00501582,0.00501608,0.00492594,0.00492781,0.00483942,0.0048429,
00081 0.00475622,0.00476127,0.00467625,0.00468287,0.00459947,0.00451785,0.00452581,0.00444573,
00082 0.00445522,0.00437664,0.00438768,0.00431057,0.00432316,0.00424745,0.0042616,0.00418726,
00083 0.004203, 0.00413, 0.00405869,0.00407563,0.00400561,0.00402414,0.00395536,0.00397553,
00084 0.00390795,0.00392975,0.00386339,0.00379862,0.00382167,0.00375805,0.00378276,0.00372031,
00085 0.00374678,0.00368538,0.00371363,0.00365335,0.00359463,0.0036242, 0.00356653,0.003598,
00086 0.00354139,0.00357481,0.00351921,0.00355464,0.00350005,0.0034471,0.00348403,0.00343208,
00087 0.0034712, 0.00342026,0.00346165,0.00341172,0.00345548,0.00340657,0.00335944,0.00340491,
00088 0.00335885,0.00340692,0.00336191,0.00341273,0.00336879,0.00342249,0.00337962,0.00333889,
00089 0.00339463,0.00335506,0.00341401,0.00337558,0.00343797,0.00340067,0.00336584,0.00343059,
00090 0.0033969, 0.00346557,0.00343302,0.00350594,0.00347448,0.00344563,0.00352163,0.00349383,
00091 0.00357485,0.00354807,0.00352395,0.00360885,0.00358571,0.00367661,0.00365446,0.00375194,
00092 0.00373078,0.00371234,0.00381532,0.00379787,0.00390882,0.00389241,0.00387881,0.00399675,
00093 0.00398425,0.00411183,0.00410042,0.00409197,0.00422843,0.00422123,0.00436974,0.0043637,
00094 0.00436082,0.00452075,0.00451934,0.00452125,0.00469406,0.00469756,0.00488741,0.00489221,
00095 0.00490102,0.00510782,0.00511801,0.00513271,0.0053589, 0.00537524,0.00562643,0.00564452,
00096 0.0056677, 0.00594482,0.00596999,0.0059999, 0.00630758,0.00634014,0.00637849,0.00672136,
00097 0.00676236,0.00680914,0.00719407,0.0072439, 0.00730063,0.0077349, 0.00779494,0.00786293,
00098 0.00835577,0.0084276, 0.00850759,0.00907162,0.00915592,0.00924925,0.00935226,0.00999779,
00099 0.0101059, 0.0102249, 0.0109763, 0.0111003, 0.0112367, 0.0113862, 0.0122637, 0.0124196,
00100 0.0125898, 0.0136311, 0.0138081, 0.0140011, 0.0142112, 0.0154536, 0.0156723, 0.0159099,
00101 0.016168, 0.0176664, 0.0179339, 0.0182246, 0.0185396, 0.020381, 0.0207026, 0.0210558,
00102 0.0214374, 0.0237377, 0.0241275, 0.0245528, 0.0250106, 0.0255038, 0.0284158, 0.0289213,
00103 0.0294621, 0.0300526, 0.0338619, 0.0344537, 0.0351108, 0.0358099, 0.036554, 0.0416399
00104 };
00105
00106 G4double G4Generator2BN::ctab[320] =
00107 { 0.482253, 0.482253, 0.489021, 0.489021, 0.489021, 0.489021, 0.495933,
00108 0.495933, 0.495933, 0.495933, 0.502993, 0.502993, 0.502993, 0.510204,
00109 0.510204, 0.510204, 0.517572, 0.517572, 0.517572, 0.5251, 0.5251,
00110 0.5251, 0.532793, 0.532793, 0.532793, 0.540657, 0.540657, 0.548697,
00111 0.548697, 0.548697, 0.556917, 0.556917, 0.565323, 0.565323, 0.573921,
00112 0.573921, 0.582717, 0.582717, 0.591716, 0.591716, 0.600925, 0.600925,
00113 0.610352, 0.610352, 0.620001, 0.620001, 0.629882, 0.629882, 0.64,
00114 0.64, 0.650364, 0.650364, 0.660982, 0.660982, 0.671862, 0.683013,
00115 0.683013, 0.694444, 0.694444, 0.706165, 0.718184, 0.718184, 0.730514,
00116 0.730514, 0.743163, 0.743163, 0.756144, 0.769468, 0.769468, 0.783147,
00117 0.783147, 0.797194, 0.797194, 0.811622, 0.826446, 0.826446, 0.84168,
00118 0.84168, 0.857339, 0.857339, 0.873439, 0.889996, 0.889996, 0.907029,
00119 0.907029, 0.924556, 0.924556, 0.942596, 0.942596, 0.961169, 0.980296,
00120 0.980296, 1, 1, 1.0203, 1.0203, 1.04123, 1.04123,
00121 1.06281, 1.06281, 1.08507, 1.08507, 1.10803, 1.10803, 1.13173,
00122 1.13173, 1.1562, 1.1562, 1.18147, 1.18147, 1.20758, 1.20758,
00123 1.23457, 1.23457, 1.26247, 1.26247, 1.29132, 1.29132, 1.32118,
00124 1.32118, 1.35208, 1.35208, 1.38408, 1.38408, 1.41723, 1.41723,
00125 1.45159, 1.45159, 1.45159, 1.48721, 1.48721, 1.52416, 1.52416,
00126 1.5625, 1.5625, 1.60231, 1.60231, 1.64366, 1.64366, 1.68663,
00127 1.68663, 1.68663, 1.7313, 1.7313, 1.77778, 1.77778, 1.82615,
00128 1.82615, 1.87652, 1.87652, 1.92901, 1.92901, 1.98373, 1.98373,
00129 1.98373, 2.04082, 2.04082, 2.1004, 2.1004, 2.16263, 2.16263,
00130 2.22767, 2.22767, 2.22767, 2.29568, 2.29568, 2.36686, 2.36686,
00131 2.44141, 2.44141, 2.51953, 2.51953, 2.51953, 2.60146, 2.60146,
00132 2.68745, 2.68745, 2.77778, 2.77778, 2.87274, 2.87274, 2.87274,
00133 2.97265, 2.97265, 3.07787, 3.07787, 3.18878, 3.18878, 3.30579,
00134 3.30579, 3.30579, 3.42936, 3.42936, 3.55999, 3.55999, 3.69822,
00135 3.69822, 3.84468, 3.84468, 3.84468, 4, 4, 4.16493,
00136 4.16493, 4.34028, 4.34028, 4.34028, 4.52694, 4.52694, 4.7259,
00137 4.7259, 4.93827, 4.93827, 4.93827, 5.16529, 5.16529, 5.40833,
00138 5.40833, 5.40833, 5.66893, 5.66893, 5.94884, 5.94884, 6.25,
00139 6.25, 6.25, 6.57462, 6.57462, 6.92521, 6.92521, 6.92521,
00140 7.3046, 7.3046, 7.71605, 7.71605, 7.71605, 8.16327, 8.16327,
00141 8.65052, 8.65052, 8.65052, 9.18274, 9.18274, 9.18274, 9.76562,
00142 9.76562, 10.4058, 10.4058, 10.4058, 11.1111, 11.1111, 11.1111,
00143 11.8906, 11.8906, 12.7551, 12.7551, 12.7551, 13.7174, 13.7174,
00144 13.7174, 14.7929, 14.7929, 14.7929, 16, 16, 16,
00145 17.3611, 17.3611, 17.3611, 18.9036, 18.9036, 18.9036, 20.6612,
00146 20.6612, 20.6612, 22.6757, 22.6757, 22.6757, 22.6757, 25,
00147 25, 25, 27.7008, 27.7008, 27.7008, 27.7008, 30.8642,
00148 30.8642, 30.8642, 34.6021, 34.6021, 34.6021, 34.6021, 39.0625,
00149 39.0625, 39.0625, 39.0625, 44.4444, 44.4444, 44.4444, 44.4444,
00150 51.0204, 51.0204, 51.0204, 51.0204, 59.1716, 59.1716, 59.1716,
00151 59.1716, 59.1716, 69.4444, 69.4444, 69.4444, 69.4444, 82.6446,
00152 82.6446, 82.6446, 82.6446, 82.6446, 100
00153 };
00154
00155
00156 G4Generator2BN::G4Generator2BN(const G4String&)
00157 : G4VEmAngularDistribution("AngularGen2BN")
00158 {
00159 b = 1.2;
00160 index_min = -300;
00161 index_max = 319;
00162
00163
00164 kmin = 100*eV;
00165 Ekmin = 250*eV;
00166 kcut = 100*eV;
00167
00168
00169 dtheta = 0.1*rad;
00170
00171 nwarn = 0;
00172
00173
00174
00175
00176 }
00177
00178
00179
00180 G4Generator2BN::~G4Generator2BN()
00181 {}
00182
00183 G4ThreeVector& G4Generator2BN::SampleDirection(const G4DynamicParticle* dp,
00184 G4double out_energy,
00185 G4int,
00186 const G4Material*)
00187 {
00188 G4double Ek = dp->GetKineticEnergy();
00189 G4double Eel = dp->GetTotalEnergy();
00190 if(Eel > 2*MeV) {
00191 return fGenerator2BS.SampleDirection(dp, out_energy, 0);
00192 }
00193
00194 G4double k = Eel - out_energy;
00195
00196 G4double t;
00197 G4double cte2;
00198 G4double y, u;
00199 G4double ds;
00200 G4double dmax;
00201
00202 G4int trials = 0;
00203
00204
00205 G4int index = G4int(std::log10(Ek/MeV)*100) - index_min;
00206 if(index > index_max) { index = index_max; }
00207 else if(index < 0) { index = 0; }
00208
00209
00210
00211
00212 G4double c = ctab[index];
00213 G4double A = Atab[index];
00214 if(index < index_max) { A = std::max(A, Atab[index+1]); }
00215
00216 do{
00217
00218 trials++;
00219
00220
00221
00222
00223
00224
00225
00226
00227 cte2 = 2*c/std::log(1+c*pi2);
00228
00229 y = G4UniformRand();
00230 t = std::sqrt((std::exp(2*c*y/cte2)-1)/c);
00231 u = G4UniformRand();
00232
00233
00234
00235
00236 dmax = A*std::pow(k,-b)*t/(1+c*t*t);
00237 ds = Calculatedsdkdt(k,t,Eel);
00238
00239
00240 if(ds > dmax && nwarn >= 20) {
00241 ++nwarn;
00242 G4cout << "### WARNING in G4Generator2BN: Ekin(MeV)= " << Ek/MeV
00243 << " D(Ekin,k)/Dmax-1= " << (ds/dmax - 1)
00244 << " results are not reliable!"
00245 << G4endl;
00246 if(20 == nwarn) {
00247 G4cout << " WARNING in G4Generator2BN is closed" << G4endl;
00248 }
00249 }
00250
00251 } while(u*dmax > ds || t > CLHEP::pi);
00252
00253 G4double sint = std::sin(t);
00254 G4double phi = twopi*G4UniformRand();
00255
00256 fLocalDirection.set(sint*std::cos(phi), sint*std::sin(phi),std::cos(t));
00257 fLocalDirection.rotateUz(dp->GetMomentumDirection());
00258
00259 return fLocalDirection;
00260 }
00261
00262 G4double G4Generator2BN::CalculateFkt(G4double k, G4double theta, G4double A, G4double c) const
00263 {
00264 G4double Fkt_value = 0;
00265 Fkt_value = A*std::pow(k,-b)*theta/(1+c*theta*theta);
00266 return Fkt_value;
00267 }
00268
00269 G4double G4Generator2BN::Calculatedsdkdt(G4double kout, G4double theta, G4double Eel) const
00270 {
00271
00272 G4double dsdkdt_value = 0.;
00273 G4double Z = 1;
00274
00275 G4double r0 = 2.82E-13;
00276
00277 G4double r02 = r0*r0*1.0E+24;
00278
00279
00280 if(kout > (Eel-electron_mass_c2)){
00281 dsdkdt_value = 0;
00282 return dsdkdt_value;
00283 }
00284
00285 G4double E0 = Eel/electron_mass_c2;
00286 G4double k = kout/electron_mass_c2;
00287 G4double E = E0 - k;
00288
00289 if(E <= 1*MeV ){
00290 dsdkdt_value = 0;
00291 return dsdkdt_value;
00292 }
00293
00294
00295 G4double p0 = std::sqrt(E0*E0-1);
00296 G4double p = std::sqrt(E*E-1);
00297 G4double L = std::log((E*E0-1+p*p0)/(E*E0-1-p*p0));
00298 G4double delta0 = E0 - p0*std::cos(theta);
00299 G4double epsilon = std::log((E+p)/(E-p));
00300 G4double Z2 = Z*Z;
00301 G4double sintheta2 = std::sin(theta)*std::sin(theta);
00302 G4double E02 = E0*E0;
00303 G4double E2 = E*E;
00304 G4double p02 = E0*E0-1;
00305 G4double k2 = k*k;
00306 G4double delta02 = delta0*delta0;
00307 G4double delta04 = delta02* delta02;
00308 G4double Q = std::sqrt(p02+k2-2*k*p0*std::cos(theta));
00309 G4double Q2 = Q*Q;
00310 G4double epsilonQ = std::log((Q+p)/(Q-p));
00311
00312
00313 dsdkdt_value = Z2 * (r02/(8*pi*137)) * (1/k) * (p/p0) *
00314 ( (8 * (sintheta2*(2*E02+1))/(p02*delta04)) -
00315 ((2*(5*E02+2*E*E0+3))/(p02 * delta02)) -
00316 ((2*(p02-k2))/((Q2*delta02))) +
00317 ((4*E)/(p02*delta0)) +
00318 (L/(p*p0))*(
00319 ((4*E0*sintheta2*(3*k-p02*E))/(p02*delta04)) +
00320 ((4*E02*(E02+E2))/(p02*delta02)) +
00321 ((2-2*(7*E02-3*E*E0+E2))/(p02*delta02)) +
00322 (2*k*(E02+E*E0-1))/((p02*delta0))
00323 ) -
00324 ((4*epsilon)/(p*delta0)) +
00325 ((epsilonQ)/(p*Q))*
00326 (4/delta02-(6*k/delta0)-(2*k*(p02-k2))/(Q2*delta0))
00327 );
00328
00329
00330
00331 dsdkdt_value = dsdkdt_value*std::sin(theta);
00332 return dsdkdt_value;
00333 }
00334
00335 void G4Generator2BN::ConstructMajorantSurface()
00336 {
00337
00338 G4double Eel;
00339 G4int vmax;
00340 G4double Ek;
00341 G4double k, theta, k0, theta0, thetamax;
00342 G4double ds, df, dsmax, dk, dt;
00343 G4double ratmin;
00344 G4double ratio = 0;
00345 G4double Vds, Vdf;
00346 G4double A, c;
00347
00348 G4cout << "**** Constructing Majorant Surface for 2BN Distribution ****" << G4endl;
00349
00350 if(kcut > kmin) kmin = kcut;
00351
00352 G4int i = 0;
00353
00354 for(G4int index = index_min; index < index_max; index++){
00355
00356 G4double fraction = index/100.;
00357 Ek = std::pow(10.,fraction);
00358 Eel = Ek + electron_mass_c2;
00359
00360
00361 dsmax = 0.;
00362 thetamax = 0.;
00363
00364 for(theta = 0.; theta < pi; theta = theta + dtheta){
00365
00366 ds = Calculatedsdkdt(kmin, theta, Eel);
00367
00368 if(ds > dsmax){
00369 dsmax = ds;
00370 thetamax = theta;
00371 }
00372 }
00373
00374
00375 if(Ek < kmin || thetamax == 0){
00376 c = 0;
00377 A = 0;
00378 }else{
00379 c = 1/(thetamax*thetamax);
00380 A = 2*std::sqrt(c)*dsmax/(std::pow(kmin,-b));
00381 }
00382
00383
00384 ratmin = 1.;
00385
00386
00387 Vds = 0.;
00388 Vdf = 0.;
00389 k0 = 0.;
00390 theta0 = 0.;
00391
00392 vmax = G4int(100.*std::log10(Ek/kmin));
00393
00394 for(G4int v = 0; v < vmax; v++){
00395 G4double fractionLocal = (v/100.);
00396 k = std::pow(10.,fractionLocal)*kmin;
00397
00398 for(theta = 0.; theta < pi; theta = theta + dtheta){
00399 dk = k - k0;
00400 dt = theta - theta0;
00401 ds = Calculatedsdkdt(k,theta, Eel);
00402 Vds = Vds + ds*dk*dt;
00403 df = CalculateFkt(k,theta, A, c);
00404 Vdf = Vdf + df*dk*dt;
00405
00406 if(ds != 0.){
00407 if(df != 0.) ratio = df/ds;
00408 }
00409
00410 if(ratio < ratmin && ratio != 0.){
00411 ratmin = ratio;
00412 }
00413 }
00414 }
00415
00416
00417
00418 Atab[i] = A/ratmin * 1.04;
00419 ctab[i] = c;
00420
00421
00422 i++;
00423 }
00424
00425 }
00426
00427 void G4Generator2BN::PrintGeneratorInformation() const
00428 {
00429 G4cout << "\n" << G4endl;
00430 G4cout << "Bremsstrahlung Angular Generator is 2BN Generator from 2BN Koch & Motz distribution (Rev Mod Phys 31(4), 920 (1959))" << G4endl;
00431 G4cout << "\n" << G4endl;
00432 }
00433