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