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 // Author: Sebastien Incerti 27 // 22 January 2012 28 // on base of G4LivermoreNuclearGammaC 29 // and G4LivermoreRayleighModel (MT ve 30 31 #include "G4LivermoreNuclearGammaConversionMod 32 #include "G4PhysicalConstants.hh" 33 #include "G4SystemOfUnits.hh" 34 #include "G4Log.hh" 35 #include "G4Exp.hh" 36 #include "G4AutoLock.hh" 37 38 //....oooOO0OOooo........oooOO0OOooo........oo 39 40 using namespace std; 41 namespace { G4Mutex LivermoreNuclearGammaConve 42 43 //....oooOO0OOooo........oooOO0OOooo........oo 44 45 G4PhysicsFreeVector* G4LivermoreNuclearGammaCo 46 47 G4LivermoreNuclearGammaConversionModel::G4Live 48 (const G4ParticleDefinition*, const G4String& 49 :G4VEmModel(nam),smallEnergy(2.*MeV), 50 isInitialised(false) 51 { 52 fParticleChange = nullptr; 53 54 lowEnergyLimit = 2.0*electron_mass_c2; 55 56 verboseLevel= 0; 57 // Verbosity scale for debugging purposes: 58 // 0 = nothing 59 // 1 = calculation of cross sections, file o 60 // 2 = entering in methods 61 62 if(verboseLevel > 0) 63 { 64 G4cout << "G4LivermoreNuclearGammaConversi 65 } 66 } 67 68 //....oooOO0OOooo........oooOO0OOooo........oo 69 70 G4LivermoreNuclearGammaConversionModel::~G4Liv 71 { 72 if(IsMaster()) { 73 for(G4int i=0; i<maxZ; ++i) { 74 if(data[i]) { 75 delete data[i]; 76 data[i] = 0; 77 } 78 } 79 } 80 } 81 82 //....oooOO0OOooo........oooOO0OOooo........oo 83 84 void G4LivermoreNuclearGammaConversionModel::I 85 const G4Partic 86 const G4DataVector& cuts) 87 { 88 if (verboseLevel > 1) 89 { 90 G4cout << "Calling Initialise() of G4Liver 91 << G4endl 92 << "Energy range: " 93 << LowEnergyLimit() / MeV << " MeV - " 94 << HighEnergyLimit() / GeV << " GeV" 95 << G4endl; 96 } 97 98 if(IsMaster()) 99 { 100 101 // Initialise element selector 102 InitialiseElementSelectors(particle, cuts) 103 104 // Access to elements 105 const char* path = G4FindDataDir("G4LEDATA 106 107 G4ProductionCutsTable* theCoupleTable = 108 G4ProductionCutsTable::GetProductionCuts 109 110 G4int numOfCouples = (G4int)theCoupleTable 111 112 for(G4int i=0; i<numOfCouples; ++i) 113 { 114 const G4Material* material = 115 theCoupleTable->GetMaterialCutsCouple( 116 const G4ElementVector* theElementVector 117 std::size_t nelm = material->GetNumberOf 118 119 for (std::size_t j=0; j<nelm; ++j) 120 { 121 G4int Z = (G4int)(*theElementVector)[j 122 if(Z < 1) { Z = 1; } 123 else if(Z > maxZ) { Z = maxZ; } 124 if(!data[Z]) { ReadData(Z, path); } 125 } 126 } 127 } 128 if(isInitialised) { return; } 129 fParticleChange = GetParticleChangeForGamma( 130 isInitialised = true; 131 } 132 133 //....oooOO0OOooo........oooOO0OOooo........oo 134 135 void G4LivermoreNuclearGammaConversionModel::I 136 const G4ParticleDefinition*, G4VEmModel* 137 { 138 SetElementSelectors(masterModel->GetElementS 139 } 140 141 //....oooOO0OOooo........oooOO0OOooo........oo 142 143 G4double 144 G4LivermoreNuclearGammaConversionModel::MinPri 145 const G4ParticleDefinition*, 146 G4double) 147 { 148 return lowEnergyLimit; 149 } 150 151 //....oooOO0OOooo........oooOO0OOooo........oo 152 153 void G4LivermoreNuclearGammaConversionModel::R 154 { 155 if (verboseLevel > 1) 156 { 157 G4cout << "Calling ReadData() of G4Livermo 158 << G4endl; 159 } 160 161 162 if(data[Z]) { return; } 163 164 const char* datadir = path; 165 166 if(!datadir) 167 { 168 datadir = G4FindDataDir("G4LEDATA"); 169 if(!datadir) 170 { 171 G4Exception("G4LivermoreNuclearGammaConv 172 "em0006",FatalException, 173 "Environment variable G4LEDATA not defin 174 return; 175 } 176 } 177 178 data[Z] = new G4PhysicsFreeVector(0,/*spline 179 180 std::ostringstream ost; 181 ost << datadir << "/livermore/pairdata/pp-pa 182 std::ifstream fin(ost.str().c_str()); 183 184 if( !fin.is_open()) 185 { 186 G4ExceptionDescription ed; 187 ed << "G4LivermoreNuclearGammaConversionMo 188 << "> is not opened!" << G4endl; 189 G4Exception("G4LivermoreNuclearGammaConver 190 "em0003",FatalException, 191 ed,"G4LEDATA version should be G4EMLOW8.0 192 return; 193 } 194 else 195 { 196 197 if(verboseLevel > 3) { G4cout << "File " 198 << " is opened by G4LivermoreNucle 199 200 data[Z]->Retrieve(fin, true); 201 } 202 203 // Activation of spline interpolation 204 data[Z] ->FillSecondDerivatives(); 205 } 206 207 //....oooOO0OOooo........oooOO0OOooo........oo 208 209 G4double 210 G4LivermoreNuclearGammaConversionModel::Comput 211 G4double GammaEnergy, 212 G4double Z, G4double, 213 G4double, G4double) 214 { 215 if (verboseLevel > 1) 216 { 217 G4cout << "Calling ComputeCrossSectionPerA 218 << G4endl; 219 } 220 221 if (GammaEnergy < lowEnergyLimit) { return 0 222 223 G4double xs = 0.0; 224 225 G4int intZ=G4int(Z); 226 227 if(intZ < 1 || intZ > maxZ) { return xs; } 228 229 G4PhysicsFreeVector* pv = data[intZ]; 230 231 // if element was not initialised 232 // do initialisation safely for MT mode 233 if(!pv) 234 { 235 InitialiseForElement(0, intZ); 236 pv = data[intZ]; 237 if(!pv) { return xs; } 238 } 239 // x-section is taken from the table 240 xs = pv->Value(GammaEnergy); 241 242 if(verboseLevel > 0) 243 { 244 std::size_t n = pv->GetVectorLength() - 1; 245 G4cout << "****** DEBUG: tcs value for Z 246 << GammaEnergy/MeV << G4endl; 247 G4cout << " cs (Geant4 internal unit)=" 248 G4cout << " -> first cs value in EADL 249 G4cout << " -> last cs value in EADL 250 G4cout << "***************************** 251 } 252 return xs; 253 } 254 255 //....oooOO0OOooo........oooOO0OOooo........oo 256 257 void G4LivermoreNuclearGammaConversionModel::S 258 std::vector<G 259 const G4MaterialCutsCouple* couple, 260 const G4DynamicParticle* aDynamicGamm 261 G4double, G4double) 262 { 263 // The energies of the e+ e- secondaries are 264 // cross sections with Coulomb correction. A 265 // number techniques of Butcher & Messel is 266 267 // Note 1 : Effects due to the breakdown of 268 // energy are ignored. 269 // Note 2 : The differential cross section i 270 // pair creation in both nuclear and atomic 271 // prodution is not generated. 272 273 if (verboseLevel > 1) { 274 G4cout << "Calling SampleSecondaries() of 275 << G4endl; 276 } 277 278 G4double photonEnergy = aDynamicGamma->GetKi 279 G4ParticleMomentum photonDirection = aDynami 280 281 G4double epsilon ; 282 G4double epsilon0Local = electron_mass_c2 / 283 284 // Do it fast if photon energy < 2. MeV 285 if (photonEnergy < smallEnergy ) 286 { 287 epsilon = epsilon0Local + (0.5 - epsilon0L 288 } 289 else 290 { 291 // Select randomly one element in the curr 292 const G4ParticleDefinition* particle = aD 293 const G4Element* element = SelectRandomAto 294 295 if (element == nullptr) 296 { 297 G4cout << "G4LivermoreNuclearGammaConversion 298 << G4endl; 299 return; 300 } 301 G4IonisParamElm* ionisation = element->Get 302 if (ionisation == nullptr) 303 { 304 G4cout << "G4LivermoreNuclearGammaConversion 305 << G4endl; 306 return; 307 } 308 309 // Extract Coulomb factor for this Element 310 G4double fZ = 8. * (ionisation->GetlogZ3() 311 if (photonEnergy > 50. * MeV) fZ += 8. * ( 312 313 // Limits of the screening variable 314 G4double screenFactor = 136. * epsilon0Loc 315 G4double screenMax = G4Exp ((42.24 - fZ)/8 316 G4double screenMin = std::min(4.*screenFac 317 318 // Limits of the energy sampling 319 G4double epsilon1 = 0.5 - 0.5 * std::sqrt( 320 G4double epsilonMin = std::max(epsilon0Loc 321 G4double epsilonRange = 0.5 - epsilonMin ; 322 323 // Sample the energy rate of the created e 324 G4double screen; 325 G4double gReject ; 326 327 G4double f10 = ScreenFunction1(screenMin) 328 G4double f20 = ScreenFunction2(screenMin) 329 G4double normF1 = std::max(f10 * epsilonRa 330 G4double normF2 = std::max(1.5 * f20,0.); 331 332 do 333 { 334 if (normF1 / (normF1 + normF2) > G4UniformRa 335 { 336 epsilon = 0.5 - epsilonRange * std::pow( 337 screen = screenFactor / (epsilon * (1. - 338 gReject = (ScreenFunction1(screen) - fZ) 339 } 340 else 341 { 342 epsilon = epsilonMin + epsilonRange * G4 343 screen = screenFactor / (epsilon * (1 - 344 gReject = (ScreenFunction2(screen) - fZ) 345 } 346 } while ( gReject < G4UniformRand() ); 347 } // End of epsilon sampling 348 349 // Fix charges randomly 350 G4double electronTotEnergy; 351 G4double positronTotEnergy; 352 353 if (G4UniformRand() > 0.5) 354 { 355 electronTotEnergy = (1. - epsilon) * pho 356 positronTotEnergy = epsilon * photonEner 357 } 358 else 359 { 360 positronTotEnergy = (1. - epsilon) * pho 361 electronTotEnergy = epsilon * photonEner 362 } 363 364 // Scattered electron (positron) angles. ( Z 365 // Universal distribution suggested by L. Ur 366 // derived from Tsai distribution (Rev. Mod. 367 368 G4double u; 369 const G4double a1 = 0.625; 370 G4double a2 = 3. * a1; 371 372 if (0.25 > G4UniformRand()) 373 { 374 u = - G4Log(G4UniformRand() * G4UniformR 375 } 376 else 377 { 378 u = - G4Log(G4UniformRand() * G4UniformR 379 } 380 381 G4double thetaEle = u*electron_mass_c2/elect 382 G4double thetaPos = u*electron_mass_c2/posit 383 G4double phi = twopi * G4UniformRand(); 384 385 G4double dxEle= std::sin(thetaEle)*std::cos( 386 G4double dxPos=-std::sin(thetaPos)*std::cos( 387 388 // Kinematics of the created pair: 389 // the electron and positron are assumed to 390 // distribution with respect to the Z axis a 391 392 G4double electronKineEnergy = std::max(0.,el 393 394 G4ThreeVector electronDirection (dxEle, dyEl 395 electronDirection.rotateUz(photonDirection); 396 397 G4DynamicParticle* particle1 = new G4Dynamic 398 electronDirection, 399 electronKineEnergy); 400 401 // The e+ is always created 402 G4double positronKineEnergy = std::max(0.,po 403 404 G4ThreeVector positronDirection (dxPos, dyPo 405 positronDirection.rotateUz(photonDirection); 406 407 // Create G4DynamicParticle object for the p 408 G4DynamicParticle* particle2 = new G4Dynamic 409 positronDirection, 410 positronKineEnergy); 411 // Fill output vector 412 fvect->push_back(particle1); 413 fvect->push_back(particle2); 414 415 // kill incident photon 416 fParticleChange->SetProposedKineticEnergy(0. 417 fParticleChange->ProposeTrackStatus(fStopAnd 418 419 } 420 421 //....oooOO0OOooo........oooOO0OOooo........oo 422 423 G4double 424 G4LivermoreNuclearGammaConversionModel::Screen 425 { 426 // Compute the value of the screening functi 427 428 G4double value; 429 430 if (screenVariable > 1.) 431 value = 42.24 - 8.368 * G4Log(screenVariab 432 else 433 value = 42.392 - screenVariable * (7.796 - 434 435 return value; 436 } 437 438 //....oooOO0OOooo........oooOO0OOooo........oo 439 440 G4double 441 G4LivermoreNuclearGammaConversionModel::Screen 442 { 443 // Compute the value of the screening functi 444 G4double value; 445 446 if (screenVariable > 1.) 447 value = 42.24 - 8.368 * G4Log(screenVariab 448 else 449 value = 41.405 - screenVariable * (5.828 - 450 451 return value; 452 } 453 454 //....oooOO0OOooo........oooOO0OOooo........oo 455 456 void G4LivermoreNuclearGammaConversionModel::I 457 const G4ParticleDefinition 458 G4int Z) 459 { 460 G4AutoLock l(&LivermoreNuclearGammaConversio 461 if(!data[Z]) { ReadData(Z); } 462 l.unlock(); 463 } 464 465 //....oooOO0OOooo........oooOO0OOooo........oo 466