Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // ------------------------------------------- 28 // 29 // GEANT4 Class file 30 // 31 // 32 // File name: G4Generator2BN 33 // 34 // Author: Andreia Trindade (andreia@li 35 // Pedro Rodrigues (psilva@lip 36 // Luis Peralta (luis@lip.p 37 // 38 // Creation date: 21 June 2003 39 // 40 // Modifications: 41 // 21 Jun 2003 42 // 05 Nov 2003 MGP 43 // 14 Mar 2004 44 // 45 // Class Description: 46 // 47 // Concrete base class for Bremsstrahlung Angu 48 // 2BN Distribution 49 // 50 // Class Description: End 51 // 52 // ------------------------------------------- 53 // 54 55 #include "G4Generator2BN.hh" 56 #include "Randomize.hh" 57 #include "G4PhysicalConstants.hh" 58 #include "G4SystemOfUnits.hh" 59 #include "G4Exp.hh" 60 61 // 62 63 G4double G4Generator2BN::Atab[320] = 64 { 0.0264767, 0.0260006, 0.0257112, 0. 65 0.023209, 0.0227892, 0.0225362, 0. 66 0.0204935, 0.0201227, 0.0197588, 0. 67 0.0179763, 0.0177866, 0.0174649, 0. 68 0.0160283, 0.0157384, 0.0155801, 0. 69 0.01432, 0.0140607, 0.0139267, 0. 70 0.0128205, 0.0125881, 0.012475, 0. 71 0.0115032, 0.0114067, 0.0111995, 0. 72 0.0104547, 0.0102646, 0.0101865, 0. 73 0.00943171,0.00936643,0.00930328,0. 74 0.00863611,0.00858428,0.00842757,0. 75 0.00784058,0.00779958,0.00776046,0. 76 0.00724201,0.00710969,0.00708004,0. 77 0.00662484,0.00650396,0.00648286,0. 78 0.00608412,0.00597331,0.00595991,0. 79 0.0056119, 0.00551005,0.00550377,0. 80 0.0051091, 0.00510777,0.00501582,0. 81 0.00475622,0.00476127,0.00467625,0. 82 0.00445522,0.00437664,0.00438768,0. 83 0.004203, 0.00413, 0.00405869,0. 84 0.00390795,0.00392975,0.00386339,0. 85 0.00374678,0.00368538,0.00371363,0. 86 0.00354139,0.00357481,0.00351921,0. 87 0.0034712, 0.00342026,0.00346165,0. 88 0.00335885,0.00340692,0.00336191,0. 89 0.00339463,0.00335506,0.00341401,0. 90 0.0033969, 0.00346557,0.00343302,0. 91 0.00357485,0.00354807,0.00352395,0. 92 0.00373078,0.00371234,0.00381532,0. 93 0.00398425,0.00411183,0.00410042,0. 94 0.00436082,0.00452075,0.00451934,0. 95 0.00490102,0.00510782,0.00511801,0. 96 0.0056677, 0.00594482,0.00596999,0. 97 0.00676236,0.00680914,0.00719407,0. 98 0.00835577,0.0084276, 0.00850759,0. 99 0.0101059, 0.0102249, 0.0109763, 0. 100 0.0125898, 0.0136311, 0.0138081, 0. 101 0.016168, 0.0176664, 0.0179339, 0. 102 0.0214374, 0.0237377, 0.0241275, 0. 103 0.0294621, 0.0300526, 0.0338619, 0. 104 }; 105 106 G4double G4Generator2BN::ctab[320] = 107 { 0.482253, 0.482253, 0.489021, 0.489021, 108 0.495933, 0.495933, 0.495933, 0.502993, 109 0.510204, 0.510204, 0.517572, 0.517572, 110 0.5251, 0.532793, 0.532793, 0.532793, 111 0.548697, 0.548697, 0.556917, 0.556917, 112 0.573921, 0.582717, 0.582717, 0.591716, 113 0.610352, 0.610352, 0.620001, 0.620001, 114 0.64, 0.650364, 0.650364, 0.660982, 115 0.683013, 0.694444, 0.694444, 0.706165, 116 0.730514, 0.743163, 0.743163, 0.756144, 117 0.783147, 0.797194, 0.797194, 0.811622, 118 0.84168, 0.857339, 0.857339, 0.873439, 119 0.907029, 0.924556, 0.924556, 0.942596, 120 0.980296, 1, 1, 1.0203, 121 1.06281, 1.06281, 1.08507, 1.08507, 122 1.13173, 1.1562, 1.1562, 1.18147, 123 1.23457, 1.23457, 1.26247, 1.26247, 124 1.32118, 1.35208, 1.35208, 1.38408, 125 1.45159, 1.45159, 1.45159, 1.48721, 126 1.5625, 1.5625, 1.60231, 1.60231, 127 1.68663, 1.68663, 1.7313, 1.7313, 128 1.82615, 1.87652, 1.87652, 1.92901, 129 1.98373, 2.04082, 2.04082, 2.1004, 130 2.22767, 2.22767, 2.22767, 2.29568, 131 2.44141, 2.44141, 2.51953, 2.51953, 132 2.68745, 2.68745, 2.77778, 2.77778, 133 2.97265, 2.97265, 3.07787, 3.07787, 134 3.30579, 3.30579, 3.42936, 3.42936, 135 3.69822, 3.84468, 3.84468, 3.84468, 136 4.16493, 4.34028, 4.34028, 4.34028, 137 4.7259, 4.93827, 4.93827, 4.93827, 138 5.40833, 5.40833, 5.66893, 5.66893, 139 6.25, 6.25, 6.57462, 6.57462, 140 7.3046, 7.3046, 7.71605, 7.71605, 141 8.65052, 8.65052, 8.65052, 9.18274, 142 9.76562, 10.4058, 10.4058, 10.4058, 143 11.8906, 11.8906, 12.7551, 12.7551, 144 13.7174, 14.7929, 14.7929, 14.7929, 145 17.3611, 17.3611, 17.3611, 18.9036, 146 20.6612, 20.6612, 22.6757, 22.6757, 147 25, 25, 27.7008, 27.7008, 148 30.8642, 30.8642, 34.6021, 34.6021, 149 39.0625, 39.0625, 39.0625, 44.4444, 150 51.0204, 51.0204, 51.0204, 51.0204, 151 59.1716, 59.1716, 69.4444, 69.4444, 152 82.6446, 82.6446, 82.6446, 82.6446, 153 }; 154 155 //....oooOO0OOooo........oooOO0OOooo........oo 156 157 G4Generator2BN::G4Generator2BN(const G4String& 158 : G4VEmAngularDistribution("AngularGen2BN") 159 { 160 b = 1.2; 161 index_min = -300; 162 index_max = 319; 163 164 // Set parameters minimum limits Ekelectron 165 kmin = 100*eV; 166 Ekmin = 250*eV; 167 kcut = 100*eV; 168 169 // Increment Theta value for surface interpo 170 dtheta = 0.1*rad; 171 172 nwarn = 0; 173 174 // Construct Majorant Surface to 2BN cross-s 175 // (we dont need this. Pre-calculated values 176 // ConstructMajorantSurface(); 177 } 178 179 //....oooOO0OOooo........oooOO0OOooo........oo 180 181 G4ThreeVector& G4Generator2BN::SampleDirection 182 G4double out_energy, 183 G4int, 184 const G4Material*) 185 { 186 G4double Ek = dp->GetKineticEnergy(); 187 G4double Eel = dp->GetTotalEnergy(); 188 if(Eel > 2*MeV) { 189 return fGenerator2BS.SampleDirection(dp, o 190 } 191 192 G4double k = Eel - out_energy; 193 194 G4double t; 195 G4double cte2; 196 G4double y, u; 197 G4double ds; 198 G4double dmax; 199 200 // find table index 201 G4int index = G4int(std::log10(Ek/MeV)*100) 202 if(index > index_max) { index = index_max; } 203 else if(index < 0) { index = 0; } 204 205 G4double c = ctab[index]; 206 G4double A = Atab[index]; 207 if(index < index_max) { A = std::max(A, Atab 208 209 do{ 210 // generate k accordimg to std::pow(k,-b) 211 212 // normalization constant 213 // cte1 = (1-b)/(std::pow(kmax,(1-b))-std 214 // y = G4UniformRand(); 215 // k = std::pow(((1-b)/cte1*y+std::pow(km 216 217 // generate theta accordimg to theta/(1+c* 218 // Normalization constant 219 cte2 = 2*c/std::log(1+c*pi2); 220 221 y = G4UniformRand(); 222 t = std::sqrt((G4Exp(2*c*y/cte2)-1)/c); 223 u = G4UniformRand(); 224 225 // point acceptance algorithm 226 //fk = std::pow(k,-b); 227 //ft = t/(1+c*t*t); 228 dmax = A*std::pow(k,-b)*t/(1+c*t*t); 229 ds = Calculatedsdkdt(k,t,Eel); 230 231 // violation point 232 if(ds > dmax && nwarn >= 20) { 233 ++nwarn; 234 G4cout << "### WARNING in G4Generator2BN 235 << " D(Ekin,k)/Dmax-1= " << (ds/dmax - 236 << " results are not reliable!" 237 << G4endl; 238 if(20 == nwarn) { 239 G4cout << " WARNING in G4Generator2BN is c 240 } 241 } 242 243 } while(u*dmax > ds || t > CLHEP::pi); 244 245 G4double sint = std::sin(t); 246 G4double phi = twopi*G4UniformRand(); 247 248 fLocalDirection.set(sint*std::cos(phi), sint 249 fLocalDirection.rotateUz(dp->GetMomentumDire 250 251 return fLocalDirection; 252 } 253 254 //....oooOO0OOooo........oooOO0OOooo........oo 255 256 G4double G4Generator2BN::CalculateFkt(G4double 257 { 258 G4double Fkt_value = 0; 259 Fkt_value = A*std::pow(k,-b)*theta/(1+c*thet 260 return Fkt_value; 261 } 262 263 //....oooOO0OOooo........oooOO0OOooo........oo 264 265 G4double G4Generator2BN::Calculatedsdkdt(G4dou 266 { 267 G4double dsdkdt_value = 0.; 268 G4double Z = 1; 269 // classic radius (in cm) 270 G4double r0 = 2.82E-13; 271 //squared classic radius (in barn) 272 G4double r02 = r0*r0*1.0E+24; 273 274 // Photon energy cannot be greater than elec 275 if(kout > (Eel-electron_mass_c2)){ 276 dsdkdt_value = 0; 277 return dsdkdt_value; 278 } 279 280 G4double E0 = Eel/electron_mass_c2; 281 G4double k = kout/electron_mass_c2; 282 G4double E = E0 - k; 283 284 if(E <= 1*MeV ){ 285 dsdkdt_value = 0; 286 return dsdkdt_value; 287 } 288 289 G4double p0 = std::sqrt(E0*E0-1); 290 G4double p = std::sqrt(E*E-1); 291 G4double LL = std::log((E*E0-1+p*p0)/(E*E0-1 292 G4double delta0 = E0 - p0*std::cos(theta); 293 G4double epsilon = std::log((E+p)/(E-p)); 294 G4double Z2 = Z*Z; 295 G4double sintheta2 = std::sin(theta)*std::si 296 G4double E02 = E0*E0; 297 G4double E2 = E*E; 298 G4double p02 = E0*E0-1; 299 G4double k2 = k*k; 300 G4double delta02 = delta0*delta0; 301 G4double delta04 = delta02* delta02; 302 G4double Q = std::sqrt(p02+k2-2*k*p0*std::co 303 G4double Q2 = Q*Q; 304 G4double epsilonQ = std::log((Q+p)/(Q-p)); 305 306 307 dsdkdt_value = Z2 * (r02/(8*pi*137)) * (1/k) 308 ( (8 * (sintheta2*(2*E02+1))/(p02*delta04) 309 ((2*(5*E02+2*E*E0+3))/(p02 * delta02)) - 310 ((2*(p02-k2))/((Q2*delta02))) + 311 ((4*E)/(p02*delta0)) + 312 (LL/(p*p0))*( 313 ((4*E0*sintheta2*(3*k-p02*E))/(p02*delt 314 ((4*E02*(E02+E2))/(p02*delta02)) + 315 ((2-2*(7*E02-3*E*E0+E2))/(p02*delta02)) 316 (2*k*(E02+E*E0-1))/((p02*delta0)) 317 ) - 318 ((4*epsilon)/(p*delta0)) + 319 ((epsilonQ)/(p*Q))* 320 (4/delta02-(6*k/delta0)-(2*k*(p02-k2))/( 321 ); 322 323 dsdkdt_value = dsdkdt_value*std::sin(theta); 324 return dsdkdt_value; 325 } 326 327 //....oooOO0OOooo........oooOO0OOooo........oo 328 329 void G4Generator2BN::ConstructMajorantSurface( 330 { 331 G4double Eel; 332 G4int vmax; 333 G4double Ek; 334 G4double k, theta, k0, theta0, thetamax; 335 G4double ds, df, dsmax, dk, dt; 336 G4double ratmin; 337 G4double ratio = 0; 338 G4double Vds, Vdf; 339 G4double A, c; 340 341 G4cout << "**** Constructing Majorant Surfac 342 343 if(kcut > kmin) kmin = kcut; 344 345 G4int i = 0; 346 // Loop on energy index 347 for(G4int index = index_min; index < index_m 348 349 G4double fraction = index/100.; 350 Ek = std::pow(10.,fraction); 351 Eel = Ek + electron_mass_c2; 352 353 // find x-section maximum at k=kmin 354 dsmax = 0.; 355 thetamax = 0.; 356 357 for(theta = 0.; theta < pi; theta = theta + 358 359 ds = Calculatedsdkdt(kmin, theta, Eel); 360 361 if(ds > dsmax){ 362 dsmax = ds; 363 thetamax = theta; 364 } 365 } 366 367 // Compute surface paremeters at kmin 368 if(Ek < kmin || thetamax == 0){ 369 c = 0; 370 A = 0; 371 }else{ 372 c = 1/(thetamax*thetamax); 373 A = 2*std::sqrt(c)*dsmax/(std::pow(kmin,-b 374 } 375 376 // look for correction factor to normalizati 377 ratmin = 1.; 378 379 // Volume under surfaces 380 Vds = 0.; 381 Vdf = 0.; 382 k0 = 0.; 383 theta0 = 0.; 384 385 vmax = G4int(100.*std::log10(Ek/kmin)); 386 387 for(G4int v = 0; v < vmax; v++){ 388 G4double fractionLocal = (v/100.); 389 k = std::pow(10.,fractionLocal)*kmin; 390 391 for(theta = 0.; theta < pi; theta = theta 392 dk = k - k0; 393 dt = theta - theta0; 394 ds = Calculatedsdkdt(k,theta, Eel); 395 Vds = Vds + ds*dk*dt; 396 df = CalculateFkt(k,theta, A, c); 397 Vdf = Vdf + df*dk*dt; 398 399 if(ds != 0.){ 400 if(df != 0.) ratio = df/ds; 401 } 402 403 if(ratio < ratmin && ratio != 0.){ 404 ratmin = ratio; 405 } 406 } 407 } 408 409 // sampling efficiency 410 Atab[i] = A/ratmin * 1.04; 411 ctab[i] = c; 412 413 // G4cout << Ek << " " << i << " " << index < 414 i++; 415 } 416 } 417 418 //....oooOO0OOooo........oooOO0OOooo........oo 419 420 void G4Generator2BN::PrintGeneratorInformation 421 { 422 G4cout << "\n" << G4endl; 423 G4cout << "Bremsstrahlung Angular Generator 424 G4cout << "\n" << G4endl; 425 } 426