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