Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // ------------------------------------------- 27 // 28 // GEANT4 Class file 29 // 30 // 31 // File name: G4NeutronElasticXS 32 // 33 // Author Ivantchenko, Geant4, 3-Aug-09 34 // 35 // Modifications: 36 // 37 38 #include "G4NeutronElasticXS.hh" 39 #include "G4Neutron.hh" 40 #include "G4DynamicParticle.hh" 41 #include "G4ElementTable.hh" 42 #include "G4Material.hh" 43 #include "G4Element.hh" 44 #include "G4PhysicsLogVector.hh" 45 #include "G4CrossSectionDataSetRegistry.hh" 46 #include "G4ComponentGGHadronNucleusXsc.hh" 47 #include "G4HadronicParameters.hh" 48 #include "Randomize.hh" 49 #include "G4SystemOfUnits.hh" 50 #include "G4IsotopeList.hh" 51 #include "G4AutoLock.hh" 52 53 #include <fstream> 54 #include <sstream> 55 56 G4PhysicsVector* G4NeutronElasticXS::data[] = 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 66 G4NeutronElasticXS::G4NeutronElasticXS() 67 : G4VCrossSectionDataSet(Default_Name()), 68 neutron(G4Neutron::Neutron()) 69 { 70 // verboseLevel = 0; 71 if (verboseLevel > 0){ 72 G4cout << "G4NeutronElasticXS::G4NeutronE 73 << MAXZEL << G4endl; 74 } 75 ggXsection = 76 G4CrossSectionDataSetRegistry::Instance()- 77 if (ggXsection == nullptr) 78 ggXsection = new G4ComponentGGHadronNucleu 79 SetForAllAtomsAndEnergies(true); 80 FindDirectoryPath(); 81 } 82 83 G4NeutronElasticXS::~G4NeutronElasticXS() 84 { 85 if (isFirst) { 86 for(G4int i=0; i<MAXZEL; ++i) { 87 delete data[i]; 88 data[i] = nullptr; 89 } 90 } 91 } 92 93 void G4NeutronElasticXS::CrossSectionDescripti 94 { 95 outFile << "G4NeutronElasticXS calculates th 96 << "cross section on nuclei using da 97 << "neutron database. These data ar 98 << "the resonance region in order to 99 << "For high energies Glauber-Gribiv 100 } 101 102 G4bool 103 G4NeutronElasticXS::IsElementApplicable(const 104 G4int, const G4Material*) 105 { 106 return true; 107 } 108 109 G4bool G4NeutronElasticXS::IsIsoApplicable(con 110 G4i 111 con 112 { 113 return false; 114 } 115 116 G4double 117 G4NeutronElasticXS::GetElementCrossSection(con 118 G4int Z, const G4Material*) 119 { 120 return ElementCrossSection(aParticle->GetKin 121 aParticle->GetLogKineticEnergy(), Z 122 } 123 124 G4double 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 138 G4double xs = (ekin <= pv->GetMaxEnergy()) ? 139 : coeff[Z]*ggXsection->GetElasticElementCr 140 141 142 #ifdef G4VERBOSE 143 if(verboseLevel > 1) { 144 G4cout << "Z= " << Z << " Ekin(MeV)= " << 145 << ", nElmXSel(b)= " << xs/CLHEP::barn 146 << G4endl; 147 } 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 173 const G4Isotope* G4NeutronElasticXS::SelectIso 174 const G4Element* anElement, G4double, G4 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 } 194 return iso; 195 } 196 197 void 198 G4NeutronElasticXS::BuildPhysicsTable(const G4 199 { 200 if(verboseLevel > 0){ 201 G4cout << "G4NeutronElasticXS::BuildPhysic 202 << p.GetParticleName() << G4endl; 203 } 204 if(p.GetParticleName() != "neutron") { 205 G4ExceptionDescription ed; 206 ed << p.GetParticleName() << " is a wrong 207 << " only neutron is allowed"; 208 G4Exception("G4NeutronElasticXS::BuildPhys 209 FatalException, ed, ""); 210 return; 211 } 212 if (fLock || isFirst) { 213 G4AutoLock l(&nElasticXSMutex); 214 if (fLock) { 215 isFirst = true; 216 fLock = false; 217 FindDirectoryPath(); 218 } 219 220 // Access to elements 221 const G4ElementTable* table = G4Element::G 222 for ( auto & elm : *table ) { 223 G4int Z = std::max( 1, std::min( elm->Ge 224 if ( nullptr == data[Z] ) { Initialise(Z 225 } 226 l.unlock(); 227 } 228 } 229 230 const G4String& G4NeutronElasticXS::FindDirect 231 { 232 // build the complete string identifying the 233 if (gDataDirectory.empty()) { 234 std::ostringstream ost; 235 ost << G4HadronicParameters::Instance()->G 236 gDataDirectory = ost.str(); 237 } 238 return gDataDirectory; 239 } 240 241 void G4NeutronElasticXS::InitialiseOnFly(G4int 242 { 243 G4AutoLock l(&nElasticXSMutex); 244 Initialise(Z); 245 l.unlock(); 246 } 247 248 void G4NeutronElasticXS::Initialise(G4int Z) 249 { 250 if(data[Z] != nullptr) { return; } 251 252 // upload data from file 253 data[Z] = new G4PhysicsLogVector(); 254 255 std::ostringstream ost; 256 ost << FindDirectoryPath() << Z ; 257 std::ifstream filein(ost.str().c_str()); 258 if (!filein.is_open()) { 259 G4ExceptionDescription ed; 260 ed << "Data file <" << ost.str().c_str() 261 << "> is not opened!"; 262 G4Exception("G4NeutronElasticXS::Initialis 263 FatalException, ed, "Check G4P 264 return; 265 } 266 if(verboseLevel > 1) { 267 G4cout << "file " << ost.str() 268 << " is opened by G4NeutronElasticXS" << 269 } 270 271 // retrieve data from DB 272 if(!data[Z]->Retrieve(filein, true)) { 273 G4ExceptionDescription ed; 274 ed << "Data file <" << ost.str().c_str() 275 << "> is not retrieved!"; 276 G4Exception("G4NeutronElasticXS::Initialis 277 FatalException, ed, "Check G4PARTICLEXSDAT 278 return; 279 } 280 // smooth transition 281 G4double sig1 = (*(data[Z]))[data[Z]->GetVe 282 G4double ehigh = data[Z]->GetMaxEnergy(); 283 G4double sig2 = ggXsection->GetElasticEleme 284 ehigh, Z, aeff[ 285 coeff[Z] = (sig2 > 0.) ? sig1/sig2 : 1.0; 286 } 287