Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // Author: Sebastien Incerti 26 // Author: Sebastien Incerti 27 // 22 January 2012 27 // 22 January 2012 28 // on base of G4LivermoreGammaConversi << 28 // on base of G4LivermoreGammaConversionModel 29 // and G4LivermoreRayleighModel (MT ve << 30 // << 31 // Modifications: Zhuxin Li@CENBG << 32 // 11 March 2020 << 33 // derives from G4PairProductio << 34 // ------------------------------------------- << 35 29 36 #include "G4LivermoreGammaConversionModel.hh" 30 #include "G4LivermoreGammaConversionModel.hh" 37 << 38 #include "G4AutoLock.hh" << 39 #include "G4Electron.hh" << 40 #include "G4EmParameters.hh" << 41 #include "G4Exp.hh" << 42 #include "G4ParticleChangeForGamma.hh" << 43 #include "G4PhysicalConstants.hh" 31 #include "G4PhysicalConstants.hh" 44 #include "G4PhysicsFreeVector.hh" << 45 #include "G4SystemOfUnits.hh" 32 #include "G4SystemOfUnits.hh" 46 33 47 namespace << 48 { << 49 G4Mutex LivermoreGammaConversionModelMutex = G << 50 } << 51 << 52 //....oooOO0OOooo........oooOO0OOooo........oo << 53 //....oooOO0OOooo........oooOO0OOooo........oo 34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 54 35 55 G4PhysicsFreeVector* G4LivermoreGammaConversio << 36 using namespace std; 56 G4String G4LivermoreGammaConversionModel::gDat << 37 >> 38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 57 39 58 G4LivermoreGammaConversionModel::G4LivermoreGa << 40 G4LivermoreGammaConversionModel::G4LivermoreGammaConversionModel(const G4ParticleDefinition*, 59 << 41 const G4String& nam) 60 : G4PairProductionRelModel(p, nam) << 42 :G4VEmModel(nam),smallEnergy(2.*MeV),isInitialised(false),maxZ(99) 61 { 43 { 62 fParticleChange = nullptr; << 44 fParticleChange = 0; 63 lowEnergyLimit = 2. * CLHEP::electron_mass_c << 45 64 verboseLevel = 0; << 46 lowEnergyLimit = 2.0*electron_mass_c2; >> 47 data.resize(maxZ+1,0); >> 48 >> 49 verboseLevel= 0; 65 // Verbosity scale for debugging purposes: 50 // Verbosity scale for debugging purposes: 66 // 0 = nothing << 51 // 0 = nothing 67 // 1 = calculation of cross sections, file o 52 // 1 = calculation of cross sections, file openings... 68 // 2 = entering in methods 53 // 2 = entering in methods 69 if (verboseLevel > 0) { << 54 >> 55 if(verboseLevel > 0) >> 56 { 70 G4cout << "G4LivermoreGammaConversionModel 57 G4cout << "G4LivermoreGammaConversionModel is constructed " << G4endl; 71 } 58 } 72 } 59 } 73 60 74 //....oooOO0OOooo........oooOO0OOooo........oo 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 75 62 76 G4LivermoreGammaConversionModel::~G4LivermoreG 63 G4LivermoreGammaConversionModel::~G4LivermoreGammaConversionModel() 77 { << 64 { 78 if (IsMaster()) { << 65 for(G4int i=0; i<=maxZ; ++i) { delete data[i]; } 79 for (G4int i = 0; i <= maxZ; ++i) { << 80 if (data[i]) { << 81 delete data[i]; << 82 data[i] = nullptr; << 83 } << 84 } << 85 } << 86 } 66 } 87 67 88 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 89 69 90 void G4LivermoreGammaConversionModel::Initiali << 70 void 91 << 71 G4LivermoreGammaConversionModel::Initialise(const G4ParticleDefinition* particle, >> 72 const G4DataVector& cuts) 92 { 73 { 93 G4PairProductionRelModel::Initialise(particl << 74 if (verboseLevel > 1) 94 if (verboseLevel > 1) { << 75 { 95 G4cout << "Calling Initialise() of G4Liver 76 G4cout << "Calling Initialise() of G4LivermoreGammaConversionModel." << G4endl 96 << "Energy range: " << LowEnergyLim << 77 << "Energy range: " 97 << " GeV isMater: " << IsMaster() < << 78 << LowEnergyLimit() / MeV << " MeV - " 98 } << 79 << HighEnergyLimit() / GeV << " GeV" 99 << 80 << G4endl; 100 if (IsMaster()) { << 81 } 101 // Initialise element selector << 82 102 InitialiseElementSelectors(particle, cuts) << 83 // Initialise element selector 103 << 84 104 // Access to elements << 85 InitialiseElementSelectors(particle, cuts); 105 const G4ElementTable* elemTable = G4Elemen << 86 106 std::size_t numElems = (*elemTable).size() << 87 // Access to elements 107 for (std::size_t ie = 0; ie < numElems; ++ << 88 108 const G4Element* elem = (*elemTable)[ie] << 89 char* path = getenv("G4LEDATA"); 109 const G4int Z = std::min(maxZ, elem->Get << 90 110 if (data[Z] == nullptr) { << 91 G4ProductionCutsTable* theCoupleTable = 111 ReadData(Z); << 92 G4ProductionCutsTable::GetProductionCutsTable(); 112 } << 93 G4int numOfCouples = theCoupleTable->GetTableSize(); >> 94 >> 95 for(G4int i=0; i<numOfCouples; ++i) >> 96 { >> 97 const G4Material* material = >> 98 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial(); >> 99 const G4ElementVector* theElementVector = material->GetElementVector(); >> 100 G4int nelm = material->GetNumberOfElements(); >> 101 >> 102 for (G4int j=0; j<nelm; ++j) >> 103 { >> 104 >> 105 G4int Z = (G4int)(*theElementVector)[j]->GetZ(); >> 106 if(Z < 1) { Z = 1; } >> 107 else if(Z > maxZ) { Z = maxZ; } >> 108 if(!data[Z]) { ReadData(Z, path); } 113 } 109 } 114 } 110 } 115 if (isInitialised) { << 111 // 116 return; << 112 117 } << 113 if(isInitialised) { return; } 118 fParticleChange = GetParticleChangeForGamma( 114 fParticleChange = GetParticleChangeForGamma(); 119 isInitialised = true; 115 isInitialised = true; 120 } 116 } 121 117 122 //....oooOO0OOooo........oooOO0OOooo........oo 118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 123 119 124 const G4String& G4LivermoreGammaConversionMode << 120 void G4LivermoreGammaConversionModel::ReadData(size_t Z, const char* path) 125 { 121 { 126 // no check in this method - environment var << 122 if (verboseLevel > 1) 127 if (gDataDirectory.empty()) { << 123 { 128 auto param = G4EmParameters::Instance(); << 124 G4cout << "Calling ReadData() of G4LivermoreGammaConversionModel" 129 std::ostringstream ost; << 125 << G4endl; 130 if (param->LivermoreDataDir() == "livermor << 126 } 131 ost << param->GetDirLEDATA() << "/liverm << 127 132 useSpline = true; << 128 if(data[Z]) { return; } >> 129 >> 130 const char* datadir = path; >> 131 >> 132 if(!datadir) >> 133 { >> 134 datadir = getenv("G4LEDATA"); >> 135 if(!datadir) >> 136 { >> 137 G4Exception("G4LivermoreGammaConversionModel::ReadData()", >> 138 "em0006",FatalException, >> 139 "Environment variable G4LEDATA not defined"); >> 140 return; 133 } 141 } 134 else { << 135 ost << param->GetDirLEDATA() << "/epics2 << 136 } << 137 gDataDirectory = ost.str(); << 138 } << 139 return gDataDirectory; << 140 } << 141 << 142 //....oooOO0OOooo........oooOO0OOooo........oo << 143 << 144 void G4LivermoreGammaConversionModel::ReadData << 145 { << 146 if (verboseLevel > 1) { << 147 G4cout << "Calling ReadData() of G4Livermo << 148 } << 149 << 150 if (data[Z] != nullptr) { << 151 return; << 152 } 142 } 153 143 >> 144 // >> 145 >> 146 data[Z] = new G4LPhysicsFreeVector(); >> 147 >> 148 // Activation of spline interpolation >> 149 data[Z] ->SetSpline(true); >> 150 // >> 151 154 std::ostringstream ost; 152 std::ostringstream ost; 155 ost << FindDirectoryPath() << "pp-cs-" << Z << 153 ost << datadir << "/livermore/pair/pp-cs-" << Z <<".dat"; 156 << 157 data[Z] = new G4PhysicsFreeVector(useSpline) << 158 << 159 std::ifstream fin(ost.str().c_str()); 154 std::ifstream fin(ost.str().c_str()); 160 << 155 161 if (!fin.is_open()) { << 156 if( !fin.is_open()) >> 157 { 162 G4ExceptionDescription ed; 158 G4ExceptionDescription ed; 163 ed << "G4LivermoreGammaConversionModel dat << 159 ed << "G4LivermoreGammaConversionModel data file <" << ost.str().c_str() 164 << G4endl; << 160 << "> is not opened!" << G4endl; 165 G4Exception("G4LivermoreGammaConversionMod << 161 G4Exception("G4LivermoreGammaConversionModel::ReadData()", 166 "G4LEDATA version should be G4 << 162 "em0003",FatalException, >> 163 ed,"G4LEDATA version should be G4EMLOW6.27 or later."); 167 return; 164 return; 168 } << 165 } 169 else { << 166 170 if (verboseLevel > 1) { << 167 else 171 G4cout << "File " << ost.str() << " is o << 168 { 172 } << 169 173 << 170 if(verboseLevel > 3) { G4cout << "File " << ost.str() >> 171 << " is opened by G4LivermoreGammaConversionModel" << G4endl;} >> 172 174 data[Z]->Retrieve(fin, true); 173 data[Z]->Retrieve(fin, true); 175 } << 174 } 176 // Activation of spline interpolation << 175 177 if (useSpline) data[Z]->FillSecondDerivative << 176 178 } 177 } 179 178 180 //....oooOO0OOooo........oooOO0OOooo........oo 179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 181 180 182 G4double << 181 G4double 183 G4LivermoreGammaConversionModel::ComputeCrossS << 182 G4LivermoreGammaConversionModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*, 184 << 183 G4double GammaEnergy, 185 << 184 G4double Z, G4double, >> 185 G4double, G4double) 186 { 186 { 187 if (verboseLevel > 1) { << 187 if (verboseLevel > 1) 188 G4cout << "G4LivermoreGammaConversionModel << 188 { >> 189 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreGammaConversionModel" >> 190 << G4endl; 189 } 191 } 190 192 191 if (GammaEnergy < lowEnergyLimit) { << 193 if (GammaEnergy < lowEnergyLimit) { return 0.0; } 192 return 0.0; << 193 } << 194 194 195 G4double xs = 0.0; 195 G4double xs = 0.0; >> 196 >> 197 G4int intZ=G4int(Z); 196 198 197 G4int intZ = std::max(1, std::min(G4lrint(Z) << 199 if(intZ < 1 || intZ > maxZ) { return xs; } 198 200 199 G4PhysicsFreeVector* pv = data[intZ]; << 201 G4LPhysicsFreeVector* pv = data[intZ]; 200 202 201 // if element was not initialised << 203 // element was not initialised 202 // do initialisation safely for MT mode << 204 if(!pv) 203 if (pv == nullptr) { << 205 { 204 InitialiseForElement(particle, intZ); << 206 char* path = getenv("G4LEDATA"); >> 207 ReadData(intZ, path); 205 pv = data[intZ]; 208 pv = data[intZ]; 206 if (pv == nullptr) { << 209 if(!pv) { return xs; } 207 return xs; << 208 } << 209 } 210 } 210 // x-section is taken from the table 211 // x-section is taken from the table 211 xs = pv->Value(GammaEnergy); << 212 xs = pv->Value(GammaEnergy); 212 213 213 if (verboseLevel > 0) { << 214 if(verboseLevel > 0) 214 G4cout << "*** Gamma conversion xs for Z=" << 215 { 215 << " cs=" << xs / millibarn << " m << 216 G4int n = pv->GetVectorLength() - 1; >> 217 G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)=" << GammaEnergy/MeV << G4endl; >> 218 G4cout << " cs (Geant4 internal unit)=" << xs << G4endl; >> 219 G4cout << " -> first cs value in EADL data file (iu) =" << (*pv)[0] << G4endl; >> 220 G4cout << " -> last cs value in EADL data file (iu) =" << (*pv)[n] << G4endl; >> 221 G4cout << "*********************************************************" << G4endl; 216 } 222 } >> 223 217 return xs; 224 return xs; >> 225 218 } 226 } 219 227 220 //....oooOO0OOooo........oooOO0OOooo........oo << 228 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 221 229 222 void G4LivermoreGammaConversionModel::Initiali << 230 void G4LivermoreGammaConversionModel::SampleSecondaries( >> 231 std::vector<G4DynamicParticle*>* fvect, >> 232 const G4MaterialCutsCouple* couple, >> 233 const G4DynamicParticle* aDynamicGamma, >> 234 G4double, G4double) 223 { 235 { 224 G4AutoLock l(&LivermoreGammaConversionModelM << 236 225 if (data[Z] == nullptr) { << 237 // The energies of the e+ e- secondaries are sampled using the Bethe - Heitler 226 ReadData(Z); << 238 // cross sections with Coulomb correction. A modified version of the random 227 } << 239 // number techniques of Butcher & Messel is used (Nuc Phys 20(1960),15). 228 l.unlock(); << 240 >> 241 // Note 1 : Effects due to the breakdown of the Born approximation at low >> 242 // energy are ignored. >> 243 // Note 2 : The differential cross section implicitly takes account of >> 244 // pair creation in both nuclear and atomic electron fields. However triplet >> 245 // prodution is not generated. >> 246 >> 247 if (verboseLevel > 1) >> 248 G4cout << "Calling SampleSecondaries() of G4LivermoreGammaConversionModel" << G4endl; >> 249 >> 250 G4double photonEnergy = aDynamicGamma->GetKineticEnergy(); >> 251 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection(); >> 252 >> 253 G4double epsilon ; >> 254 G4double epsilon0Local = electron_mass_c2 / photonEnergy ; >> 255 >> 256 // Do it fast if photon energy < 2. MeV >> 257 if (photonEnergy < smallEnergy ) >> 258 { >> 259 epsilon = epsilon0Local + (0.5 - epsilon0Local) * G4UniformRand(); >> 260 } >> 261 else >> 262 { >> 263 // Select randomly one element in the current material >> 264 >> 265 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition(); >> 266 const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy); >> 267 >> 268 if (element == 0) >> 269 { >> 270 G4cout << "G4LivermoreGammaConversionModel::SampleSecondaries - element = 0" >> 271 << G4endl; >> 272 return; >> 273 } >> 274 G4IonisParamElm* ionisation = element->GetIonisation(); >> 275 if (ionisation == 0) >> 276 { >> 277 G4cout << "G4LivermoreGammaConversionModel::SampleSecondaries - ionisation = 0" >> 278 << G4endl; >> 279 return; >> 280 } >> 281 >> 282 // Extract Coulomb factor for this Element >> 283 G4double fZ = 8. * (ionisation->GetlogZ3()); >> 284 if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb()); >> 285 >> 286 // Limits of the screening variable >> 287 G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3()) ; >> 288 G4double screenMax = std::exp ((42.24 - fZ)/8.368) - 0.952 ; >> 289 G4double screenMin = std::min(4.*screenFactor,screenMax) ; >> 290 >> 291 // Limits of the energy sampling >> 292 G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ; >> 293 G4double epsilonMin = std::max(epsilon0Local,epsilon1); >> 294 G4double epsilonRange = 0.5 - epsilonMin ; >> 295 >> 296 // Sample the energy rate of the created electron (or positron) >> 297 G4double screen; >> 298 G4double gReject ; >> 299 >> 300 G4double f10 = ScreenFunction1(screenMin) - fZ; >> 301 G4double f20 = ScreenFunction2(screenMin) - fZ; >> 302 G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.); >> 303 G4double normF2 = std::max(1.5 * f20,0.); >> 304 >> 305 do >> 306 { >> 307 if (normF1 / (normF1 + normF2) > G4UniformRand() ) >> 308 { >> 309 epsilon = 0.5 - epsilonRange * std::pow(G4UniformRand(), 0.333333) ; >> 310 screen = screenFactor / (epsilon * (1. - epsilon)); >> 311 gReject = (ScreenFunction1(screen) - fZ) / f10 ; >> 312 } >> 313 else >> 314 { >> 315 epsilon = epsilonMin + epsilonRange * G4UniformRand(); >> 316 screen = screenFactor / (epsilon * (1 - epsilon)); >> 317 gReject = (ScreenFunction2(screen) - fZ) / f20 ; >> 318 } >> 319 } while ( gReject < G4UniformRand() ); >> 320 >> 321 } // End of epsilon sampling >> 322 >> 323 // Fix charges randomly >> 324 >> 325 G4double electronTotEnergy; >> 326 G4double positronTotEnergy; >> 327 >> 328 if (G4UniformRand() > 0.5) >> 329 { >> 330 electronTotEnergy = (1. - epsilon) * photonEnergy; >> 331 positronTotEnergy = epsilon * photonEnergy; >> 332 } >> 333 else >> 334 { >> 335 positronTotEnergy = (1. - epsilon) * photonEnergy; >> 336 electronTotEnergy = epsilon * photonEnergy; >> 337 } >> 338 >> 339 // Scattered electron (positron) angles. ( Z - axis along the parent photon) >> 340 // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211), >> 341 // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977) >> 342 >> 343 G4double u; >> 344 const G4double a1 = 0.625; >> 345 G4double a2 = 3. * a1; >> 346 // G4double d = 27. ; >> 347 >> 348 // if (9. / (9. + d) > G4UniformRand()) >> 349 if (0.25 > G4UniformRand()) >> 350 { >> 351 u = - std::log(G4UniformRand() * G4UniformRand()) / a1 ; >> 352 } >> 353 else >> 354 { >> 355 u = - std::log(G4UniformRand() * G4UniformRand()) / a2 ; >> 356 } >> 357 >> 358 G4double thetaEle = u*electron_mass_c2/electronTotEnergy; >> 359 G4double thetaPos = u*electron_mass_c2/positronTotEnergy; >> 360 G4double phi = twopi * G4UniformRand(); >> 361 >> 362 G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle); >> 363 G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos); >> 364 >> 365 >> 366 // Kinematics of the created pair: >> 367 // the electron and positron are assumed to have a symetric angular >> 368 // distribution with respect to the Z axis along the parent photon >> 369 >> 370 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ; >> 371 >> 372 G4ThreeVector electronDirection (dxEle, dyEle, dzEle); >> 373 electronDirection.rotateUz(photonDirection); >> 374 >> 375 G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(), >> 376 electronDirection, >> 377 electronKineEnergy); >> 378 >> 379 // The e+ is always created >> 380 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ; >> 381 >> 382 G4ThreeVector positronDirection (dxPos, dyPos, dzPos); >> 383 positronDirection.rotateUz(photonDirection); >> 384 >> 385 // Create G4DynamicParticle object for the particle2 >> 386 G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(), >> 387 positronDirection, >> 388 positronKineEnergy); >> 389 // Fill output vector >> 390 fvect->push_back(particle1); >> 391 fvect->push_back(particle2); >> 392 >> 393 // kill incident photon >> 394 fParticleChange->SetProposedKineticEnergy(0.); >> 395 fParticleChange->ProposeTrackStatus(fStopAndKill); >> 396 229 } 397 } 230 398 231 //....oooOO0OOooo........oooOO0OOooo........oo << 399 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 400 >> 401 G4double >> 402 G4LivermoreGammaConversionModel::ScreenFunction1(G4double screenVariable) >> 403 { >> 404 // Compute the value of the screening function 3*phi1 - phi2 >> 405 >> 406 G4double value; >> 407 >> 408 if (screenVariable > 1.) >> 409 value = 42.24 - 8.368 * std::log(screenVariable + 0.952); >> 410 else >> 411 value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable); >> 412 >> 413 return value; >> 414 } >> 415 >> 416 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 417 >> 418 G4double >> 419 G4LivermoreGammaConversionModel::ScreenFunction2(G4double screenVariable) >> 420 { >> 421 // Compute the value of the screening function 1.5*phi1 - 0.5*phi2 >> 422 >> 423 G4double value; >> 424 >> 425 if (screenVariable > 1.) >> 426 value = 42.24 - 8.368 * std::log(screenVariable + 0.952); >> 427 else >> 428 value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable); >> 429 >> 430 return value; >> 431 } >> 432 232 433