| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // -------------------------------------------------------------------
27 //
28 // GEANT4 Class file
29 //
30 //
31 // File name: G4NeutronInelasticXS
32 //
33 // Author Ivantchenko, Geant4, 3-Aug-09
34 //
35
36 #include "G4NeutronInelasticXS.hh"
37 #include "G4Neutron.hh"
38 #include "G4DynamicParticle.hh"
39 #include "G4ElementTable.hh"
40 #include "G4Material.hh"
41 #include "G4Element.hh"
42 #include "G4PhysicsLogVector.hh"
43 #include "G4CrossSectionDataSetRegistry.hh"
44 #include "G4ComponentGGHadronNucleusXsc.hh"
45 #include "G4HadronicParameters.hh"
46 #include "Randomize.hh"
47 #include "G4Neutron.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "G4IsotopeList.hh"
50 #include "G4NuclearRadii.hh"
51 #include "G4AutoLock.hh"
52
53 #include <fstream>
54 #include <sstream>
55 #include <thread>
56
57 G4double G4NeutronInelasticXS::coeff[] = {1.0};
58 G4double G4NeutronInelasticXS::lowcoeff[] = {1.0};
59 G4ElementData* G4NeutronInelasticXS::data = nullptr;
60 G4String G4NeutronInelasticXS::gDataDirectory = "";
61
62 static std::once_flag applyOnce;
63
64 namespace
65 {
66 G4Mutex nInelasticXSMutex = G4MUTEX_INITIALIZER;
67 }
68
69 G4NeutronInelasticXS::G4NeutronInelasticXS()
70 : G4VCrossSectionDataSet(Default_Name()),
71 neutron(G4Neutron::Neutron()),
72 elimit(20*CLHEP::MeV),
73 lowElimit(1.0e-7*CLHEP::eV)
74 {
75 verboseLevel = 0;
76 if (verboseLevel > 0) {
77 G4cout << "G4NeutronInelasticXS::G4NeutronInelasticXS Initialise for Z < "
78 << MAXZINEL << G4endl;
79 }
80 loglowElimit = G4Log(lowElimit);
81 if (nullptr == data) {
82 data = new G4ElementData(MAXZINEL);
83 data->SetName("nInelastic");
84 FindDirectoryPath();
85 }
86 ggXsection =
87 G4CrossSectionDataSetRegistry::Instance()->GetComponentCrossSection("Glauber-Gribov");
88 if(ggXsection == nullptr)
89 ggXsection = new G4ComponentGGHadronNucleusXsc();
90
91 SetForAllAtomsAndEnergies(true);
92 }
93
94 void G4NeutronInelasticXS::CrossSectionDescription(std::ostream& outFile) const
95 {
96 outFile << "G4NeutronInelasticXS calculates the neutron inelastic scattering\n"
97 << "cross section on nuclei using data from the high precision\n"
98 << "neutron database. These data are simplified and smoothed over\n"
99 << "the resonance region in order to reduce CPU time.\n"
100 << "For high energy Glauber-Gribov cross section model is used\n";
101 }
102
103 G4bool
104 G4NeutronInelasticXS::IsElementApplicable(const G4DynamicParticle*,
105 G4int, const G4Material*)
106 {
107 return true;
108 }
109
110 G4bool
111 G4NeutronInelasticXS::IsIsoApplicable(const G4DynamicParticle*,
112 G4int, G4int,
113 const G4Element*, const G4Material*)
114 {
115 return true;
116 }
117
118 G4double
119 G4NeutronInelasticXS::GetElementCrossSection(const G4DynamicParticle* aParticle,
120 G4int Z, const G4Material*)
121 {
122 return ElementCrossSection(aParticle->GetKineticEnergy(),
123 aParticle->GetLogKineticEnergy(), Z);
124 }
125
126 G4double
127 G4NeutronInelasticXS::ComputeCrossSectionPerElement(G4double ekin, G4double loge,
128 const G4ParticleDefinition*,
129 const G4Element* elm,
130 const G4Material*)
131 {
132 return ElementCrossSection(ekin, loge, elm->GetZasInt());
133 }
134
135 G4double
136 G4NeutronInelasticXS::ElementCrossSection(G4double eKin, G4double logE, G4int ZZ)
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);
150
151 const G4double e0 = pv->Energy(0);
152 if (ekin <= e0) {
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 {
158 xs = coeff[Z]*ggXsection->GetInelasticElementCrossSection(neutron, ekin,
159 Z, aeff[Z]);
160 }
161
162 #ifdef G4VERBOSE
163 if(verboseLevel > 1) {
164 G4cout << "G4NeutronInelasticXS::ElementCrossSection Z= " << Z
165 << " Ekin(MeV)= " << ekin/CLHEP::MeV
166 << ", ElmXSinel(b)= " << xs/CLHEP::barn
167 << G4endl;
168 }
169 #endif
170 return xs;
171 }
172
173 G4double
174 G4NeutronInelasticXS::ComputeIsoCrossSection(G4double ekin, G4double loge,
175 const G4ParticleDefinition*,
176 G4int Z, G4int A,
177 const G4Isotope*, const G4Element*,
178 const G4Material*)
179 {
180 return IsoCrossSection(ekin, loge, Z, A);
181 }
182
183 G4double
184 G4NeutronInelasticXS::GetIsoCrossSection(const G4DynamicParticle* aParticle,
185 G4int Z, G4int A,
186 const G4Isotope*, const G4Element*,
187 const G4Material*)
188 {
189 return IsoCrossSection(aParticle->GetKineticEnergy(),
190 aParticle->GetLogKineticEnergy(), Z, A);
191 }
192
193 G4double
194 G4NeutronInelasticXS::IsoCrossSection(G4double eKin, G4double logE,
195 G4int ZZ, G4int A)
196 {
197 G4double xs;
198 G4int Z = std::min(ZZ, MAXZINEL-1);
199 G4double ekin = eKin;
200 G4double loge = logE;
201
202 /*
203 G4cout << "G4NeutronInelasticXS::IsoCrossSection Z= "
204 << Z << " A= " << A << G4endl;
205 G4cout << " Amin= " << amin[Z] << " Amax= " << amax[Z]
206 << " E(MeV)= " << ekin << " Ncomp="
207 << data->GetNumberOfComponents(Z) << G4endl;
208 */
209 GetPhysicsVector(Z);
210
211 // use isotope cross section if applicable
212 if (ekin <= elimit && data->GetNumberOfComponents(Z) > 0) {
213 auto pviso = data->GetComponentDataByID(Z, A);
214 if (nullptr != pviso) {
215 const G4double e0 = pviso->Energy(0);
216 if (ekin > e0) {
217 xs = pviso->LogVectorValue(ekin, loge);
218 } else {
219 xs = (*pviso)[0];
220 if (xs > 0.0) { xs *= std::sqrt(e0/ekin); }
221 }
222 #ifdef G4VERBOSE
223 if(verboseLevel > 1) {
224 G4cout << "G4NeutronInelasticXS::IsoXS: Ekin(MeV)= "
225 << ekin/CLHEP::MeV
226 << " xs(b)= " << xs/CLHEP::barn
227 << " Z= " << Z << " A= " << A << G4endl;
228 }
229 #endif
230 return xs;
231 }
232 }
233
234 // use element x-section
235 xs = ElementCrossSection(ekin, loge, Z)*A/aeff[Z];
236
237 #ifdef G4VERBOSE
238 if(verboseLevel > 1) {
239 G4cout << "G4NeutronInelasticXS::IsoXS: Z= " << Z << " A= " << A
240 << " Ekin(MeV)= " << ekin/CLHEP::MeV
241 << ", ElmXS(b)= " << xs/CLHEP::barn << G4endl;
242 }
243 #endif
244 return xs;
245 }
246
247 const G4Isotope* G4NeutronInelasticXS::SelectIsotope(
248 const G4Element* anElement, G4double kinEnergy, G4double logE)
249 {
250 std::size_t nIso = anElement->GetNumberOfIsotopes();
251 const G4Isotope* iso = anElement->GetIsotope(0);
252 if(1 == nIso) { return iso; }
253
254 // more than 1 isotope
255 G4int Z = anElement->GetZasInt();
256 if (nullptr == data->GetElementData(Z)) { InitialiseOnFly(Z); }
257
258 const G4double* abundVector = anElement->GetRelativeAbundanceVector();
259 G4double q = G4UniformRand();
260 G4double sum = 0.0;
261 std::size_t j;
262
263 // isotope wise cross section not available
264 if (Z >= MAXZINEL || 0 == data->GetNumberOfComponents(Z)) {
265 for (j=0; j<nIso; ++j) {
266 sum += abundVector[j];
267 if(q <= sum) {
268 iso = anElement->GetIsotope((G4int)j);
269 break;
270 }
271 }
272 return iso;
273 }
274
275 // use isotope cross sections
276 auto nn = temp.size();
277 if(nn < nIso) { temp.resize(nIso, 0.); }
278
279 for (j=0; j<nIso; ++j) {
280 // G4cout << j << "-th isotope " << anElement->GetIsotope(j)->GetN()
281 // << " abund= " << abundVector[j] << G4endl;
282 sum += abundVector[j]*IsoCrossSection(kinEnergy, logE, Z,
283 anElement->GetIsotope((G4int)j)->GetN());
284 temp[j] = sum;
285 }
286 sum *= q;
287 for (j = 0; j<nIso; ++j) {
288 if (temp[j] >= sum) {
289 iso = anElement->GetIsotope((G4int)j);
290 break;
291 }
292 }
293 return iso;
294 }
295
296 void
297 G4NeutronInelasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
298 {
299 if (verboseLevel > 0) {
300 G4cout << "G4NeutronInelasticXS::BuildPhysicsTable for "
301 << p.GetParticleName() << G4endl;
302 }
303 if (p.GetParticleName() != "neutron") {
304 G4ExceptionDescription ed;
305 ed << p.GetParticleName() << " is a wrong particle type -"
306 << " only neutron is allowed";
307 G4Exception("G4NeutronInelasticXS::BuildPhysicsTable(..)","had012",
308 FatalException, ed, "");
309 return;
310 }
311 // it is possible re-initialisation for the new run
312 const G4ElementTable* table = G4Element::GetElementTable();
313
314 // initialise static tables only once
315 std::call_once(applyOnce, [this]() { isInitializer = true; });
316
317 if (isInitializer) {
318 G4AutoLock l(&nInelasticXSMutex);
319
320 // Upload data for elements used in geometry
321 for ( auto const & elm : *table ) {
322 G4int Z = std::max( 1, std::min( elm->GetZasInt(), MAXZINEL-1) );
323 if ( nullptr == data->GetElementData(Z) ) { Initialise(Z); }
324 }
325 l.unlock();
326 }
327 // prepare isotope selection
328 std::size_t nIso = temp.size();
329 for ( auto const & elm : *table ) {
330 std::size_t n = elm->GetNumberOfIsotopes();
331 if (n > nIso) { nIso = n; }
332 }
333 temp.resize(nIso, 0.0);
334 }
335
336 const G4String& G4NeutronInelasticXS::FindDirectoryPath()
337 {
338 // build the complete string identifying the file with the data set
339 if (gDataDirectory.empty()) {
340 std::ostringstream ost;
341 ost << G4HadronicParameters::Instance()->GetDirPARTICLEXS() << "/neutron/inel";
342 gDataDirectory = ost.str();
343 }
344 return gDataDirectory;
345 }
346
347 void G4NeutronInelasticXS::InitialiseOnFly(G4int Z)
348 {
349 G4AutoLock l(&nInelasticXSMutex);
350 Initialise(Z);
351 l.unlock();
352 }
353
354 void G4NeutronInelasticXS::Initialise(G4int Z)
355 {
356 if (nullptr != data->GetElementData(Z)) { return; }
357
358 // upload element data
359 std::ostringstream ost;
360 ost << FindDirectoryPath() << Z;
361 G4PhysicsVector* v = RetrieveVector(ost, true);
362 data->InitialiseForElement(Z, v);
363 if (verboseLevel > 1) {
364 G4cout << "G4NeutronInelasticXS::Initialise for Z= " << Z
365 << " A= " << aeff[Z] << " Amin= " << amin[Z]
366 << " Amax= " << amax[Z] << G4endl;
367 }
368 // upload isotope data
369 G4bool noComp = true;
370 if (amin[Z] < amax[Z]) {
371
372 for (G4int A=amin[Z]; A<=amax[Z]; ++A) {
373 std::ostringstream ost1;
374 ost1 << gDataDirectory << Z << "_" << A;
375 G4PhysicsVector* v1 = RetrieveVector(ost1, false);
376 if (nullptr != 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 }
385 }
386 // no components case
387 if (noComp) { data->InitialiseForComponent(Z, 0); }
388
389 // smooth transition
390 G4double sig1 = (*v)[v->GetVectorLength()-1];
391 G4double ehigh= v->GetMaxEnergy();
392 G4double sig2 = ggXsection->GetInelasticElementCrossSection(neutron,
393 ehigh, Z, aeff[Z]);
394 coeff[Z] = (sig2 > 0.) ? sig1/sig2 : 1.0;
395 }
396
397 G4PhysicsVector*
398 G4NeutronInelasticXS::RetrieveVector(std::ostringstream& ost, G4bool warn)
399 {
400 G4PhysicsLogVector* v = nullptr;
401 std::ifstream filein(ost.str().c_str());
402 if (!filein.is_open()) {
403 if(warn) {
404 G4ExceptionDescription ed;
405 ed << "Data file <" << ost.str().c_str()
406 << "> is not opened!";
407 G4Exception("G4NeutronInelasticXS::RetrieveVector(..)","had014",
408 FatalException, ed, "Check G4PARTICLEXSDATA");
409 }
410 } else {
411 if(verboseLevel > 1) {
412 G4cout << "File " << ost.str()
413 << " is opened by G4NeutronInelasticXS" << G4endl;
414 }
415 // retrieve data from DB
416 v = new G4PhysicsLogVector();
417 if(!v->Retrieve(filein, true)) {
418 G4ExceptionDescription ed;
419 ed << "Data file <" << ost.str().c_str()
420 << "> is not retrieved!";
421 G4Exception("G4NeutronInelasticXS::RetrieveVector(..)","had015",
422 FatalException, ed, "Check G4PARTICLEXSDATA");
423 }
424 }
425 return v;
426 }
427