Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // neutron_hp -- source file 27 // J.P. Wellisch, Nov-1996 28 // A prototype of the low energy neutron transport model. 29 // 30 // 070523 add neglecting doppler broadening on the fly. T. Koi 31 // 070613 fix memory leaking by T. Koi 32 // 071002 enable cross section dump by T. Koi 33 // 080428 change checking point of "neglecting doppler broadening" flag 34 // from GetCrossSection to BuildPhysicsTable by T. Koi 35 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties:: 36 // 37 // P. Arce, June-2014 Conversion neutron_hp to particle_hp 38 // 39 #include "G4ParticleHPElasticData.hh" 40 41 #include "G4Electron.hh" 42 #include "G4ElementTable.hh" 43 #include "G4HadronicParameters.hh" 44 #include "G4Neutron.hh" 45 #include "G4NucleiProperties.hh" 46 #include "G4Nucleus.hh" 47 #include "G4ParticleHPData.hh" 48 #include "G4ParticleHPManager.hh" 49 #include "G4PhysicalConstants.hh" 50 #include "G4Pow.hh" 51 #include "G4SystemOfUnits.hh" 52 53 G4ParticleHPElasticData::G4ParticleHPElasticData() : G4VCrossSectionDataSet("NeutronHPElasticXS") 54 { 55 SetMinKinEnergy(0 * MeV); 56 SetMaxKinEnergy(20 * MeV); 57 58 theCrossSections = nullptr; 59 instanceOfWorker = false; 60 if (G4Threading::IsWorkerThread()) { 61 instanceOfWorker = true; 62 } 63 element_cache = nullptr; 64 material_cache = nullptr; 65 ke_cache = 0.0; 66 xs_cache = 0.0; 67 // BuildPhysicsTable( *G4Neutron::Neutron() ); 68 } 69 70 G4ParticleHPElasticData::~G4ParticleHPElasticData() 71 { 72 if (theCrossSections != nullptr && !instanceOfWorker) { 73 theCrossSections->clearAndDestroy(); 74 delete theCrossSections; 75 theCrossSections = nullptr; 76 } 77 } 78 79 G4bool G4ParticleHPElasticData::IsIsoApplicable(const G4DynamicParticle* dp, G4int /*Z*/, 80 G4int /*A*/, const G4Element* /*elm*/, 81 const G4Material* /*mat*/) 82 { 83 G4double eKin = dp->GetKineticEnergy(); 84 return eKin <= GetMaxKinEnergy() && eKin >= GetMinKinEnergy() 85 && dp->GetDefinition() == G4Neutron::Neutron(); 86 } 87 88 G4double G4ParticleHPElasticData::GetIsoCrossSection(const G4DynamicParticle* dp, G4int /*Z*/, 89 G4int /*A*/, const G4Isotope* /*iso*/, 90 const G4Element* element, 91 const G4Material* material) 92 { 93 if (dp->GetKineticEnergy() == ke_cache && element == element_cache && material == material_cache) 94 return xs_cache; 95 96 ke_cache = dp->GetKineticEnergy(); 97 element_cache = element; 98 material_cache = material; 99 G4double xs = GetCrossSection(dp, element, material->GetTemperature()); 100 xs_cache = xs; 101 return xs; 102 } 103 104 void G4ParticleHPElasticData::BuildPhysicsTable(const G4ParticleDefinition& aP) 105 { 106 if (&aP != G4Neutron::Neutron()) 107 throw G4HadronicException(__FILE__, __LINE__, 108 "Attempt to use NeutronHP data for particles other than neutrons!!!"); 109 110 if (G4Threading::IsWorkerThread()) { 111 theCrossSections = G4ParticleHPManager::GetInstance()->GetElasticCrossSections(); 112 return; 113 } 114 115 std::size_t numberOfElements = G4Element::GetNumberOfElements(); 116 // TKDB 117 // if ( theCrossSections == 0 ) theCrossSections = new G4PhysicsTable( numberOfElements ); 118 if (theCrossSections == nullptr) 119 theCrossSections = new G4PhysicsTable(numberOfElements); 120 else 121 theCrossSections->clearAndDestroy(); 122 123 // make a PhysicsVector for each element 124 125 auto theElementTable = G4Element::GetElementTable(); 126 for (std::size_t i = 0; i < numberOfElements; ++i) { 127 G4PhysicsVector* physVec = G4ParticleHPData::Instance(G4Neutron::Neutron()) 128 ->MakePhysicsVector((*theElementTable)[i], this); 129 theCrossSections->push_back(physVec); 130 } 131 132 G4ParticleHPManager::GetInstance()->RegisterElasticCrossSections(theCrossSections); 133 } 134 135 void G4ParticleHPElasticData::DumpPhysicsTable(const G4ParticleDefinition&) 136 { 137 #ifdef G4VERBOSE 138 if (G4HadronicParameters::Instance()->GetVerboseLevel() == 0) return; 139 140 // 141 // Dump element based cross section 142 // range 10e-5 eV to 20 MeV 143 // 10 point per decade 144 // in barn 145 // 146 147 G4cout << G4endl; 148 G4cout << G4endl; 149 G4cout << "Elastic Cross Section of Neutron HP" << G4endl; 150 G4cout << "(Pointwise cross-section at 0 Kelvin.)" << G4endl; 151 G4cout << G4endl; 152 G4cout << "Name of Element" << G4endl; 153 G4cout << "Energy[eV] XS[barn]" << G4endl; 154 G4cout << G4endl; 155 156 std::size_t numberOfElements = G4Element::GetNumberOfElements(); 157 auto theElementTable = G4Element::GetElementTable(); 158 159 for (std::size_t i = 0; i < numberOfElements; ++i) { 160 G4cout << (*theElementTable)[i]->GetName() << G4endl; 161 G4int ie = 0; 162 163 for (ie = 0; ie < 130; ++ie) { 164 G4double eKinetic = 1.0e-5 * G4Pow::GetInstance()->powA(10.0, ie / 10.0) * eV; 165 G4bool outOfRange = false; 166 167 if (eKinetic < 20 * MeV) { 168 G4cout << eKinetic / eV << " " 169 << (*((*theCrossSections)(i))).GetValue(eKinetic, outOfRange) / barn << G4endl; 170 } 171 } 172 G4cout << G4endl; 173 } 174 #endif 175 } 176 177 G4double G4ParticleHPElasticData::GetCrossSection(const G4DynamicParticle* aP, const G4Element* anE, 178 G4double aT) 179 { 180 G4double result = 0; 181 G4bool outOfRange; 182 auto index = (G4int)anE->GetIndex(); 183 184 // prepare neutron 185 G4double eKinetic = aP->GetKineticEnergy(); 186 187 if (G4ParticleHPManager::GetInstance()->GetNeglectDoppler()) { 188 // NEGLECT_DOPPLER_B. 189 G4double factor = 1.0; 190 if (eKinetic < aT * k_Boltzmann) { 191 // below 0.1 eV neutrons 192 // Have to do some, but now just igonre. 193 // Will take care after performance check. 194 // factor = factor * targetV; 195 } 196 return ((*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange)) * factor; 197 } 198 199 G4ReactionProduct theNeutron(aP->GetDefinition()); 200 theNeutron.SetMomentum(aP->GetMomentum()); 201 theNeutron.SetKineticEnergy(eKinetic); 202 203 // prepare thermal nucleus 204 G4Nucleus aNuc; 205 G4double eps = 0.0001; 206 G4double theA = anE->GetN(); 207 G4double theZ = anE->GetZ(); 208 G4double eleMass; 209 210 eleMass = (G4NucleiProperties::GetNuclearMass(G4int(theA + eps), G4int(theZ + eps))) 211 / G4Neutron::Neutron()->GetPDGMass(); 212 213 G4ReactionProduct boosted; 214 G4double aXsection; 215 216 // MC integration loop 217 G4int counter = 0; 218 G4double buffer = 0; 219 G4int size = G4int(std::max(10., aT / 60 * kelvin)); 220 G4ThreeVector neutronVelocity = 221 1. / G4Neutron::Neutron()->GetPDGMass() * theNeutron.GetMomentum(); 222 G4double neutronVMag = neutronVelocity.mag(); 223 224 while (counter == 0 225 || std::abs(buffer - result / std::max(1, counter)) 226 > 0.03 * buffer) // Loop checking, 11.05.2015, T. Koi 227 { 228 if (counter != 0) buffer = result / counter; 229 while (counter < size) // Loop checking, 11.05.2015, T. Koi 230 { 231 counter++; 232 G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus(eleMass, aT); 233 boosted.Lorentz(theNeutron, aThermalNuc); 234 G4double theEkin = boosted.GetKineticEnergy(); 235 aXsection = (*((*theCrossSections)(index))).GetValue(theEkin, outOfRange); 236 // velocity correction. 237 G4ThreeVector targetVelocity = 1. / aThermalNuc.GetMass() * aThermalNuc.GetMomentum(); 238 aXsection *= (targetVelocity - neutronVelocity).mag() / neutronVMag; 239 result += aXsection; 240 } 241 size += size; 242 } 243 result /= counter; 244 /* 245 // Checking impact of G4NEUTRONHP_NEGLECT_DOPPLER 246 G4cout << " result " << result << " " 247 << (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) << " " 248 << (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) /result << G4endl; 249 */ 250 return result; 251 } 252 253 G4int G4ParticleHPElasticData::GetVerboseLevel() const 254 { 255 return G4ParticleHPManager::GetInstance()->GetVerboseLevel(); 256 } 257 258 void G4ParticleHPElasticData::SetVerboseLevel(G4int newValue) 259 { 260 G4ParticleHPManager::GetInstance()->SetVerboseLevel(newValue); 261 } 262 263 void G4ParticleHPElasticData::CrossSectionDescription(std::ostream& outFile) const 264 { 265 outFile << "High Precision cross data based on Evaluated Nuclear Data Files (ENDF) for elastic " 266 "reaction of neutrons below 20MeV\n"; 267 } 268