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 // ******************************************* 26 // ********************************************************************* 27 // | 27 // | | 28 // | G4LowEPComptonModel-- Geant4 28 // | G4LowEPComptonModel-- Geant4 Monash University | 29 // | low energy Compton scat 29 // | low energy Compton scattering model. | 30 // | J. M. C. Brown, Monash Univer 30 // | J. M. C. Brown, Monash University, Australia | 31 // | ## Unpolarised photons 31 // | ## Unpolarised photons only ## | 32 // | 32 // | | 33 // | 33 // | | 34 // ******************************************* 34 // ********************************************************************* 35 // | 35 // | | 36 // | The following is a Geant4 class to simula 36 // | The following is a Geant4 class to simulate the process of | 37 // | bound electron Compton scattering. Genera 37 // | bound electron Compton scattering. General code structure is | 38 // | based on G4LowEnergyCompton.cc and G4Live 38 // | based on G4LowEnergyCompton.cc and G4LivermoreComptonModel.cc. | 39 // | Algorithms for photon energy, and ejected 39 // | Algorithms for photon energy, and ejected Compton electron | 40 // | direction taken from: 40 // | direction taken from: | 41 // | 41 // | | 42 // | J. M. C. Brown, M. R. Dimmock, J. E. Gill 42 // | J. M. C. Brown, M. R. Dimmock, J. E. Gillam and D. M. Paganin, | 43 // | "A low energy bound atomic electron Compt 43 // | "A low energy bound atomic electron Compton scattering model | 44 // | for Geant4", NIMB, Vol. 338, 77-88, 2014 44 // | for Geant4", NIMB, Vol. 338, 77-88, 2014. | 45 // | 45 // | | 46 // | The author acknowledges the work of the G 46 // | The author acknowledges the work of the Geant4 collaboration | 47 // | in developing the following algorithms th 47 // | in developing the following algorithms that have been employed | 48 // | or adapeted for the present software: 48 // | or adapeted for the present software: | 49 // | 49 // | | 50 // | # sampling of photon scattering angle, 50 // | # sampling of photon scattering angle, | 51 // | # target element selection in composite 51 // | # target element selection in composite materials, | 52 // | # target shell selection in element, 52 // | # target shell selection in element, | 53 // | # and sampling of bound electron momentu 53 // | # and sampling of bound electron momentum from Compton profiles. | 54 // | 54 // | | 55 // ******************************************* 55 // ********************************************************************* 56 // | 56 // | | 57 // | History: 57 // | History: | 58 // | -------- 58 // | -------- | 59 // | 59 // | | 60 // | Nov. 2011 JMCB - First version 60 // | Nov. 2011 JMCB - First version | 61 // | Feb. 2012 JMCB - Migration to Geant 61 // | Feb. 2012 JMCB - Migration to Geant4 9.5 | 62 // | Sep. 2012 JMCB - Final fixes for Ge 62 // | Sep. 2012 JMCB - Final fixes for Geant4 9.6 | 63 // | Feb. 2013 JMCB - Geant4 9.6 FPE fix 63 // | Feb. 2013 JMCB - Geant4 9.6 FPE fix for bug 1426 | 64 // | Jan. 2015 JMCB - Migration to MT (B 64 // | Jan. 2015 JMCB - Migration to MT (Based on Livermore | 65 // | implementation) 65 // | implementation) | 66 // | Feb. 2016 JMCB - Geant4 10.2 FPE fi 66 // | Feb. 2016 JMCB - Geant4 10.2 FPE fix for bug 1676 | 67 // | 67 // | | 68 // ******************************************* 68 // ********************************************************************* 69 69 70 #include <limits> 70 #include <limits> 71 #include "G4LowEPComptonModel.hh" 71 #include "G4LowEPComptonModel.hh" 72 #include "G4PhysicalConstants.hh" 72 #include "G4PhysicalConstants.hh" 73 #include "G4SystemOfUnits.hh" 73 #include "G4SystemOfUnits.hh" 74 #include "G4Electron.hh" 74 #include "G4Electron.hh" 75 #include "G4ParticleChangeForGamma.hh" 75 #include "G4ParticleChangeForGamma.hh" 76 #include "G4LossTableManager.hh" 76 #include "G4LossTableManager.hh" 77 #include "G4VAtomDeexcitation.hh" 77 #include "G4VAtomDeexcitation.hh" 78 #include "G4AtomicShell.hh" 78 #include "G4AtomicShell.hh" 79 #include "G4AutoLock.hh" << 80 #include "G4Gamma.hh" 79 #include "G4Gamma.hh" 81 #include "G4ShellData.hh" 80 #include "G4ShellData.hh" 82 #include "G4DopplerProfile.hh" 81 #include "G4DopplerProfile.hh" 83 #include "G4Log.hh" 82 #include "G4Log.hh" 84 #include "G4Exp.hh" 83 #include "G4Exp.hh" 85 84 86 //******************************************** 85 //**************************************************************************** 87 86 88 using namespace std; 87 using namespace std; 89 namespace { G4Mutex LowEPComptonModelMutex = G << 90 88 91 G4PhysicsFreeVector* G4LowEPComptonModel::data << 89 G4int G4LowEPComptonModel::maxZ = 99; 92 G4ShellData* G4LowEPComptonModel::shellD << 90 G4LPhysicsFreeVector* G4LowEPComptonModel::data[] = {0}; 93 G4DopplerProfile* G4LowEPComptonModel::profil << 91 G4ShellData* G4LowEPComptonModel::shellData = 0; >> 92 G4DopplerProfile* G4LowEPComptonModel::profileData = 0; 94 93 95 static const G4double ln10 = G4Log(10.); 94 static const G4double ln10 = G4Log(10.); 96 95 97 G4LowEPComptonModel::G4LowEPComptonModel(const 96 G4LowEPComptonModel::G4LowEPComptonModel(const G4ParticleDefinition*, 98 97 const G4String& nam) 99 : G4VEmModel(nam),isInitialised(false) 98 : G4VEmModel(nam),isInitialised(false) 100 { 99 { 101 verboseLevel=1 ; 100 verboseLevel=1 ; 102 // Verbosity scale: 101 // Verbosity scale: 103 // 0 = nothing 102 // 0 = nothing 104 // 1 = warning for energy non-conservation 103 // 1 = warning for energy non-conservation 105 // 2 = details of energy budget 104 // 2 = details of energy budget 106 // 3 = calculation of cross sections, file o 105 // 3 = calculation of cross sections, file openings, sampling of atoms 107 // 4 = entering in methods 106 // 4 = entering in methods 108 107 109 if( verboseLevel>1 ) { 108 if( verboseLevel>1 ) { 110 G4cout << "Low energy photon Compton model 109 G4cout << "Low energy photon Compton model is constructed " << G4endl; 111 } 110 } 112 111 113 //Mark this model as "applicable" for atomic 112 //Mark this model as "applicable" for atomic deexcitation 114 SetDeexcitationFlag(true); 113 SetDeexcitationFlag(true); 115 114 116 fParticleChange = nullptr; << 115 fParticleChange = 0; 117 fAtomDeexcitation = nullptr; << 116 fAtomDeexcitation = 0; 118 } 117 } 119 118 120 //******************************************** 119 //**************************************************************************** 121 120 122 G4LowEPComptonModel::~G4LowEPComptonModel() 121 G4LowEPComptonModel::~G4LowEPComptonModel() 123 { 122 { 124 if(IsMaster()) { 123 if(IsMaster()) { 125 delete shellData; 124 delete shellData; 126 shellData = nullptr; << 125 shellData = 0; 127 delete profileData; 126 delete profileData; 128 profileData = nullptr; << 127 profileData = 0; 129 } 128 } 130 } 129 } 131 130 132 //******************************************** 131 //**************************************************************************** 133 132 134 void G4LowEPComptonModel::Initialise(const G4P 133 void G4LowEPComptonModel::Initialise(const G4ParticleDefinition* particle, 135 const G4D << 134 const G4DataVector& cuts) 136 { 135 { 137 if (verboseLevel > 1) { 136 if (verboseLevel > 1) { 138 G4cout << "Calling G4LowEPComptonModel::In 137 G4cout << "Calling G4LowEPComptonModel::Initialise()" << G4endl; 139 } 138 } 140 139 141 // Initialise element selector 140 // Initialise element selector >> 141 142 if(IsMaster()) { 142 if(IsMaster()) { 143 143 144 // Access to elements 144 // Access to elements 145 const char* path = G4FindDataDir("G4LEDATA << 145 >> 146 char* path = getenv("G4LEDATA"); 146 147 147 G4ProductionCutsTable* theCoupleTable = 148 G4ProductionCutsTable* theCoupleTable = 148 G4ProductionCutsTable::GetProductionCuts 149 G4ProductionCutsTable::GetProductionCutsTable(); 149 G4int numOfCouples = (G4int)theCoupleTable << 150 G4int numOfCouples = theCoupleTable->GetTableSize(); 150 151 151 for(G4int i=0; i<numOfCouples; ++i) { 152 for(G4int i=0; i<numOfCouples; ++i) { 152 const G4Material* material = 153 const G4Material* material = 153 theCoupleTable->GetMaterialCutsCouple( 154 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial(); 154 const G4ElementVector* theElementVector 155 const G4ElementVector* theElementVector = material->GetElementVector(); 155 std::size_t nelm = material->GetNumberOf << 156 G4int nelm = material->GetNumberOfElements(); 156 157 157 for (std::size_t j=0; j<nelm; ++j) { << 158 for (G4int j=0; j<nelm; ++j) { 158 G4int Z = G4lrint((*theElementVector)[ 159 G4int Z = G4lrint((*theElementVector)[j]->GetZ()); 159 if(Z < 1) { Z = 1; } 160 if(Z < 1) { Z = 1; } 160 else if(Z > maxZ){ Z = maxZ; } 161 else if(Z > maxZ){ Z = maxZ; } 161 162 162 if( (!data[Z]) ) { ReadData(Z, path); 163 if( (!data[Z]) ) { ReadData(Z, path); } 163 } 164 } 164 } 165 } 165 166 166 // For Doppler broadening 167 // For Doppler broadening 167 if(!shellData) { 168 if(!shellData) { 168 shellData = new G4ShellData(); 169 shellData = new G4ShellData(); 169 shellData->SetOccupancyData(); 170 shellData->SetOccupancyData(); 170 G4String file = "/doppler/shell-doppler" 171 G4String file = "/doppler/shell-doppler"; 171 shellData->LoadData(file); 172 shellData->LoadData(file); 172 } 173 } 173 if(!profileData) { profileData = new G4Dop 174 if(!profileData) { profileData = new G4DopplerProfile(); } 174 175 175 InitialiseElementSelectors(particle, cuts) 176 InitialiseElementSelectors(particle, cuts); 176 } 177 } 177 178 178 if (verboseLevel > 2) { 179 if (verboseLevel > 2) { 179 G4cout << "Loaded cross section files" << 180 G4cout << "Loaded cross section files" << G4endl; 180 } 181 } 181 182 182 if( verboseLevel>1 ) { 183 if( verboseLevel>1 ) { 183 G4cout << "G4LowEPComptonModel is initiali 184 G4cout << "G4LowEPComptonModel is initialized " << G4endl 184 << "Energy range: " 185 << "Energy range: " 185 << LowEnergyLimit() / eV << " eV - 186 << LowEnergyLimit() / eV << " eV - " 186 << HighEnergyLimit() / GeV << " GeV 187 << HighEnergyLimit() / GeV << " GeV" 187 << G4endl; 188 << G4endl; 188 } 189 } 189 190 190 if(isInitialised) { return; } 191 if(isInitialised) { return; } 191 192 192 fParticleChange = GetParticleChangeForGamma( 193 fParticleChange = GetParticleChangeForGamma(); 193 fAtomDeexcitation = G4LossTableManager::Ins 194 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation(); 194 isInitialised = true; 195 isInitialised = true; 195 } 196 } 196 197 197 //******************************************** 198 //**************************************************************************** 198 199 199 void G4LowEPComptonModel::InitialiseLocal(cons 200 void G4LowEPComptonModel::InitialiseLocal(const G4ParticleDefinition*, 200 201 G4VEmModel* masterModel) 201 { 202 { 202 SetElementSelectors(masterModel->GetElementS 203 SetElementSelectors(masterModel->GetElementSelectors()); 203 } 204 } 204 205 205 //******************************************** 206 //**************************************************************************** 206 207 207 void G4LowEPComptonModel::ReadData(std::size_t << 208 void G4LowEPComptonModel::ReadData(size_t Z, const char* path) 208 { 209 { 209 if (verboseLevel > 1) 210 if (verboseLevel > 1) 210 { 211 { 211 G4cout << "G4LowEPComptonModel::ReadData() 212 G4cout << "G4LowEPComptonModel::ReadData()" 212 << G4endl; 213 << G4endl; 213 } 214 } 214 if(data[Z]) { return; } 215 if(data[Z]) { return; } 215 const char* datadir = path; 216 const char* datadir = path; 216 if(!datadir) 217 if(!datadir) 217 { 218 { 218 datadir = G4FindDataDir("G4LEDATA"); << 219 datadir = getenv("G4LEDATA"); 219 if(!datadir) 220 if(!datadir) 220 { 221 { 221 G4Exception("G4LowEPComptonModel::ReadDa 222 G4Exception("G4LowEPComptonModel::ReadData()", 222 "em0006",FatalException, 223 "em0006",FatalException, 223 "Environment variable G4LEDA 224 "Environment variable G4LEDATA not defined"); 224 return; 225 return; 225 } 226 } 226 } 227 } 227 228 228 data[Z] = new G4PhysicsFreeVector(); << 229 data[Z] = new G4LPhysicsFreeVector(); >> 230 >> 231 // Activation of spline interpolation >> 232 data[Z]->SetSpline(false); 229 233 230 std::ostringstream ost; 234 std::ostringstream ost; 231 ost << datadir << "/livermore/comp/ce-cs-" < 235 ost << datadir << "/livermore/comp/ce-cs-" << Z <<".dat"; 232 std::ifstream fin(ost.str().c_str()); 236 std::ifstream fin(ost.str().c_str()); 233 237 234 if( !fin.is_open()) 238 if( !fin.is_open()) 235 { 239 { 236 G4ExceptionDescription ed; 240 G4ExceptionDescription ed; 237 ed << "G4LowEPComptonModel data file <" 241 ed << "G4LowEPComptonModel data file <" << ost.str().c_str() 238 << "> is not opened!" << G4endl; 242 << "> is not opened!" << G4endl; 239 G4Exception("G4LowEPComptonModel::ReadData 243 G4Exception("G4LowEPComptonModel::ReadData()", 240 "em0003",FatalException, 244 "em0003",FatalException, 241 ed,"G4LEDATA version should be 245 ed,"G4LEDATA version should be G4EMLOW6.34 or later"); 242 return; 246 return; 243 } else { 247 } else { 244 if(verboseLevel > 3) { 248 if(verboseLevel > 3) { 245 G4cout << "File " << ost.str() 249 G4cout << "File " << ost.str() 246 << " is opened by G4LowEPCompto 250 << " is opened by G4LowEPComptonModel" << G4endl; 247 } 251 } 248 data[Z]->Retrieve(fin, true); 252 data[Z]->Retrieve(fin, true); 249 data[Z]->ScaleVector(MeV, MeV*barn); 253 data[Z]->ScaleVector(MeV, MeV*barn); 250 } 254 } 251 fin.close(); 255 fin.close(); 252 } 256 } 253 257 254 //******************************************** 258 //**************************************************************************** 255 259 256 260 257 G4double 261 G4double 258 G4LowEPComptonModel::ComputeCrossSectionPerAto 262 G4LowEPComptonModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*, 259 263 G4double GammaEnergy, 260 264 G4double Z, G4double, 261 265 G4double, G4double) 262 { 266 { 263 if (verboseLevel > 3) { 267 if (verboseLevel > 3) { 264 G4cout << "G4LowEPComptonModel::ComputeCro 268 G4cout << "G4LowEPComptonModel::ComputeCrossSectionPerAtom()" 265 << G4endl; 269 << G4endl; 266 } 270 } 267 G4double cs = 0.0; 271 G4double cs = 0.0; 268 272 269 if (GammaEnergy < LowEnergyLimit()) { return 273 if (GammaEnergy < LowEnergyLimit()) { return 0.0; } 270 274 271 G4int intZ = G4lrint(Z); 275 G4int intZ = G4lrint(Z); 272 if(intZ < 1 || intZ > maxZ) { return cs; } 276 if(intZ < 1 || intZ > maxZ) { return cs; } 273 277 274 G4PhysicsFreeVector* pv = data[intZ]; << 278 G4LPhysicsFreeVector* pv = data[intZ]; 275 279 276 // if element was not initialised 280 // if element was not initialised 277 // do initialisation safely for MT mode 281 // do initialisation safely for MT mode 278 if(!pv) 282 if(!pv) 279 { 283 { 280 InitialiseForElement(0, intZ); 284 InitialiseForElement(0, intZ); 281 pv = data[intZ]; 285 pv = data[intZ]; 282 if(!pv) { return cs; } 286 if(!pv) { return cs; } 283 } 287 } 284 288 285 G4int n = G4int(pv->GetVectorLength() - 1); << 289 G4int n = pv->GetVectorLength() - 1; 286 G4double e1 = pv->Energy(0); 290 G4double e1 = pv->Energy(0); 287 G4double e2 = pv->Energy(n); 291 G4double e2 = pv->Energy(n); 288 292 289 if(GammaEnergy <= e1) { cs = GammaEnerg 293 if(GammaEnergy <= e1) { cs = GammaEnergy/(e1*e1)*pv->Value(e1); } 290 else if(GammaEnergy <= e2) { cs = pv->Value( 294 else if(GammaEnergy <= e2) { cs = pv->Value(GammaEnergy)/GammaEnergy; } 291 else if(GammaEnergy > e2) { cs = pv->Value( 295 else if(GammaEnergy > e2) { cs = pv->Value(e2)/GammaEnergy; } 292 296 293 return cs; 297 return cs; 294 } 298 } 295 299 296 //******************************************** 300 //**************************************************************************** 297 301 298 void G4LowEPComptonModel::SampleSecondaries(st 302 void G4LowEPComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, 299 303 const G4MaterialCutsCouple* couple, 300 304 const G4DynamicParticle* aDynamicGamma, 301 305 G4double, G4double) 302 { 306 { 303 307 304 // The scattered gamma energy is sampled acc 308 // The scattered gamma energy is sampled according to Klein - Nishina formula. 305 // then accepted or rejected depending on th 309 // then accepted or rejected depending on the Scattering Function multiplied 306 // by factor from Klein - Nishina formula. 310 // by factor from Klein - Nishina formula. 307 // Expression of the angular distribution as 311 // Expression of the angular distribution as Klein Nishina 308 // angular and energy distribution and Scatt 312 // angular and energy distribution and Scattering fuctions is taken from 309 // D. E. Cullen "A simple model of photon tr 313 // D. E. Cullen "A simple model of photon transport" Nucl. Instr. Meth. 310 // Phys. Res. B 101 (1995). Method of sampli 314 // Phys. Res. B 101 (1995). Method of sampling with form factors is different 311 // data are interpolated while in the articl 315 // data are interpolated while in the article they are fitted. 312 // Reference to the article is from J. Stepa 316 // Reference to the article is from J. Stepanek New Photon, Positron 313 // and Electron Interaction Data for GEANT i 317 // and Electron Interaction Data for GEANT in Energy Range from 1 eV to 10 314 // TeV (draft). 318 // TeV (draft). 315 // The random number techniques of Butcher & 319 // The random number techniques of Butcher & Messel are used 316 // (Nucl Phys 20(1960),15). 320 // (Nucl Phys 20(1960),15). 317 321 >> 322 318 G4double photonEnergy0 = aDynamicGamma->GetK 323 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy()/MeV; 319 324 320 if (verboseLevel > 3) { 325 if (verboseLevel > 3) { 321 G4cout << "G4LowEPComptonModel::SampleSeco 326 G4cout << "G4LowEPComptonModel::SampleSecondaries() E(MeV)= " 322 << photonEnergy0/MeV << " in " << c 327 << photonEnergy0/MeV << " in " << couple->GetMaterial()->GetName() 323 << G4endl; 328 << G4endl; 324 } 329 } 325 // do nothing below the threshold 330 // do nothing below the threshold 326 // should never get here because the XS is z 331 // should never get here because the XS is zero below the limit 327 if (photonEnergy0 < LowEnergyLimit()) 332 if (photonEnergy0 < LowEnergyLimit()) 328 return ; 333 return ; 329 334 330 G4double e0m = photonEnergy0 / electron_mass 335 G4double e0m = photonEnergy0 / electron_mass_c2 ; 331 G4ParticleMomentum photonDirection0 = aDynam 336 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection(); 332 337 333 // Select randomly one element in the curren 338 // Select randomly one element in the current material >> 339 334 const G4ParticleDefinition* particle = aDyn 340 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition(); 335 const G4Element* elm = SelectRandomAtom(coup 341 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0); 336 G4int Z = (G4int)elm->GetZ(); 342 G4int Z = (G4int)elm->GetZ(); 337 343 338 G4double LowEPCepsilon0 = 1. / (1. + 2. * e0 344 G4double LowEPCepsilon0 = 1. / (1. + 2. * e0m); 339 G4double LowEPCepsilon0Sq = LowEPCepsilon0 * 345 G4double LowEPCepsilon0Sq = LowEPCepsilon0 * LowEPCepsilon0; 340 G4double alpha1 = -std::log(LowEPCepsilon0); 346 G4double alpha1 = -std::log(LowEPCepsilon0); 341 G4double alpha2 = 0.5 * (1. - LowEPCepsilon0 347 G4double alpha2 = 0.5 * (1. - LowEPCepsilon0Sq); 342 348 343 G4double wlPhoton = h_Planck*c_light/photonE 349 G4double wlPhoton = h_Planck*c_light/photonEnergy0; 344 350 345 // Sample the energy of the scattered photon 351 // Sample the energy of the scattered photon 346 G4double LowEPCepsilon; 352 G4double LowEPCepsilon; 347 G4double LowEPCepsilonSq; 353 G4double LowEPCepsilonSq; 348 G4double oneCosT; 354 G4double oneCosT; 349 G4double sinT2; 355 G4double sinT2; 350 G4double gReject; 356 G4double gReject; 351 357 352 if (verboseLevel > 3) { 358 if (verboseLevel > 3) { 353 G4cout << "Started loop to sample gamma en 359 G4cout << "Started loop to sample gamma energy" << G4endl; 354 } 360 } 355 361 356 do 362 do 357 { 363 { 358 if ( alpha1/(alpha1+alpha2) > G4UniformR 364 if ( alpha1/(alpha1+alpha2) > G4UniformRand()) 359 { 365 { 360 LowEPCepsilon = G4Exp(-alpha1 * G4Un 366 LowEPCepsilon = G4Exp(-alpha1 * G4UniformRand()); 361 LowEPCepsilonSq = LowEPCepsilon * Lo 367 LowEPCepsilonSq = LowEPCepsilon * LowEPCepsilon; 362 } 368 } 363 else 369 else 364 { 370 { 365 LowEPCepsilonSq = LowEPCepsilon0Sq + 371 LowEPCepsilonSq = LowEPCepsilon0Sq + (1. - LowEPCepsilon0Sq) * G4UniformRand(); 366 LowEPCepsilon = std::sqrt(LowEPCepsi 372 LowEPCepsilon = std::sqrt(LowEPCepsilonSq); 367 } 373 } 368 374 369 oneCosT = (1. - LowEPCepsilon) / ( LowEP 375 oneCosT = (1. - LowEPCepsilon) / ( LowEPCepsilon * e0m); 370 sinT2 = oneCosT * (2. - oneCosT); 376 sinT2 = oneCosT * (2. - oneCosT); 371 G4double x = std::sqrt(oneCosT/2.) / (wl 377 G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm); 372 G4double scatteringFunction = ComputeSca 378 G4double scatteringFunction = ComputeScatteringFunction(x, Z); 373 gReject = (1. - LowEPCepsilon * sinT2 / 379 gReject = (1. - LowEPCepsilon * sinT2 / (1. + LowEPCepsilonSq)) * scatteringFunction; 374 380 375 } while(gReject < G4UniformRand()*Z); 381 } while(gReject < G4UniformRand()*Z); 376 382 377 G4double cosTheta = 1. - oneCosT; 383 G4double cosTheta = 1. - oneCosT; 378 G4double sinTheta = std::sqrt(sinT2); 384 G4double sinTheta = std::sqrt(sinT2); 379 G4double phi = twopi * G4UniformRand(); 385 G4double phi = twopi * G4UniformRand(); 380 G4double dirx = sinTheta * std::cos(phi); 386 G4double dirx = sinTheta * std::cos(phi); 381 G4double diry = sinTheta * std::sin(phi); 387 G4double diry = sinTheta * std::sin(phi); 382 G4double dirz = cosTheta ; 388 G4double dirz = cosTheta ; 383 389 >> 390 384 // Scatter photon energy and Compton electro 391 // Scatter photon energy and Compton electron direction - Method based on: 385 // J. M. C. Brown, M. R. Dimmock, J. E. Gill 392 // J. M. C. Brown, M. R. Dimmock, J. E. Gillam and D. M. Paganin' 386 // "A low energy bound atomic electron Compt 393 // "A low energy bound atomic electron Compton scattering model for Geant4" 387 // NIMB, Vol. 338, 77-88, 2014. 394 // NIMB, Vol. 338, 77-88, 2014. 388 395 389 // Set constants and initialize scattering p 396 // Set constants and initialize scattering parameters >> 397 390 const G4double vel_c = c_light / (m/s); 398 const G4double vel_c = c_light / (m/s); 391 const G4double momentum_au_to_nat = halfpi* 399 const G4double momentum_au_to_nat = halfpi* hbar_Planck / Bohr_radius / (kg*m/s); 392 const G4double e_mass_kg = electron_mass_c2 400 const G4double e_mass_kg = electron_mass_c2 / c_squared / kg ; 393 401 394 const G4int maxDopplerIterations = 1000; 402 const G4int maxDopplerIterations = 1000; 395 G4double bindingE = 0.; 403 G4double bindingE = 0.; 396 G4double pEIncident = photonEnergy0 ; 404 G4double pEIncident = photonEnergy0 ; 397 G4double pERecoil = -1.; 405 G4double pERecoil = -1.; 398 G4double eERecoil = -1.; 406 G4double eERecoil = -1.; 399 G4double e_alpha =0.; 407 G4double e_alpha =0.; 400 G4double e_beta = 0.; 408 G4double e_beta = 0.; 401 409 402 G4double CE_emission_flag = 0.; 410 G4double CE_emission_flag = 0.; 403 G4double ePAU = -1; 411 G4double ePAU = -1; 404 G4int shellIdx = 0; 412 G4int shellIdx = 0; 405 G4double u_temp = 0; 413 G4double u_temp = 0; 406 G4double cosPhiE =0; 414 G4double cosPhiE =0; 407 G4double sinThetaE =0; 415 G4double sinThetaE =0; 408 G4double cosThetaE =0; 416 G4double cosThetaE =0; 409 G4int iteration = 0; 417 G4int iteration = 0; 410 418 411 if (verboseLevel > 3) { 419 if (verboseLevel > 3) { 412 G4cout << "Started loop to sample photon e 420 G4cout << "Started loop to sample photon energy and electron direction" << G4endl; 413 } 421 } 414 422 415 do{ 423 do{ >> 424 >> 425 416 // ************************************* 426 // ****************************************** 417 // | Determine scatter photon energy 427 // | Determine scatter photon energy | 418 // ************************************* 428 // ****************************************** 419 do << 420 { << 421 iteration++; << 422 // ***************************************** << 423 // | Sample bound electron information << 424 // ***************************************** << 425 << 426 // Select shell based on shell occupancy << 427 shellIdx = shellData->SelectRandomShell(Z); << 428 bindingE = shellData->BindingEnergy(Z,shellI << 429 << 430 // Randomly sample bound electron momentum ( << 431 ePAU = profileData->RandomSelectMomentum(Z,s << 432 << 433 // Convert to SI units << 434 G4double ePSI = ePAU * momentum_au_to_nat; << 435 << 436 //Calculate bound electron velocity and norm << 437 u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / << 438 << 439 // Sample incident electron direction, amorp << 440 e_alpha = pi*G4UniformRand(); << 441 e_beta = twopi*G4UniformRand(); << 442 << 443 // Total energy of system << 444 << 445 G4double eEIncident = electron_mass_c2 / sqr << 446 G4double systemE = eEIncident + pEIncident; << 447 G4double gamma_temp = 1.0 / sqrt( 1 - (u_tem << 448 G4double numerator = gamma_temp*electron_mas << 449 G4double subdenom1 = u_temp*cosTheta*std::c << 450 G4double subdenom2 = u_temp*sinTheta*std::si << 451 G4double denominator = (1.0 - cosTheta) + ( << 452 pERecoil = (numerator/denominator); << 453 eERecoil = systemE - pERecoil; << 454 CE_emission_flag = pEIncident - pERecoil; << 455 } while ( (iteration <= maxDopplerIterat << 456 << 457 // End of recalculation of photon energy w << 458 << 459 // *************************************** << 460 // | Determine ejected Compton electro << 461 // *************************************** << 462 << 463 // Calculate velocity of ejected Compton e << 464 << 465 G4double a_temp = eERecoil / electron_mass << 466 G4double u_p_temp = sqrt(1 - (1 / (a_temp* << 467 << 468 // Coefficients and terms from simulatenou << 469 G4double sinAlpha = std::sin(e_alpha); << 470 G4double cosAlpha = std::cos(e_alpha); << 471 G4double sinBeta = std::sin(e_beta); << 472 G4double cosBeta = std::cos(e_beta); << 473 << 474 G4double gamma = 1.0 / sqrt(1 - (u_temp*u_ << 475 G4double gamma_p = 1.0 / sqrt(1 - (u_p_tem << 476 << 477 G4double var_A = pERecoil*u_p_temp*sinThet << 478 G4double var_B = u_p_temp* (pERecoil*cosTh << 479 G4double var_C = (pERecoil-pEIncident) - ( << 480 << 481 G4double var_D1 = gamma*electron_mass_c2*p << 482 G4double var_D2 = (1 - (u_temp*cosTheta*co << 483 G4double var_D3 = ((electron_mass_c2*elect << 484 G4double var_D = var_D1*var_D2 + var_D3; << 485 << 486 G4double var_E1 = ((gamma*gamma_p)*(electr << 487 G4double var_E2 = gamma_p*electron_mass_c2 << 488 G4double var_E = var_E1 - var_E2; << 489 << 490 G4double var_F1 = ((gamma*gamma_p)*(electr << 491 G4double var_F2 = (gamma_p*electron_mass_c << 492 G4double var_F = var_F1 - var_F2; << 493 << 494 G4double var_G = (gamma*gamma_p)*(electron << 495 << 496 // Two equations form a quadratic form of << 497 // Coefficents and solution to quadratic << 498 G4double var_W1 = (var_F*var_B - var_E*var << 499 G4double var_W2 = (var_G*var_G)*(var_A*var << 500 G4double var_W = var_W1 + var_W2; << 501 << 502 G4double var_Y = 2.0*(((var_A*var_D-var_F* << 503 << 504 G4double var_Z1 = (var_A*var_D - var_F*var << 505 G4double var_Z2 = (var_G*var_G)*(var_C*var << 506 G4double var_Z = var_Z1 + var_Z2; << 507 G4double diff1 = var_Y*var_Y; << 508 G4double diff2 = 4*var_W*var_Z; << 509 G4double diff = diff1 - diff2; << 510 429 511 // Check if diff is less than zero, if so << 430 do >> 431 { >> 432 iteration++; >> 433 >> 434 >> 435 // ******************************************** >> 436 // | Sample bound electron information | >> 437 // ******************************************** >> 438 >> 439 // Select shell based on shell occupancy >> 440 >> 441 shellIdx = shellData->SelectRandomShell(Z); >> 442 bindingE = shellData->BindingEnergy(Z,shellIdx)/MeV; >> 443 >> 444 >> 445 // Randomly sample bound electron momentum (memento: the data set is in Atomic Units) >> 446 ePAU = profileData->RandomSelectMomentum(Z,shellIdx); >> 447 >> 448 // Convert to SI units >> 449 G4double ePSI = ePAU * momentum_au_to_nat; >> 450 >> 451 //Calculate bound electron velocity and normalise to natural units >> 452 u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)) )/vel_c; >> 453 >> 454 // Sample incident electron direction, amorphous material, to scattering photon scattering plane >> 455 >> 456 e_alpha = pi*G4UniformRand(); >> 457 e_beta = twopi*G4UniformRand(); >> 458 >> 459 // Total energy of system >> 460 >> 461 G4double eEIncident = electron_mass_c2 / sqrt( 1 - (u_temp*u_temp)); >> 462 G4double systemE = eEIncident + pEIncident; >> 463 >> 464 >> 465 G4double gamma_temp = 1.0 / sqrt( 1 - (u_temp*u_temp)); >> 466 G4double numerator = gamma_temp*electron_mass_c2*(1 - u_temp * std::cos(e_alpha)); >> 467 G4double subdenom1 = u_temp*cosTheta*std::cos(e_alpha); >> 468 G4double subdenom2 = u_temp*sinTheta*std::sin(e_alpha)*std::cos(e_beta); >> 469 G4double denominator = (1.0 - cosTheta) + (gamma_temp*electron_mass_c2*(1 - subdenom1 - subdenom2) / pEIncident); >> 470 pERecoil = (numerator/denominator); >> 471 eERecoil = systemE - pERecoil; >> 472 CE_emission_flag = pEIncident - pERecoil; >> 473 } while ( (iteration <= maxDopplerIterations) && (CE_emission_flag < bindingE)); >> 474 >> 475 // End of recalculation of photon energy with Doppler broadening >> 476 >> 477 >> 478 >> 479 // ******************************************************* >> 480 // | Determine ejected Compton electron direction | >> 481 // ******************************************************* >> 482 >> 483 // Calculate velocity of ejected Compton electron >> 484 >> 485 G4double a_temp = eERecoil / electron_mass_c2; >> 486 G4double u_p_temp = sqrt(1 - (1 / (a_temp*a_temp))); >> 487 >> 488 // Coefficients and terms from simulatenous equations >> 489 >> 490 G4double sinAlpha = std::sin(e_alpha); >> 491 G4double cosAlpha = std::cos(e_alpha); >> 492 G4double sinBeta = std::sin(e_beta); >> 493 G4double cosBeta = std::cos(e_beta); >> 494 >> 495 G4double gamma = 1.0 / sqrt(1 - (u_temp*u_temp)); >> 496 G4double gamma_p = 1.0 / sqrt(1 - (u_p_temp*u_p_temp)); >> 497 >> 498 G4double var_A = pERecoil*u_p_temp*sinTheta; >> 499 G4double var_B = u_p_temp* (pERecoil*cosTheta-pEIncident); >> 500 G4double var_C = (pERecoil-pEIncident) - ( (pERecoil*pEIncident) / (gamma_p*electron_mass_c2))*(1 - cosTheta); >> 501 >> 502 G4double var_D1 = gamma*electron_mass_c2*pERecoil; >> 503 G4double var_D2 = (1 - (u_temp*cosTheta*cosAlpha) - (u_temp*sinTheta*cosBeta*sinAlpha)); >> 504 G4double var_D3 = ((electron_mass_c2*electron_mass_c2)*(gamma*gamma_p - 1)) - (gamma_p*electron_mass_c2*pERecoil); >> 505 G4double var_D = var_D1*var_D2 + var_D3; >> 506 >> 507 G4double var_E1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosAlpha); >> 508 G4double var_E2 = gamma_p*electron_mass_c2*pERecoil*u_p_temp*cosTheta; >> 509 G4double var_E = var_E1 - var_E2; >> 510 >> 511 G4double var_F1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosBeta*sinAlpha); >> 512 G4double var_F2 = (gamma_p*electron_mass_c2*pERecoil*u_p_temp*sinTheta); >> 513 G4double var_F = var_F1 - var_F2; >> 514 >> 515 G4double var_G = (gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*sinBeta*sinAlpha; >> 516 >> 517 // Two equations form a quadratic form of Wx^2 + Yx + Z = 0 >> 518 // Coefficents and solution to quadratic >> 519 >> 520 G4double var_W1 = (var_F*var_B - var_E*var_A)*(var_F*var_B - var_E*var_A); >> 521 G4double var_W2 = (var_G*var_G)*(var_A*var_A) + (var_G*var_G)*(var_B*var_B); >> 522 G4double var_W = var_W1 + var_W2; >> 523 >> 524 G4double var_Y = 2.0*(((var_A*var_D-var_F*var_C)*(var_F*var_B-var_E*var_A)) - ((var_G*var_G)*var_B*var_C)); >> 525 >> 526 G4double var_Z1 = (var_A*var_D - var_F*var_C)*(var_A*var_D - var_F*var_C); >> 527 G4double var_Z2 = (var_G*var_G)*(var_C*var_C) - (var_G*var_G)*(var_A*var_A); >> 528 G4double var_Z = var_Z1 + var_Z2; >> 529 G4double diff1 = var_Y*var_Y; >> 530 G4double diff2 = 4*var_W*var_Z; >> 531 G4double diff = diff1 - diff2; >> 532 >> 533 >> 534 // Check if diff is less than zero, if so ensure it is due to FPE >> 535 512 //Determine number of digits (in decimal b 536 //Determine number of digits (in decimal base) that G4double can accurately represent 513 G4double g4d_order = G4double(numeric_limi << 537 G4double g4d_order = G4double(numeric_limits<G4double>::digits10); 514 G4double g4d_limit = std::pow(10.,-g4d_ord << 538 G4double g4d_limit = std::pow(10.,-g4d_order); 515 //Confirm that diff less than zero is due << 539 //Confirm that diff less than zero is due FPE, i.e if abs of diff / diff1 and diff/ diff2 is less 516 //than 10^(-g4d_order), then set diff to z << 540 //than 10^(-g4d_order), then set diff to zero 517 << 541 518 if ((diff < 0.0) && (abs(diff / diff1) < g << 542 if ((diff < 0.0) && (abs(diff / diff1) < g4d_limit) && (abs(diff / diff2) < g4d_limit) ) 519 { << 543 { 520 diff = 0.0; << 544 diff = 0.0; 521 } << 545 } 522 << 546 523 << 547 524 // Plus and minus of quadratic << 548 // Plus and minus of quadratic 525 G4double X_p = (-var_Y + sqrt (diff))/(2*v << 549 G4double X_p = (-var_Y + sqrt (diff))/(2*var_W); 526 G4double X_m = (-var_Y - sqrt (diff))/(2*v << 550 G4double X_m = (-var_Y - sqrt (diff))/(2*var_W); 527 << 551 528 // Floating point precision protection << 552 529 // Check if X_p and X_m are greater than o << 553 // Floating point precision protection 530 // Issue due to propagation of FPE and onl << 554 // Check if X_p and X_m are greater than or less than 1 or -1, if so clean up FPE 531 if(X_p >1){X_p=1;} if(X_p<-1){X_p=-1;} << 555 // Issue due to propagation of FPE and only impacts 8th sig fig onwards 532 if(X_m >1){X_m=1;} if(X_m<-1){X_m=-1;} << 556 533 // End of FP protection << 557 if(X_p >1){X_p=1;} if(X_p<-1){X_p=-1;} 534 << 558 if(X_m >1){X_m=1;} if(X_m<-1){X_m=-1;} 535 G4double ThetaE = 0.; << 559 536 // Randomly sample one of the two possible << 560 // End of FP protection 537 G4double sol_select = G4UniformRand(); << 561 538 << 562 G4double ThetaE = 0.; 539 if (sol_select < 0.5) << 563 >> 564 >> 565 // Randomly sample one of the two possible solutions and determin theta angle of ejected Compton electron >> 566 G4double sol_select = G4UniformRand(); >> 567 >> 568 if (sol_select < 0.5) 540 { 569 { 541 ThetaE = std::acos(X_p); << 570 ThetaE = std::acos(X_p); 542 } 571 } 543 if (sol_select > 0.5) << 572 if (sol_select > 0.5) 544 { 573 { 545 ThetaE = std::acos(X_m); << 574 ThetaE = std::acos(X_m); 546 } 575 } 547 cosThetaE = std::cos(ThetaE); << 576 548 sinThetaE = std::sin(ThetaE); << 577 549 G4double Theta = std::acos(cosTheta); << 578 cosThetaE = std::cos(ThetaE); 550 << 579 sinThetaE = std::sin(ThetaE); 551 //Calculate electron Phi << 580 G4double Theta = std::acos(cosTheta); 552 G4double iSinThetaE = std::sqrt(1+std::tan << 581 553 G4double iSinTheta = std::sqrt(1+std::tan( << 582 //Calculate electron Phi 554 G4double ivar_A = iSinTheta/ (pERecoil*u_p << 583 G4double iSinThetaE = std::sqrt(1+std::tan((pi/2.0)-ThetaE)*std::tan((pi/2.0)-ThetaE)); 555 // Trigs << 584 G4double iSinTheta = std::sqrt(1+std::tan((pi/2.0)-Theta)*std::tan((pi/2.0)-Theta)); 556 cosPhiE = (var_C - var_B*cosThetaE)*(ivar_ << 585 G4double ivar_A = iSinTheta/ (pERecoil*u_p_temp); 557 << 586 // Trigs 558 // End of calculation of ejection Compton << 587 cosPhiE = (var_C - var_B*cosThetaE)*(ivar_A*iSinThetaE); 559 //Fix for floating point errors << 588 560 << 589 // End of calculation of ejection Compton electron direction 561 } while ( (iteration <= maxDopplerIterations << 590 562 << 591 //Fix for floating point errors 563 // Revert to original if maximum number of i << 592 >> 593 } while ( (iteration <= maxDopplerIterations) && (abs(cosPhiE) > 1)); >> 594 >> 595 // Revert to original if maximum number of iterations threshold has been reached 564 if (iteration >= maxDopplerIterations) 596 if (iteration >= maxDopplerIterations) 565 { 597 { 566 pERecoil = photonEnergy0 ; 598 pERecoil = photonEnergy0 ; 567 bindingE = 0.; 599 bindingE = 0.; 568 dirx=0.0; 600 dirx=0.0; 569 diry=0.0; 601 diry=0.0; 570 dirz=1.0; 602 dirz=1.0; 571 } 603 } 572 604 573 // Set "scattered" photon direction and ener 605 // Set "scattered" photon direction and energy >> 606 574 G4ThreeVector photonDirection1(dirx,diry,dir 607 G4ThreeVector photonDirection1(dirx,diry,dirz); 575 photonDirection1.rotateUz(photonDirection0); 608 photonDirection1.rotateUz(photonDirection0); 576 fParticleChange->ProposeMomentumDirection(ph 609 fParticleChange->ProposeMomentumDirection(photonDirection1) ; 577 610 578 if (pERecoil > 0.) 611 if (pERecoil > 0.) 579 { 612 { 580 fParticleChange->SetProposedKineticEnergy 613 fParticleChange->SetProposedKineticEnergy(pERecoil) ; 581 614 582 // Set ejected Compton electron direction 615 // Set ejected Compton electron direction and energy 583 G4double PhiE = std::acos(cosPhiE); 616 G4double PhiE = std::acos(cosPhiE); 584 G4double eDirX = sinThetaE * std::cos(phi 617 G4double eDirX = sinThetaE * std::cos(phi+PhiE); 585 G4double eDirY = sinThetaE * std::sin(phi 618 G4double eDirY = sinThetaE * std::sin(phi+PhiE); 586 G4double eDirZ = cosThetaE; 619 G4double eDirZ = cosThetaE; 587 620 588 G4double eKineticEnergy = pEIncident - pE 621 G4double eKineticEnergy = pEIncident - pERecoil - bindingE; 589 622 590 G4ThreeVector eDirection(eDirX,eDirY,eDir 623 G4ThreeVector eDirection(eDirX,eDirY,eDirZ); 591 eDirection.rotateUz(photonDirection0); 624 eDirection.rotateUz(photonDirection0); 592 G4DynamicParticle* dp = new G4DynamicPart 625 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(), 593 626 eDirection,eKineticEnergy) ; 594 fvect->push_back(dp); 627 fvect->push_back(dp); 595 628 596 } 629 } 597 else 630 else 598 { 631 { 599 fParticleChange->SetProposedKineticEnerg 632 fParticleChange->SetProposedKineticEnergy(0.); 600 fParticleChange->ProposeTrackStatus(fSto 633 fParticleChange->ProposeTrackStatus(fStopAndKill); 601 } 634 } 602 635 603 // sample deexcitation 636 // sample deexcitation 604 // 637 // >> 638 605 if (verboseLevel > 3) { 639 if (verboseLevel > 3) { 606 G4cout << "Started atomic de-excitation " 640 G4cout << "Started atomic de-excitation " << fAtomDeexcitation << G4endl; 607 } 641 } 608 642 609 if(fAtomDeexcitation && iteration < maxDoppl 643 if(fAtomDeexcitation && iteration < maxDopplerIterations) { 610 G4int index = couple->GetIndex(); 644 G4int index = couple->GetIndex(); 611 if(fAtomDeexcitation->CheckDeexcitationAct 645 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) { 612 std::size_t nbefore = fvect->size(); << 646 size_t nbefore = fvect->size(); 613 G4AtomicShellEnumerator as = G4AtomicShe 647 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIdx); 614 const G4AtomicShell* shell = fAtomDeexci 648 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as); 615 fAtomDeexcitation->GenerateParticles(fve 649 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index); 616 std::size_t nafter = fvect->size(); << 650 size_t nafter = fvect->size(); 617 if(nafter > nbefore) { 651 if(nafter > nbefore) { 618 for (std::size_t i=nbefore; i<nafter; << 652 for (size_t i=nbefore; i<nafter; ++i) { 619 //Check if there is enough residual 653 //Check if there is enough residual energy 620 if (bindingE >= ((*fvect)[i])->GetKi 654 if (bindingE >= ((*fvect)[i])->GetKineticEnergy()) 621 { 655 { 622 //Ok, this is a valid secondary: 656 //Ok, this is a valid secondary: keep it 623 bindingE -= ((*fvect)[i])->GetKin 657 bindingE -= ((*fvect)[i])->GetKineticEnergy(); 624 } 658 } 625 else 659 else 626 { 660 { 627 //Invalid secondary: not enough e 661 //Invalid secondary: not enough energy to create it! 628 //Keep its energy in the local de 662 //Keep its energy in the local deposit 629 delete (*fvect)[i]; 663 delete (*fvect)[i]; 630 (*fvect)[i]=nullptr; << 664 (*fvect)[i]=0; 631 } 665 } 632 } 666 } 633 } 667 } 634 } 668 } 635 } 669 } 636 670 637 //This should never happen 671 //This should never happen 638 if(bindingE < 0.0) 672 if(bindingE < 0.0) 639 G4Exception("G4LowEPComptonModel::SampleS 673 G4Exception("G4LowEPComptonModel::SampleSecondaries()", 640 "em2051",FatalException,"Nega 674 "em2051",FatalException,"Negative local energy deposit"); 641 675 642 fParticleChange->ProposeLocalEnergyDeposit(b 676 fParticleChange->ProposeLocalEnergyDeposit(bindingE); >> 677 643 } 678 } 644 679 645 //******************************************** 680 //**************************************************************************** 646 681 647 G4double 682 G4double 648 G4LowEPComptonModel::ComputeScatteringFunction 683 G4LowEPComptonModel::ComputeScatteringFunction(G4double x, G4int Z) 649 { 684 { 650 G4double value = Z; 685 G4double value = Z; 651 if (x <= ScatFuncFitParam[Z][2]) { 686 if (x <= ScatFuncFitParam[Z][2]) { 652 687 653 G4double lgq = G4Log(x)/ln10; 688 G4double lgq = G4Log(x)/ln10; 654 689 655 if (lgq < ScatFuncFitParam[Z][1]) { 690 if (lgq < ScatFuncFitParam[Z][1]) { 656 value = ScatFuncFitParam[Z][3] + lgq*Sca 691 value = ScatFuncFitParam[Z][3] + lgq*ScatFuncFitParam[Z][4]; 657 } else { 692 } else { 658 value = ScatFuncFitParam[Z][5] + lgq*Sca 693 value = ScatFuncFitParam[Z][5] + lgq*ScatFuncFitParam[Z][6] + 659 lgq*lgq*ScatFuncFitParam[Z][7] + lgq*l 694 lgq*lgq*ScatFuncFitParam[Z][7] + lgq*lgq*lgq*ScatFuncFitParam[Z][8]; 660 } 695 } 661 value = G4Exp(value*ln10); 696 value = G4Exp(value*ln10); 662 } 697 } 663 return value; 698 return value; 664 } 699 } 665 700 666 701 667 //******************************************** 702 //**************************************************************************** >> 703 >> 704 #include "G4AutoLock.hh" >> 705 namespace { G4Mutex LowEPComptonModelMutex = G4MUTEX_INITIALIZER; } 668 706 669 void 707 void 670 G4LowEPComptonModel::InitialiseForElement(cons 708 G4LowEPComptonModel::InitialiseForElement(const G4ParticleDefinition*, 671 709 G4int Z) 672 { 710 { 673 G4AutoLock l(&LowEPComptonModelMutex); 711 G4AutoLock l(&LowEPComptonModelMutex); 674 if(!data[Z]) { ReadData(Z); } 712 if(!data[Z]) { ReadData(Z); } 675 l.unlock(); 713 l.unlock(); 676 } 714 } 677 715 678 //******************************************** 716 //**************************************************************************** 679 717 680 //Fitting data to compute scattering function 718 //Fitting data to compute scattering function 681 719 682 const G4double G4LowEPComptonModel::ScatFuncFi 720 const G4double G4LowEPComptonModel::ScatFuncFitParam[101][9] = { 683 { 0, 0., 0., 0., 0., 721 { 0, 0., 0., 0., 0., 0., 0., 0., 0.}, 684 { 1, 6.673, 1.49968E+08, -14.352, 1.999, -143 722 { 1, 6.673, 1.49968E+08, -14.352, 1.999, -143.374, 50.787, -5.951, 0.2304418}, 685 { 2, 6.500, 2.50035E+08, -14.215, 1.970, -53. 723 { 2, 6.500, 2.50035E+08, -14.215, 1.970, -53.649, 13.892, -0.948, 0.006996759}, 686 { 3, 6.551, 3.99945E+08, -13.555, 1.993, -62. 724 { 3, 6.551, 3.99945E+08, -13.555, 1.993, -62.090, 21.462, -2.453, 0.093416}, 687 { 4, 6.500, 5.00035E+08, -13.746, 1.998, -127 725 { 4, 6.500, 5.00035E+08, -13.746, 1.998, -127.906, 46.491, -5.614, 0.2262103}, 688 { 5, 6.500, 5.99791E+08, -13.800, 1.998, -131 726 { 5, 6.500, 5.99791E+08, -13.800, 1.998, -131.153, 47.132, -5.619, 0.2233819}, 689 { 6, 6.708, 6.99842E+08, -13.885, 1.999, -128 727 { 6, 6.708, 6.99842E+08, -13.885, 1.999, -128.143, 45.379, -5.325, 0.2083009}, 690 { 7, 6.685, 7.99834E+08, -13.885, 2.000, -131 728 { 7, 6.685, 7.99834E+08, -13.885, 2.000, -131.048, 46.314, -5.421, 0.2114925}, 691 { 8, 6.669, 7.99834E+08, -13.962, 2.001, -128 729 { 8, 6.669, 7.99834E+08, -13.962, 2.001, -128.225, 44.818, -5.183, 0.1997155}, 692 { 9, 6.711, 7.99834E+08, -13.999, 2.000, -122 730 { 9, 6.711, 7.99834E+08, -13.999, 2.000, -122.112, 42.103, -4.796, 0.1819099}, 693 { 10, 6.702, 7.99834E+08, -14.044, 1.999, -110 731 { 10, 6.702, 7.99834E+08, -14.044, 1.999, -110.143, 37.225, -4.143, 0.1532094}, 694 { 11, 6.425, 1.00000E+09, -13.423, 1.993, -41. 732 { 11, 6.425, 1.00000E+09, -13.423, 1.993, -41.137, 12.313, -1.152, 0.03384553}, 695 { 12, 6.542, 1.00000E+09, -13.389, 1.997, -53. 733 { 12, 6.542, 1.00000E+09, -13.389, 1.997, -53.549, 17.420, -1.840, 0.06431849}, 696 { 13, 6.570, 1.49968E+09, -13.401, 1.997, -66. 734 { 13, 6.570, 1.49968E+09, -13.401, 1.997, -66.243, 22.297, -2.460, 0.09045854}, 697 { 14, 6.364, 1.49968E+09, -13.452, 1.999, -78. 735 { 14, 6.364, 1.49968E+09, -13.452, 1.999, -78.271, 26.757, -3.008, 0.1128195}, 698 { 15, 6.500, 1.49968E+09, -13.488, 1.998, -85. 736 { 15, 6.500, 1.49968E+09, -13.488, 1.998, -85.069, 29.164, -3.291, 0.1239113}, 699 { 16, 6.500, 1.49968E+09, -13.532, 1.998, -93. 737 { 16, 6.500, 1.49968E+09, -13.532, 1.998, -93.640, 32.274, -3.665, 0.1388633}, 700 { 17, 6.500, 1.49968E+09, -13.584, 2.000, -98. 738 { 17, 6.500, 1.49968E+09, -13.584, 2.000, -98.534, 33.958, -3.857, 0.1461557}, 701 { 18, 6.500, 1.49968E+09, -13.618, 1.999, -100 739 { 18, 6.500, 1.49968E+09, -13.618, 1.999, -100.077, 34.379, -3.891, 0.1468902}, 702 { 19, 6.500, 1.99986E+09, -13.185, 1.992, -53. 740 { 19, 6.500, 1.99986E+09, -13.185, 1.992, -53.819, 17.528, -1.851, 0.0648722}, 703 { 20, 6.490, 1.99986E+09, -13.123, 1.993, -52. 741 { 20, 6.490, 1.99986E+09, -13.123, 1.993, -52.221, 17.169, -1.832, 0.06502094}, 704 { 21, 6.498, 1.99986E+09, -13.157, 1.994, -55. 742 { 21, 6.498, 1.99986E+09, -13.157, 1.994, -55.365, 18.276, -1.961, 0.07002778}, 705 { 22, 6.495, 1.99986E+09, -13.183, 1.994, -57. 743 { 22, 6.495, 1.99986E+09, -13.183, 1.994, -57.412, 18.957, -2.036, 0.07278856}, 706 { 23, 6.487, 1.99986E+09, -13.216, 1.995, -58. 744 { 23, 6.487, 1.99986E+09, -13.216, 1.995, -58.478, 19.270, -2.065, 0.07362722}, 707 { 24, 6.500, 1.99986E+09, -13.330, 1.997, -62. 745 { 24, 6.500, 1.99986E+09, -13.330, 1.997, -62.192, 20.358, -2.167, 0.07666583}, 708 { 25, 6.488, 1.99986E+09, -13.277, 1.997, -58. 746 { 25, 6.488, 1.99986E+09, -13.277, 1.997, -58.007, 18.924, -2.003, 0.0704305}, 709 { 26, 6.500, 5.00035E+09, -13.292, 1.997, -61. 747 { 26, 6.500, 5.00035E+09, -13.292, 1.997, -61.176, 20.067, -2.141, 0.0760269}, 710 { 27, 6.500, 5.00035E+09, -13.321, 1.998, -61. 748 { 27, 6.500, 5.00035E+09, -13.321, 1.998, -61.909, 20.271, -2.159, 0.07653559}, 711 { 28, 6.500, 5.00035E+09, -13.340, 1.998, -62. 749 { 28, 6.500, 5.00035E+09, -13.340, 1.998, -62.402, 20.391, -2.167, 0.07664243}, 712 { 29, 6.500, 5.00035E+09, -13.439, 1.998, -67. 750 { 29, 6.500, 5.00035E+09, -13.439, 1.998, -67.305, 21.954, -2.331, 0.0823267}, 713 { 30, 6.500, 5.00035E+09, -13.383, 1.999, -62. 751 { 30, 6.500, 5.00035E+09, -13.383, 1.999, -62.064, 20.136, -2.122, 0.07437589}, 714 { 31, 6.500, 5.00035E+09, -13.349, 1.997, -61. 752 { 31, 6.500, 5.00035E+09, -13.349, 1.997, -61.068, 19.808, -2.086, 0.07307488}, 715 { 32, 6.500, 5.00035E+09, -13.373, 1.999, -63. 753 { 32, 6.500, 5.00035E+09, -13.373, 1.999, -63.126, 20.553, -2.175, 0.07660222}, 716 { 33, 6.500, 5.00035E+09, -13.395, 1.999, -65. 754 { 33, 6.500, 5.00035E+09, -13.395, 1.999, -65.674, 21.445, -2.278, 0.08054694}, 717 { 34, 6.500, 5.00035E+09, -13.417, 1.999, -69. 755 { 34, 6.500, 5.00035E+09, -13.417, 1.999, -69.457, 22.811, -2.442, 0.08709536}, 718 { 35, 6.500, 5.00035E+09, -13.442, 2.000, -72. 756 { 35, 6.500, 5.00035E+09, -13.442, 2.000, -72.283, 23.808, -2.558, 0.09156808}, 719 { 36, 6.500, 5.00035E+09, -13.451, 1.998, -74. 757 { 36, 6.500, 5.00035E+09, -13.451, 1.998, -74.696, 24.641, -2.653, 0.09516597}, 720 { 37, 6.500, 5.00035E+09, -13.082, 1.991, -46. 758 { 37, 6.500, 5.00035E+09, -13.082, 1.991, -46.235, 14.519, -1.458, 0.04837659}, 721 { 38, 6.465, 5.00035E+09, -13.022, 1.993, -41. 759 { 38, 6.465, 5.00035E+09, -13.022, 1.993, -41.784, 13.065, -1.300, 0.04267703}, 722 { 39, 6.492, 5.00035E+09, -13.043, 1.994, -44. 760 { 39, 6.492, 5.00035E+09, -13.043, 1.994, -44.609, 14.114, -1.429, 0.0479348}, 723 { 40, 6.499, 5.00035E+09, -13.064, 1.994, -47. 761 { 40, 6.499, 5.00035E+09, -13.064, 1.994, -47.142, 15.019, -1.536, 0.0521347}, 724 { 41, 6.384, 5.00035E+09, -13.156, 1.996, -53. 762 { 41, 6.384, 5.00035E+09, -13.156, 1.996, -53.114, 17.052, -1.766, 0.06079426}, 725 { 42, 6.500, 5.00035E+09, -13.176, 1.996, -54. 763 { 42, 6.500, 5.00035E+09, -13.176, 1.996, -54.590, 17.550, -1.822, 0.06290335}, 726 { 43, 6.500, 5.00035E+09, -13.133, 1.997, -51. 764 { 43, 6.500, 5.00035E+09, -13.133, 1.997, -51.272, 16.423, -1.694, 0.05806108}, 727 { 44, 6.500, 5.00035E+09, -13.220, 1.996, -58. 765 { 44, 6.500, 5.00035E+09, -13.220, 1.996, -58.314, 18.839, -1.969, 0.0684608}, 728 { 45, 6.500, 5.00035E+09, -13.246, 1.998, -59. 766 { 45, 6.500, 5.00035E+09, -13.246, 1.998, -59.674, 19.295, -2.020, 0.07037294}, 729 { 46, 6.500, 5.00035E+09, -13.407, 1.999, -72. 767 { 46, 6.500, 5.00035E+09, -13.407, 1.999, -72.228, 23.693, -2.532, 0.09017969}, 730 { 47, 6.500, 5.00035E+09, -13.277, 1.998, -60. 768 { 47, 6.500, 5.00035E+09, -13.277, 1.998, -60.890, 19.647, -2.053, 0.07138694}, 731 { 48, 6.500, 5.00035E+09, -13.222, 1.998, -56. 769 { 48, 6.500, 5.00035E+09, -13.222, 1.998, -56.152, 18.002, -1.863, 0.06410123}, 732 { 49, 6.500, 5.00035E+09, -13.199, 1.997, -56. 770 { 49, 6.500, 5.00035E+09, -13.199, 1.997, -56.208, 18.052, -1.872, 0.06456884}, 733 { 50, 6.500, 5.00035E+09, -13.215, 1.998, -58. 771 { 50, 6.500, 5.00035E+09, -13.215, 1.998, -58.478, 18.887, -1.973, 0.06860356}, 734 { 51, 6.500, 5.00035E+09, -13.230, 1.998, -60. 772 { 51, 6.500, 5.00035E+09, -13.230, 1.998, -60.708, 19.676, -2.066, 0.07225841}, 735 { 52, 6.500, 7.99834E+09, -13.246, 1.998, -63. 773 { 52, 6.500, 7.99834E+09, -13.246, 1.998, -63.341, 20.632, -2.180, 0.0767412}, 736 { 53, 6.500, 5.00035E+09, -13.262, 1.998, -66. 774 { 53, 6.500, 5.00035E+09, -13.262, 1.998, -66.339, 21.716, -2.310, 0.08191981}, 737 { 54, 6.500, 7.99834E+09, -13.279, 1.998, -67. 775 { 54, 6.500, 7.99834E+09, -13.279, 1.998, -67.649, 22.151, -2.357, 0.08357825}, 738 { 55, 6.500, 5.00035E+09, -12.951, 1.990, -45. 776 { 55, 6.500, 5.00035E+09, -12.951, 1.990, -45.302, 14.219, -1.423, 0.04712317}, 739 { 56, 6.425, 5.00035E+09, -12.882, 1.992, -39. 777 { 56, 6.425, 5.00035E+09, -12.882, 1.992, -39.825, 12.363, -1.214, 0.03931009}, 740 { 57, 6.466, 2.82488E+09, -12.903, 1.992, -38. 778 { 57, 6.466, 2.82488E+09, -12.903, 1.992, -38.952, 11.982, -1.160, 0.03681554}, 741 { 58, 6.451, 5.00035E+09, -12.915, 1.993, -41. 779 { 58, 6.451, 5.00035E+09, -12.915, 1.993, -41.959, 13.118, -1.302, 0.04271291}, 742 { 59, 6.434, 5.00035E+09, -12.914, 1.993, -40. 780 { 59, 6.434, 5.00035E+09, -12.914, 1.993, -40.528, 12.555, -1.230, 0.03971407}, 743 { 60, 6.444, 5.00035E+09, -12.922, 1.992, -39. 781 { 60, 6.444, 5.00035E+09, -12.922, 1.992, -39.986, 12.329, -1.200, 0.03843737}, 744 { 61, 6.414, 7.99834E+09, -12.930, 1.993, -42. 782 { 61, 6.414, 7.99834E+09, -12.930, 1.993, -42.756, 13.362, -1.327, 0.0436124}, 745 { 62, 6.420, 7.99834E+09, -12.938, 1.992, -42. 783 { 62, 6.420, 7.99834E+09, -12.938, 1.992, -42.682, 13.314, -1.319, 0.04322509}, 746 { 63, 6.416, 7.99834E+09, -12.946, 1.993, -42. 784 { 63, 6.416, 7.99834E+09, -12.946, 1.993, -42.399, 13.185, -1.301, 0.04243861}, 747 { 64, 6.443, 7.99834E+09, -12.963, 1.993, -43. 785 { 64, 6.443, 7.99834E+09, -12.963, 1.993, -43.226, 13.475, -1.335, 0.04377341}, 748 { 65, 6.449, 7.99834E+09, -12.973, 1.993, -43. 786 { 65, 6.449, 7.99834E+09, -12.973, 1.993, -43.232, 13.456, -1.330, 0.04347536}, 749 { 66, 6.419, 7.99834E+09, -12.966, 1.993, -42. 787 { 66, 6.419, 7.99834E+09, -12.966, 1.993, -42.047, 12.990, -1.270, 0.04095499}, 750 { 67, 6.406, 7.99834E+09, -12.976, 1.993, -42. 788 { 67, 6.406, 7.99834E+09, -12.976, 1.993, -42.405, 13.106, -1.283, 0.04146024}, 751 { 68, 6.424, 7.99834E+09, -12.986, 1.993, -41. 789 { 68, 6.424, 7.99834E+09, -12.986, 1.993, -41.974, 12.926, -1.259, 0.040435}, 752 { 69, 6.417, 7.99834E+09, -12.989, 1.993, -42. 790 { 69, 6.417, 7.99834E+09, -12.989, 1.993, -42.132, 12.967, -1.262, 0.04048908}, 753 { 70, 6.405, 7.99834E+09, -13.000, 1.994, -42. 791 { 70, 6.405, 7.99834E+09, -13.000, 1.994, -42.582, 13.122, -1.280, 0.04119599}, 754 { 71, 6.449, 7.99834E+09, -13.015, 1.994, -42. 792 { 71, 6.449, 7.99834E+09, -13.015, 1.994, -42.586, 13.115, -1.278, 0.04107587}, 755 { 72, 6.465, 7.99834E+09, -13.030, 1.994, -43. 793 { 72, 6.465, 7.99834E+09, -13.030, 1.994, -43.708, 13.509, -1.324, 0.04286491}, 756 { 73, 6.447, 7.99834E+09, -13.048, 1.996, -44. 794 { 73, 6.447, 7.99834E+09, -13.048, 1.996, -44.838, 13.902, -1.369, 0.04457132}, 757 { 74, 6.452, 7.99834E+09, -13.073, 1.997, -45. 795 { 74, 6.452, 7.99834E+09, -13.073, 1.997, -45.545, 14.137, -1.395, 0.04553459}, 758 { 75, 6.432, 7.99834E+09, -13.082, 1.997, -46. 796 { 75, 6.432, 7.99834E+09, -13.082, 1.997, -46.426, 14.431, -1.428, 0.04678218}, 759 { 76, 6.439, 7.99834E+09, -13.100, 1.997, -47. 797 { 76, 6.439, 7.99834E+09, -13.100, 1.997, -47.513, 14.806, -1.471, 0.04842566}, 760 { 77, 6.432, 7.99834E+09, -13.110, 1.997, -48. 798 { 77, 6.432, 7.99834E+09, -13.110, 1.997, -48.225, 15.042, -1.497, 0.04938364}, 761 { 78, 6.500, 7.99834E+09, -13.185, 1.997, -53. 799 { 78, 6.500, 7.99834E+09, -13.185, 1.997, -53.256, 16.739, -1.687, 0.05645173}, 762 { 79, 6.500, 7.99834E+09, -13.200, 1.997, -53. 800 { 79, 6.500, 7.99834E+09, -13.200, 1.997, -53.900, 16.946, -1.709, 0.05723134}, 763 { 80, 6.500, 7.99834E+09, -13.156, 1.998, -49. 801 { 80, 6.500, 7.99834E+09, -13.156, 1.998, -49.801, 15.536, -1.547, 0.05103522}, 764 { 81, 6.500, 7.99834E+09, -13.128, 1.997, -49. 802 { 81, 6.500, 7.99834E+09, -13.128, 1.997, -49.651, 15.512, -1.548, 0.05123203}, 765 { 82, 6.500, 7.99834E+09, -13.134, 1.997, -51. 803 { 82, 6.500, 7.99834E+09, -13.134, 1.997, -51.021, 16.018, -1.609, 0.05364831}, 766 { 83, 6.500, 7.99834E+09, -13.148, 1.998, -52. 804 { 83, 6.500, 7.99834E+09, -13.148, 1.998, -52.693, 16.612, -1.679, 0.05638698}, 767 { 84, 6.500, 7.99834E+09, -13.161, 1.998, -54. 805 { 84, 6.500, 7.99834E+09, -13.161, 1.998, -54.415, 17.238, -1.754, 0.05935566}, 768 { 85, 6.500, 7.99834E+09, -13.175, 1.998, -56. 806 { 85, 6.500, 7.99834E+09, -13.175, 1.998, -56.083, 17.834, -1.824, 0.06206068}, 769 { 86, 6.500, 7.99834E+09, -13.189, 1.998, -57. 807 { 86, 6.500, 7.99834E+09, -13.189, 1.998, -57.860, 18.463, -1.898, 0.0649633}, 770 { 87, 6.500, 7.99834E+09, -12.885, 1.990, -39. 808 { 87, 6.500, 7.99834E+09, -12.885, 1.990, -39.973, 12.164, -1.162, 0.0364598}, 771 { 88, 6.417, 7.99834E+09, -12.816, 1.991, -34. 809 { 88, 6.417, 7.99834E+09, -12.816, 1.991, -34.591, 10.338, -0.956, 0.0287409}, 772 { 89, 6.442, 7.99834E+09, -12.831, 1.992, -36. 810 { 89, 6.442, 7.99834E+09, -12.831, 1.992, -36.002, 10.867, -1.021, 0.03136835}, 773 { 90, 6.463, 7.99834E+09, -12.850, 1.993, -37. 811 { 90, 6.463, 7.99834E+09, -12.850, 1.993, -37.660, 11.475, -1.095, 0.03435334}, 774 { 91, 6.447, 7.99834E+09, -12.852, 1.993, -37. 812 { 91, 6.447, 7.99834E+09, -12.852, 1.993, -37.268, 11.301, -1.071, 0.0330539}, 775 { 92, 6.439, 7.99834E+09, -12.858, 1.993, -37. 813 { 92, 6.439, 7.99834E+09, -12.858, 1.993, -37.695, 11.438, -1.085, 0.03376669}, 776 { 93, 6.437, 1.00000E+10, -12.866, 1.993, -39. 814 { 93, 6.437, 1.00000E+10, -12.866, 1.993, -39.010, 11.927, -1.146, 0.03630848}, 777 { 94, 6.432, 7.99834E+09, -12.862, 1.993, -37. 815 { 94, 6.432, 7.99834E+09, -12.862, 1.993, -37.192, 11.229, -1.057, 0.0325621}, 778 { 95, 6.435, 7.99834E+09, -12.869, 1.993, -37. 816 { 95, 6.435, 7.99834E+09, -12.869, 1.993, -37.589, 11.363, -1.072, 0.03312393}, 779 { 96, 6.449, 1.00000E+10, -12.886, 1.993, -39. 817 { 96, 6.449, 1.00000E+10, -12.886, 1.993, -39.573, 12.095, -1.162, 0.03680527}, 780 { 97, 6.446, 1.00000E+10, -12.892, 1.993, -40. 818 { 97, 6.446, 1.00000E+10, -12.892, 1.993, -40.007, 12.242, -1.178, 0.03737377}, 781 { 98, 6.421, 1.00000E+10, -12.887, 1.993, -39. 819 { 98, 6.421, 1.00000E+10, -12.887, 1.993, -39.509, 12.041, -1.152, 0.03629023}, 782 { 99, 6.414, 1.00000E+10, -12.894, 1.993, -39. 820 { 99, 6.414, 1.00000E+10, -12.894, 1.993, -39.939, 12.183, -1.168, 0.03690464}, 783 {100, 6.412, 1.00000E+10, -12.900, 1.993, -39. 821 {100, 6.412, 1.00000E+10, -12.900, 1.993, -39.973, 12.180, -1.166, 0.036773} 784 }; 822 }; 785 823 786 824 787 825 788 826