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