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: G4NeutronCaptureXS.cc,v 1.8 2011-01-09 02:37:48 dennis Exp $ >> 27 // GEANT4 tag $Name: not supported by cvs2svn $ >> 28 // 26 // ------------------------------------------- 29 // ------------------------------------------------------------------- 27 // 30 // 28 // GEANT4 Class file 31 // GEANT4 Class file 29 // 32 // 30 // 33 // 31 // File name: G4NeutronCaptureXS 34 // File name: G4NeutronCaptureXS 32 // 35 // 33 // Author Ivantchenko, Geant4, 3-Aug-09 36 // Author Ivantchenko, Geant4, 3-Aug-09 34 // 37 // 35 // Modifications: 38 // Modifications: 36 // 39 // 37 40 38 #include <fstream> << 41 #include "G4HadronicException.hh" 39 #include <sstream> << 40 #include <thread> << 41 << 42 #include "G4SystemOfUnits.hh" << 43 #include "G4NeutronCaptureXS.hh" 42 #include "G4NeutronCaptureXS.hh" 44 #include "G4Material.hh" << 45 #include "G4Element.hh" 43 #include "G4Element.hh" 46 #include "G4PhysicsLogVector.hh" << 47 #include "G4DynamicParticle.hh" << 48 #include "G4ElementTable.hh" 44 #include "G4ElementTable.hh" 49 #include "G4IsotopeList.hh" << 45 #include "G4PhysicsLogVector.hh" 50 #include "G4HadronicParameters.hh" << 46 #include "G4PhysicsVector.hh" 51 #include "Randomize.hh" << 52 #include "G4Log.hh" << 53 #include "G4AutoLock.hh" << 54 << 55 G4ElementData* G4NeutronCaptureXS::data = null << 56 G4String G4NeutronCaptureXS::gDataDirectory = << 57 47 58 static std::once_flag applyOnce; << 48 #include <iostream> >> 49 #include <fstream> >> 50 #include <sstream> 59 51 60 namespace << 52 using namespace std; 61 { << 62 G4Mutex neutronCaptureXSMutex = G4MUTEX_INIT << 63 const G4int MAXZCAPTURE = 92; << 64 } << 65 53 66 G4NeutronCaptureXS::G4NeutronCaptureXS() 54 G4NeutronCaptureXS::G4NeutronCaptureXS() 67 : G4VCrossSectionDataSet(Default_Name()), << 55 : G4VCrossSectionDataSet("G4NeutronCaptureXS"), 68 emax(20*CLHEP::MeV), elimit(1.0e-5*CLHEP::e << 56 emax(20*MeV),maxZ(92) 69 { 57 { 70 verboseLevel = 0; << 58 // verboseLevel = 0; 71 if (verboseLevel > 0) { << 59 if(verboseLevel > 0){ 72 G4cout << "G4NeutronCaptureXS::G4NeutronC 60 G4cout << "G4NeutronCaptureXS::G4NeutronCaptureXS: Initialise for Z < " 73 << MAXZCAPTURE << G4endl; << 61 << maxZ + 1 << G4endl; 74 } << 75 logElimit = G4Log(elimit); << 76 if (nullptr == data) { << 77 data = new G4ElementData(MAXZCAPTURE+1); << 78 data->SetName("nCapture"); << 79 FindDirectoryPath(); << 80 } 62 } >> 63 data.resize(maxZ+1, 0); >> 64 isInitialized = false; 81 } 65 } 82 66 83 void G4NeutronCaptureXS::CrossSectionDescripti << 67 G4NeutronCaptureXS::~G4NeutronCaptureXS() 84 { 68 { 85 outFile << "G4NeutronCaptureXS calculates th << 69 for(G4int i=0; i<=maxZ; ++i) { 86 << "on nuclei using data from the hi << 70 delete data[i]; 87 << "These data are simplified and sm << 71 } 88 << "in order to reduce CPU time. G4N << 89 << "above 20 MeV for all targets. Fo << 90 << "Uranium is used.\n"; << 91 } 72 } 92 << 73 93 G4bool 74 G4bool 94 G4NeutronCaptureXS::IsElementApplicable(const << 75 G4NeutronCaptureXS::IsApplicable(const G4DynamicParticle*, 95 G4int, const G4Material*) << 76 const G4Element*) 96 { 77 { 97 return true; 78 return true; 98 } 79 } 99 80 100 G4bool 81 G4bool 101 G4NeutronCaptureXS::IsIsoApplicable(const G4Dy 82 G4NeutronCaptureXS::IsIsoApplicable(const G4DynamicParticle*, 102 G4int, G4int, << 83 G4int /*ZZ*/, G4int /*AA*/) 103 const G4Element*, const G4Material << 104 { 84 { 105 return true; << 85 return false; 106 } 86 } 107 87 108 G4double 88 G4double 109 G4NeutronCaptureXS::GetElementCrossSection(con << 89 G4NeutronCaptureXS::GetCrossSection(const G4DynamicParticle* aParticle, 110 G4int Z, const G4Material*) << 90 const G4Element* elm, >> 91 G4double) 111 { 92 { 112 G4double xs = 0.0; 93 G4double xs = 0.0; 113 G4double ekin = aParticle->GetKineticEnergy( 94 G4double ekin = aParticle->GetKineticEnergy(); 114 if (ekin < emax) { << 95 if(ekin > emax) { return xs; } 115 xs = ElementCrossSection(ekin, aParticle-> << 96 const G4double elimit = 1.0e-10*eV; 116 } << 97 if(ekin < elimit) { ekin = elimit; } 117 return xs; << 98 118 } << 99 G4int Z = G4int(elm->GetZ()); >> 100 G4PhysicsVector* pv = data[Z]; >> 101 >> 102 // element was not initialised >> 103 if(!pv) { >> 104 Initialise(Z); >> 105 pv = data[Z]; >> 106 if(!pv) { return xs; } >> 107 } >> 108 >> 109 G4int n = pv->GetVectorLength() - 1; >> 110 G4double e1 = pv->Energy(0); >> 111 G4double e2 = pv->Energy(n); >> 112 if(ekin < e1) { xs = (*pv)[0]*std::sqrt(e1/ekin); } >> 113 else if(ekin <= e2) { xs = pv->Value(ekin); } 119 114 120 G4double << 115 if(verboseLevel > 0){ 121 G4NeutronCaptureXS::ComputeCrossSectionPerElem << 116 G4cout << "ekin= " << ekin << ", xs= " << xs << G4endl; 122 const G4ParticleDefi << 123 const G4Element* elm << 124 const G4Material*) << 125 { << 126 G4double xs = 0.0; << 127 if (ekin < emax) { << 128 xs = ElementCrossSection(ekin, loge, elm-> << 129 } 117 } 130 return xs; 118 return xs; 131 } 119 } 132 120 133 G4double << 134 G4NeutronCaptureXS::ElementCrossSection(G4doub << 135 { << 136 G4int Z = std::min(ZZ, MAXZCAPTURE); << 137 G4double ekin = eKin; << 138 G4double logEkin = logE; << 139 if (ekin < elimit) { << 140 ekin = elimit; << 141 logEkin = logElimit; << 142 } << 143 << 144 auto pv = GetPhysicsVector(Z); << 145 const G4double e0 = pv->Energy(0); << 146 G4double xs = (ekin >= e0) ? pv->LogVectorVa << 147 : (*pv)[0]*std::sqrt(e0/ekin); << 148 << 149 #ifdef G4VERBOSE << 150 if (verboseLevel > 1){ << 151 G4cout << "Ekin= " << ekin/CLHEP::MeV << 152 << " ElmXScap(b)= " << xs/CLHEP::b << 153 } << 154 #endif << 155 return xs; << 156 } << 157 << 158 G4double << 159 G4NeutronCaptureXS::ComputeIsoCrossSection(G4d << 160 const G4ParticleDefinition* << 161 G4int Z, G4int A, << 162 const G4Isotope*, const G4E << 163 const G4Material*) << 164 { << 165 return IsoCrossSection(ekin, loge, Z, A); << 166 } << 167 << 168 G4double << 169 G4NeutronCaptureXS::GetIsoCrossSection(const G << 170 G4int Z, G4int A, << 171 const G4Isotope*, const G4Eleme << 172 const G4Material*) << 173 { << 174 return IsoCrossSection(aParticle->GetKinetic << 175 aParticle->GetLogKine << 176 Z, A); << 177 } << 178 << 179 G4double G4NeutronCaptureXS::IsoCrossSection(G << 180 G << 181 { << 182 G4double xs = 0.0; << 183 if (eKin > emax) { return xs; } << 184 << 185 G4int Z = std::min(ZZ, MAXZCAPTURE); << 186 G4double ekin = eKin; << 187 G4double logEkin = logE; << 188 if (ekin < elimit) { << 189 ekin = elimit; << 190 logEkin = logElimit; << 191 } << 192 << 193 auto pv = GetPhysicsVector(Z); << 194 if (pv == nullptr) { return xs; } << 195 << 196 // use isotope x-section if possible << 197 if (data->GetNumberOfComponents(Z) > 0) { << 198 G4PhysicsVector* pviso = data->GetComponen << 199 if(pviso != nullptr) { << 200 const G4double e0 = pviso->Energy(0); << 201 xs = (ekin >= e0) ? pviso->LogVectorValu << 202 : (*pviso)[0]*std::sqrt(e0/ekin); << 203 #ifdef G4VERBOSE << 204 if(verboseLevel > 0) { << 205 G4cout << "G4NeutronCaptureXS::IsoXS: Ekin(M << 206 << " xs(b)= " << xs/barn << 207 << " Z= " << Z << " A= " << A << G4 << 208 } << 209 #endif << 210 return xs; << 211 } << 212 } << 213 // isotope data are not available or applica << 214 const G4double e0 = pv->Energy(0); << 215 xs = (ekin >= e0) ? pv->LogVectorValue(ekin, << 216 : (*pv)[0]*std::sqrt(e0/ekin); << 217 #ifdef G4VERBOSE << 218 if (verboseLevel > 0) { << 219 G4cout << "G4NeutronCaptureXS::IsoXS: Ekin << 220 << " xs(b)= " << xs/barn << 221 << " Z= " << Z << " A= " << A << " no i << 222 } << 223 #endif << 224 return xs; << 225 } << 226 << 227 const G4Isotope* << 228 G4NeutronCaptureXS::SelectIsotope(const G4Elem << 229 G4double kinEnergy, G4double logE) << 230 { << 231 G4int nIso = (G4int)anElement->GetNumberOfIs << 232 const G4Isotope* iso = anElement->GetIsotope << 233 << 234 //G4cout << "SelectIsotope NIso= " << nIso < << 235 if(1 == nIso) { return iso; } << 236 << 237 // more than 1 isotope << 238 G4int Z = anElement->GetZasInt(); << 239 if (nullptr == data->GetElementData(Z)) { In << 240 << 241 const G4double* abundVector = anElement->Get << 242 G4double q = G4UniformRand(); << 243 G4double sum = 0.0; << 244 << 245 // is there isotope wise cross section? << 246 G4int j; << 247 if (Z > MAXZCAPTURE || 0 == data->GetNumberO << 248 for (j = 0; j<nIso; ++j) { << 249 sum += abundVector[j]; << 250 if(q <= sum) { << 251 iso = anElement->GetIsotope(j); << 252 break; << 253 } << 254 } << 255 return iso; << 256 } << 257 G4int nn = (G4int)temp.size(); << 258 if (nn < nIso) { temp.resize(nIso, 0.); } << 259 << 260 for (j=0; j<nIso; ++j) { << 261 sum += abundVector[j]*IsoCrossSection(kinE << 262 anElement->GetIsotope(j)->GetN()); << 263 temp[j] = sum; << 264 } << 265 sum *= q; << 266 for (j = 0; j<nIso; ++j) { << 267 if (temp[j] >= sum) { << 268 iso = anElement->GetIsotope(j); << 269 break; << 270 } << 271 } << 272 return iso; << 273 } << 274 << 275 void 121 void 276 G4NeutronCaptureXS::BuildPhysicsTable(const G4 122 G4NeutronCaptureXS::BuildPhysicsTable(const G4ParticleDefinition& p) 277 { 123 { 278 if (verboseLevel > 0){ << 124 if(verboseLevel > 0){ 279 G4cout << "G4NeutronCaptureXS::BuildPhysic 125 G4cout << "G4NeutronCaptureXS::BuildPhysicsTable for " 280 << p.GetParticleName() << G4endl; 126 << p.GetParticleName() << G4endl; 281 } 127 } 282 if (p.GetParticleName() != "neutron") { << 128 if(isInitialized || p.GetParticleName() != "neutron") { return; } 283 G4ExceptionDescription ed; << 129 isInitialized = true; 284 ed << p.GetParticleName() << " is a wrong << 285 << " only neutron is allowed"; << 286 G4Exception("G4NeutronCaptureXS::BuildPhys << 287 FatalException, ed, ""); << 288 return; << 289 } << 290 << 291 // it is possible re-initialisation for the << 292 const G4ElementTable* table = G4Element::Get << 293 130 294 // initialise static tables only once << 131 // check environment variable 295 std::call_once(applyOnce, [this]() { isIniti << 132 // Build the complete string identifying the file with the data set 296 << 133 char* path = getenv("G4NEUTRONXSDATA"); 297 if (isInitializer) { << 134 if (!path){ 298 G4AutoLock l(&neutronCaptureXSMutex); << 135 G4cout << "G4NEUTRONXSDATA environment variable not set" << G4endl; 299 // Access to elements << 136 } 300 for ( auto const & elm : *table ) { << 137 301 G4int Z = std::max( 1, std::min( elm->Ge << 138 // Access to elements 302 if ( nullptr == data->GetElementData(Z) << 139 const G4ElementTable* theElmTable = G4Element::GetElementTable(); >> 140 size_t numOfElm = G4Element::GetNumberOfElements(); >> 141 if(numOfElm > 0) { >> 142 for(size_t i=0; i<numOfElm; ++i) { >> 143 G4int Z = G4int(((*theElmTable)[i])->GetZ()); >> 144 if(Z < 1) { Z = 1; } >> 145 else if(Z > maxZ) { Z = maxZ; } >> 146 //G4cout << "Z= " << Z << G4endl; >> 147 // Initialisation >> 148 if(!data[Z]) { Initialise(Z, path); } 303 } 149 } 304 l.unlock(); << 305 } << 306 << 307 // prepare isotope selection << 308 std::size_t nIso = temp.size(); << 309 for ( auto const & elm : *table ) { << 310 std::size_t n = elm->GetNumberOfIsotopes() << 311 if (n > nIso) { nIso = n; } << 312 } 150 } 313 temp.resize(nIso, 0.0); << 314 } 151 } 315 152 316 const G4String& G4NeutronCaptureXS::FindDirect << 153 void 317 { << 154 G4NeutronCaptureXS::DumpPhysicsTable(const G4ParticleDefinition&) 318 // build the complete string identifying the << 155 {} 319 if(gDataDirectory.empty()) { << 320 std::ostringstream ost; << 321 ost << G4HadronicParameters::Instance()->G << 322 gDataDirectory = ost.str(); << 323 } << 324 return gDataDirectory; << 325 } << 326 156 327 void G4NeutronCaptureXS::InitialiseOnFly(G4int << 157 void >> 158 G4NeutronCaptureXS::Initialise(G4int Z, const char* p) 328 { 159 { 329 G4AutoLock l(&neutronCaptureXSMutex); << 160 if(data[Z]) { return; } 330 Initialise(Z); << 161 const char* path = p; 331 l.unlock(); << 162 if(!p) { 332 } << 163 // check environment variable >> 164 // Build the complete string identifying the file with the data set >> 165 path = getenv("G4NEUTRONXSDATA"); >> 166 if (!path) { >> 167 if(verboseLevel > 1) { >> 168 G4cout << "G4NEUTRONXSDATA environment variable not set" << G4endl; >> 169 } >> 170 return; >> 171 } >> 172 } 333 173 334 void G4NeutronCaptureXS::Initialise(G4int Z) << 174 // upload data from file 335 { << 175 data[Z] = new G4PhysicsLogVector(); 336 if (nullptr != data->GetElementData(Z)) { re << 337 176 338 // upload element data << 339 std::ostringstream ost; 177 std::ostringstream ost; 340 ost << FindDirectoryPath() << Z ; << 178 ost << path << "/cap" << Z ; 341 G4PhysicsVector* v = RetrieveVector(ost, tru << 342 data->InitialiseForElement(Z, v); << 343 << 344 // upload isotope data << 345 G4bool noComp = true; << 346 if (amin[Z] < amax[Z]) { << 347 for(G4int A=amin[Z]; A<=amax[Z]; ++A) { << 348 std::ostringstream ost1; << 349 ost1 << gDataDirectory << Z << "_" << A; << 350 G4PhysicsVector* v1 = RetrieveVector(ost << 351 if (nullptr != v1) { << 352 if (noComp) { << 353 G4int nmax = amax[Z] - A + 1; << 354 data->InitialiseForComponent(Z, nmax); << 355 noComp = false; << 356 } << 357 data->AddComponent(Z, A, v1); << 358 } << 359 } << 360 } << 361 // no components case << 362 if (noComp) { data->InitialiseForComponent(Z << 363 } << 364 << 365 G4PhysicsVector* << 366 G4NeutronCaptureXS::RetrieveVector(std::ostrin << 367 { << 368 G4PhysicsLogVector* v = nullptr; << 369 std::ifstream filein(ost.str().c_str()); 179 std::ifstream filein(ost.str().c_str()); 370 if (!filein.is_open()) { << 180 if (!(filein)) { 371 if (warn) { << 181 G4cout << ost.str() << " is not opened by G4NeutronCaptureXS" << G4endl; 372 G4ExceptionDescription ed; << 182 throw G4HadronicException(__FILE__, __LINE__, 373 ed << "Data file <" << ost.str().c_str() << 183 "G4NeutronCaptureXS: no data sets registered"); 374 << "> is not opened!"; << 184 return; 375 G4Exception("G4NeutronCaptureXS::Retriev << 185 }else{ 376 FatalException, ed, "Check G4PARTICLEXSD << 186 if(verboseLevel > 1) { 377 } << 378 } else { << 379 if (verboseLevel > 1) { << 380 G4cout << "File " << ost.str() 187 G4cout << "File " << ost.str() 381 << " is opened by G4NeutronCaptureXS" < 188 << " is opened by G4NeutronCaptureXS" << G4endl; 382 } 189 } 383 // retrieve data from DB 190 // retrieve data from DB 384 v = new G4PhysicsLogVector(); << 191 data[Z]->Retrieve(filein, true); 385 if (!v->Retrieve(filein, true)) { << 192 } 386 G4ExceptionDescription ed; << 387 ed << "Data file <" << ost.str().c_str() << 388 << "> is not retrieved!"; << 389 G4Exception("G4NeutronCaptureXS::Retriev << 390 FatalException, ed, "Check G4PARTICLEXSD << 391 } << 392 } << 393 return v; << 394 } 193 } >> 194 395 195