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 // $Id: G4NeutronElasticXS.cc 76889 2013-11-18 13:01:55Z gcosmo $ >> 27 // 26 // ------------------------------------------- 28 // ------------------------------------------------------------------- 27 // 29 // 28 // GEANT4 Class file 30 // GEANT4 Class file 29 // 31 // 30 // 32 // 31 // File name: G4NeutronElasticXS 33 // File name: G4NeutronElasticXS 32 // 34 // 33 // Author Ivantchenko, Geant4, 3-Aug-09 35 // Author Ivantchenko, Geant4, 3-Aug-09 34 // 36 // 35 // Modifications: 37 // Modifications: 36 // 38 // 37 39 38 #include "G4NeutronElasticXS.hh" 40 #include "G4NeutronElasticXS.hh" 39 #include "G4Neutron.hh" 41 #include "G4Neutron.hh" 40 #include "G4DynamicParticle.hh" 42 #include "G4DynamicParticle.hh" 41 #include "G4ElementTable.hh" << 42 #include "G4Material.hh" << 43 #include "G4Element.hh" 43 #include "G4Element.hh" >> 44 #include "G4ElementTable.hh" 44 #include "G4PhysicsLogVector.hh" 45 #include "G4PhysicsLogVector.hh" 45 #include "G4CrossSectionDataSetRegistry.hh" << 46 #include "G4PhysicsVector.hh" 46 #include "G4ComponentGGHadronNucleusXsc.hh" << 47 #include "G4GlauberGribovCrossSection.hh" 47 #include "G4HadronicParameters.hh" << 48 #include "G4HadronNucleonXsc.hh" 48 #include "Randomize.hh" << 49 #include "G4NistManager.hh" 49 #include "G4SystemOfUnits.hh" << 50 #include "G4Proton.hh" 50 #include "G4IsotopeList.hh" << 51 #include "G4AutoLock.hh" << 52 51 >> 52 #include <iostream> 53 #include <fstream> 53 #include <fstream> 54 #include <sstream> 54 #include <sstream> 55 55 56 G4PhysicsVector* G4NeutronElasticXS::data[] = << 56 using namespace std; 57 G4double G4NeutronElasticXS::coeff[] = {0.0}; << 58 G4String G4NeutronElasticXS::gDataDirectory = << 59 G4bool G4NeutronElasticXS::fLock = true; << 60 << 61 namespace << 62 { << 63 G4Mutex nElasticXSMutex = G4MUTEX_INITIALIZE << 64 } << 65 57 66 G4NeutronElasticXS::G4NeutronElasticXS() 58 G4NeutronElasticXS::G4NeutronElasticXS() 67 : G4VCrossSectionDataSet(Default_Name()), << 59 : G4VCrossSectionDataSet("G4NeutronElasticXS"), 68 neutron(G4Neutron::Neutron()) << 60 proton(G4Proton::Proton()), maxZ(92) 69 { 61 { 70 // verboseLevel = 0; 62 // verboseLevel = 0; 71 if (verboseLevel > 0){ << 63 if(verboseLevel > 0){ 72 G4cout << "G4NeutronElasticXS::G4NeutronE 64 G4cout << "G4NeutronElasticXS::G4NeutronElasticXS Initialise for Z < " 73 << MAXZEL << G4endl; << 65 << maxZ + 1 << G4endl; 74 } 66 } 75 ggXsection = << 67 data.resize(maxZ+1, 0); 76 G4CrossSectionDataSetRegistry::Instance()- << 68 coeff.resize(maxZ+1, 1.0); 77 if (ggXsection == nullptr) << 69 ggXsection = new G4GlauberGribovCrossSection(); 78 ggXsection = new G4ComponentGGHadronNucleu << 70 fNucleon = new G4HadronNucleonXsc(); 79 SetForAllAtomsAndEnergies(true); << 71 isInitialized = false; 80 FindDirectoryPath(); << 81 } 72 } 82 73 83 G4NeutronElasticXS::~G4NeutronElasticXS() 74 G4NeutronElasticXS::~G4NeutronElasticXS() 84 { 75 { 85 if (isFirst) { << 76 delete fNucleon; 86 for(G4int i=0; i<MAXZEL; ++i) { << 77 /* 87 delete data[i]; << 78 for(G4int i=0; i<=maxZ; ++i) { 88 data[i] = nullptr; << 79 delete data[i]; 89 } << 90 } 80 } >> 81 */ 91 } 82 } 92 83 93 void G4NeutronElasticXS::CrossSectionDescripti 84 void G4NeutronElasticXS::CrossSectionDescription(std::ostream& outFile) const 94 { 85 { 95 outFile << "G4NeutronElasticXS calculates th 86 outFile << "G4NeutronElasticXS calculates the neutron elastic scattering\n" 96 << "cross section on nuclei using da 87 << "cross section on nuclei using data from the high precision\n" 97 << "neutron database. These data ar 88 << "neutron database. These data are simplified and smoothed over\n" 98 << "the resonance region in order to 89 << "the resonance region in order to reduce CPU time.\n" 99 << "For high energies Glauber-Gribiv << 90 << "G4NeutronElasticXS is valid for energies up to 20 MeV, for all\n" >> 91 << "targets through U.\n"; 100 } 92 } 101 93 102 G4bool 94 G4bool 103 G4NeutronElasticXS::IsElementApplicable(const 95 G4NeutronElasticXS::IsElementApplicable(const G4DynamicParticle*, 104 G4int, const G4Material*) 96 G4int, const G4Material*) 105 { 97 { 106 return true; 98 return true; 107 } 99 } 108 100 109 G4bool G4NeutronElasticXS::IsIsoApplicable(con << 110 G4i << 111 con << 112 { << 113 return false; << 114 } << 115 << 116 G4double 101 G4double 117 G4NeutronElasticXS::GetElementCrossSection(con 102 G4NeutronElasticXS::GetElementCrossSection(const G4DynamicParticle* aParticle, 118 G4int Z, const G4Material*) 103 G4int Z, const G4Material*) 119 { 104 { 120 return ElementCrossSection(aParticle->GetKin << 105 G4double xs = 0.0; 121 aParticle->GetLogKineticEnergy(), Z << 106 G4double ekin = aParticle->GetKineticEnergy(); 122 } << 123 107 124 G4double << 108 if(Z < 1 || Z > maxZ) { return xs; } 125 G4NeutronElasticXS::ComputeCrossSectionPerElem << 126 const G4ParticleDefinition*, << 127 const G4Element* elm, << 128 const G4Material*) << 129 { << 130 return ElementCrossSection(ekin, loge, elm-> << 131 } << 132 << 133 G4double G4NeutronElasticXS::ElementCrossSecti << 134 { << 135 G4int Z = (ZZ >= MAXZEL) ? MAXZEL - 1 : ZZ; << 136 auto pv = GetPhysicsVector(Z); << 137 109 138 G4double xs = (ekin <= pv->GetMaxEnergy()) ? << 110 G4int Amean = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z)); 139 : coeff[Z]*ggXsection->GetElasticElementCr << 111 G4PhysicsVector* pv = data[Z]; 140 << 112 // G4cout << "G4NeutronElasticXS::GetCrossSection e= " << ekin 141 << 113 // << " Z= " << Z << G4endl; 142 #ifdef G4VERBOSE << 114 143 if(verboseLevel > 1) { << 115 // element was not initialised 144 G4cout << "Z= " << Z << " Ekin(MeV)= " << << 116 if(!pv) { 145 << ", nElmXSel(b)= " << xs/CLHEP::barn << 117 Initialise(Z); 146 << G4endl; << 118 pv = data[Z]; >> 119 if(!pv) { return xs; } >> 120 } >> 121 >> 122 G4double e1 = pv->Energy(0); >> 123 if(ekin <= e1) { return (*pv)[0]; } >> 124 >> 125 G4int n = pv->GetVectorLength() - 1; >> 126 G4double e2 = pv->Energy(n); >> 127 >> 128 if(ekin <= e2) { >> 129 xs = pv->Value(ekin); >> 130 } else if(1 == Z) { >> 131 fNucleon->GetHadronNucleonXscPDG(aParticle, proton); >> 132 xs = coeff[1]*fNucleon->GetElasticHadronNucleonXsc(); >> 133 } else { >> 134 ggXsection->GetIsoCrossSection(aParticle, Z, Amean); >> 135 xs = coeff[Z]*ggXsection->GetElasticGlauberGribovXsc(); 147 } 136 } 148 #endif << 149 return xs; << 150 } << 151 << 152 G4double << 153 G4NeutronElasticXS::ComputeIsoCrossSection(G4d << 154 const G4ParticleDefinition* << 155 G4int Z, G4int A, << 156 const G4Isotope*, const G4E << 157 const G4Material*) << 158 { << 159 return ElementCrossSection(ekin, loge, Z)*A/ << 160 } << 161 << 162 G4double << 163 G4NeutronElasticXS::GetIsoCrossSection(const G << 164 G4int Z, G4int A, << 165 const G4Isotope*, const G4Eleme << 166 const G4Material*) << 167 { << 168 return ElementCrossSection(aParticle->GetKin << 169 aParticle->GetLogKineticEnergy(), Z << 170 << 171 } << 172 137 173 const G4Isotope* G4NeutronElasticXS::SelectIso << 138 if(verboseLevel > 0){ 174 const G4Element* anElement, G4double, G4 << 139 G4cout << "ekin= " << ekin << ", XSinel= " << xs << G4endl; 175 { << 176 G4int nIso = (G4int)anElement->GetNumberOfIs << 177 const G4Isotope* iso = anElement->GetIsotope << 178 << 179 //G4cout << "SelectIsotope NIso= " << nIso < << 180 if(1 == nIso) { return iso; } << 181 << 182 const G4double* abundVector = anElement->Get << 183 G4double q = G4UniformRand(); << 184 G4double sum = 0.0; << 185 << 186 // isotope wise cross section not used << 187 for (G4int j=0; j<nIso; ++j) { << 188 sum += abundVector[j]; << 189 if(q <= sum) { << 190 iso = anElement->GetIsotope(j); << 191 break; << 192 } << 193 } 140 } 194 return iso; << 141 return xs; 195 } 142 } 196 143 197 void 144 void 198 G4NeutronElasticXS::BuildPhysicsTable(const G4 145 G4NeutronElasticXS::BuildPhysicsTable(const G4ParticleDefinition& p) 199 { 146 { >> 147 if(isInitialized) { return; } 200 if(verboseLevel > 0){ 148 if(verboseLevel > 0){ 201 G4cout << "G4NeutronElasticXS::BuildPhysic 149 G4cout << "G4NeutronElasticXS::BuildPhysicsTable for " 202 << p.GetParticleName() << G4endl; 150 << p.GetParticleName() << G4endl; 203 } 151 } 204 if(p.GetParticleName() != "neutron") { 152 if(p.GetParticleName() != "neutron") { 205 G4ExceptionDescription ed; 153 G4ExceptionDescription ed; 206 ed << p.GetParticleName() << " is a wrong 154 ed << p.GetParticleName() << " is a wrong particle type -" 207 << " only neutron is allowed"; 155 << " only neutron is allowed"; 208 G4Exception("G4NeutronElasticXS::BuildPhys 156 G4Exception("G4NeutronElasticXS::BuildPhysicsTable(..)","had012", 209 FatalException, ed, ""); 157 FatalException, ed, ""); 210 return; 158 return; 211 } 159 } 212 if (fLock || isFirst) { << 160 isInitialized = true; 213 G4AutoLock l(&nElasticXSMutex); << 214 if (fLock) { << 215 isFirst = true; << 216 fLock = false; << 217 FindDirectoryPath(); << 218 } << 219 161 220 // Access to elements << 162 // check environment variable 221 const G4ElementTable* table = G4Element::G << 163 // Build the complete string identifying the file with the data set 222 for ( auto & elm : *table ) { << 164 char* path = getenv("G4NEUTRONXSDATA"); 223 G4int Z = std::max( 1, std::min( elm->Ge << 165 224 if ( nullptr == data[Z] ) { Initialise(Z << 166 G4DynamicParticle* dynParticle = >> 167 new G4DynamicParticle(G4Neutron::Neutron(),G4ThreeVector(1,0,0),1); >> 168 >> 169 // Access to elements >> 170 const G4ElementTable* theElmTable = G4Element::GetElementTable(); >> 171 size_t numOfElm = G4Element::GetNumberOfElements(); >> 172 if(numOfElm > 0) { >> 173 for(size_t i=0; i<numOfElm; ++i) { >> 174 G4int Z = G4int(((*theElmTable)[i])->GetZ()); >> 175 if(Z < 1) { Z = 1; } >> 176 else if(Z > maxZ) { Z = maxZ; } >> 177 //G4cout << "Z= " << Z << G4endl; >> 178 // Initialisation >> 179 if(!data[Z]) { Initialise(Z, dynParticle, path); } 225 } 180 } 226 l.unlock(); << 227 } 181 } >> 182 delete dynParticle; 228 } 183 } 229 184 230 const G4String& G4NeutronElasticXS::FindDirect << 185 void >> 186 G4NeutronElasticXS::Initialise(G4int Z, G4DynamicParticle* dp, >> 187 const char* p) 231 { 188 { 232 // build the complete string identifying the << 189 if(data[Z]) { return; } 233 if (gDataDirectory.empty()) { << 190 const char* path = p; 234 std::ostringstream ost; << 191 if(!p) { 235 ost << G4HadronicParameters::Instance()->G << 192 // check environment variable 236 gDataDirectory = ost.str(); << 193 // Build the complete string identifying the file with the data set >> 194 path = getenv("G4NEUTRONXSDATA"); >> 195 if (!path) { >> 196 G4Exception("G4NeutronElasticXS::Initialise(..)","had013", >> 197 FatalException, >> 198 "Environment variable G4NEUTRONXSDATA is not defined"); >> 199 return; >> 200 } >> 201 } >> 202 G4DynamicParticle* dynParticle = dp; >> 203 if(!dp) { >> 204 dynParticle = >> 205 new G4DynamicParticle(G4Neutron::Neutron(),G4ThreeVector(1,0,0),1); 237 } 206 } 238 return gDataDirectory; << 239 } << 240 << 241 void G4NeutronElasticXS::InitialiseOnFly(G4int << 242 { << 243 G4AutoLock l(&nElasticXSMutex); << 244 Initialise(Z); << 245 l.unlock(); << 246 } << 247 207 248 void G4NeutronElasticXS::Initialise(G4int Z) << 208 G4int Amean = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z)); 249 { << 250 if(data[Z] != nullptr) { return; } << 251 209 252 // upload data from file 210 // upload data from file 253 data[Z] = new G4PhysicsLogVector(); 211 data[Z] = new G4PhysicsLogVector(); 254 212 255 std::ostringstream ost; 213 std::ostringstream ost; 256 ost << FindDirectoryPath() << Z ; << 214 ost << path << "/elast" << Z ; 257 std::ifstream filein(ost.str().c_str()); 215 std::ifstream filein(ost.str().c_str()); 258 if (!filein.is_open()) { << 216 if (!(filein)) { 259 G4ExceptionDescription ed; 217 G4ExceptionDescription ed; 260 ed << "Data file <" << ost.str().c_str() 218 ed << "Data file <" << ost.str().c_str() 261 << "> is not opened!"; 219 << "> is not opened!"; 262 G4Exception("G4NeutronElasticXS::Initialis 220 G4Exception("G4NeutronElasticXS::Initialise(..)","had014", 263 FatalException, ed, "Check G4P << 221 FatalException, ed, "Check G4NEUTRONXSDATA"); 264 return; 222 return; 265 } << 223 }else{ 266 if(verboseLevel > 1) { << 224 if(verboseLevel > 1) { 267 G4cout << "file " << ost.str() << 225 G4cout << "file " << ost.str() 268 << " is opened by G4NeutronElasticXS" << << 226 << " is opened by G4NeutronElasticXS" << G4endl; 269 } << 227 } 270 228 271 // retrieve data from DB << 229 // retrieve data from DB 272 if(!data[Z]->Retrieve(filein, true)) { << 230 if(!data[Z]->Retrieve(filein, true)) { 273 G4ExceptionDescription ed; << 231 G4ExceptionDescription ed; 274 ed << "Data file <" << ost.str().c_str() << 232 ed << "Data file <" << ost.str().c_str() 275 << "> is not retrieved!"; << 233 << "> is not retrieved!"; 276 G4Exception("G4NeutronElasticXS::Initialis << 234 G4Exception("G4NeutronElasticXS::Initialise(..)","had015", 277 FatalException, ed, "Check G4PARTICLEXSDAT << 235 FatalException, ed, "Check G4NEUTRONXSDATA"); 278 return; << 236 return; 279 } << 237 } 280 // smooth transition << 238 281 G4double sig1 = (*(data[Z]))[data[Z]->GetVe << 239 // smooth transition 282 G4double ehigh = data[Z]->GetMaxEnergy(); << 240 size_t n = data[Z]->GetVectorLength() - 1; 283 G4double sig2 = ggXsection->GetElasticEleme << 241 G4double emax = data[Z]->Energy(n); 284 ehigh, Z, aeff[ << 242 G4double sig1 = (*data[Z])[n]; 285 coeff[Z] = (sig2 > 0.) ? sig1/sig2 : 1.0; << 243 dynParticle->SetKineticEnergy(emax); >> 244 G4double sig2 = 0.0; >> 245 if(1 == Z) { >> 246 fNucleon->GetHadronNucleonXscPDG(dynParticle, proton); >> 247 sig2 = fNucleon->GetElasticHadronNucleonXsc(); >> 248 } else { >> 249 ggXsection->GetIsoCrossSection(dynParticle, Z, Amean); >> 250 sig2 = ggXsection->GetElasticGlauberGribovXsc(); >> 251 } >> 252 if(sig2 > 0.) { coeff[Z] = sig1/sig2; } >> 253 } >> 254 if(!dp) { delete dynParticle; } 286 } 255 } 287 256