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