| 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: G4NeutronCaptureXS
32 //
33 // Author Ivantchenko, Geant4, 3-Aug-09
34 //
35 // Modifications:
36 //
37
38 #include <fstream>
39 #include <sstream>
40 #include <thread>
41
42 #include "G4SystemOfUnits.hh"
43 #include "G4NeutronCaptureXS.hh"
44 #include "G4Material.hh"
45 #include "G4Element.hh"
46 #include "G4PhysicsLogVector.hh"
47 #include "G4DynamicParticle.hh"
48 #include "G4ElementTable.hh"
49 #include "G4IsotopeList.hh"
50 #include "G4HadronicParameters.hh"
51 #include "Randomize.hh"
52 #include "G4Log.hh"
53 #include "G4AutoLock.hh"
54
55 G4ElementData* G4NeutronCaptureXS::data = nullptr;
56 G4String G4NeutronCaptureXS::gDataDirectory = "";
57
58 static std::once_flag applyOnce;
59
60 namespace
61 {
62 G4Mutex neutronCaptureXSMutex = G4MUTEX_INITIALIZER;
63 const G4int MAXZCAPTURE = 92;
64 }
65
66 G4NeutronCaptureXS::G4NeutronCaptureXS()
67 : G4VCrossSectionDataSet(Default_Name()),
68 emax(20*CLHEP::MeV), elimit(1.0e-5*CLHEP::eV)
69 {
70 verboseLevel = 0;
71 if (verboseLevel > 0) {
72 G4cout << "G4NeutronCaptureXS::G4NeutronCaptureXS: Initialise for Z < "
73 << MAXZCAPTURE << G4endl;
74 }
75 logElimit = G4Log(elimit);
76 if (nullptr == data) {
77 data = new G4ElementData(MAXZCAPTURE+1);
78 data->SetName("nCapture");
79 FindDirectoryPath();
80 }
81 }
82
83 void G4NeutronCaptureXS::CrossSectionDescription(std::ostream& outFile) const
84 {
85 outFile << "G4NeutronCaptureXS calculates the neutron capture cross sections\n"
86 << "on nuclei using data from the high precision neutron database.\n"
87 << "These data are simplified and smoothed over the resonance region\n"
88 << "in order to reduce CPU time. G4NeutronCaptureXS is set to zero\n"
89 << "above 20 MeV for all targets. For Z > 92 the cross section of\n"
90 << "Uranium is used.\n";
91 }
92
93 G4bool
94 G4NeutronCaptureXS::IsElementApplicable(const G4DynamicParticle*,
95 G4int, const G4Material*)
96 {
97 return true;
98 }
99
100 G4bool
101 G4NeutronCaptureXS::IsIsoApplicable(const G4DynamicParticle*,
102 G4int, G4int,
103 const G4Element*, const G4Material*)
104 {
105 return true;
106 }
107
108 G4double
109 G4NeutronCaptureXS::GetElementCrossSection(const G4DynamicParticle* aParticle,
110 G4int Z, const G4Material*)
111 {
112 G4double xs = 0.0;
113 G4double ekin = aParticle->GetKineticEnergy();
114 if (ekin < emax) {
115 xs = ElementCrossSection(ekin, aParticle->GetLogKineticEnergy(), Z);
116 }
117 return xs;
118 }
119
120 G4double
121 G4NeutronCaptureXS::ComputeCrossSectionPerElement(G4double ekin, G4double loge,
122 const G4ParticleDefinition*,
123 const G4Element* elm,
124 const G4Material*)
125 {
126 G4double xs = 0.0;
127 if (ekin < emax) {
128 xs = ElementCrossSection(ekin, loge, elm->GetZasInt());
129 }
130 return xs;
131 }
132
133 G4double
134 G4NeutronCaptureXS::ElementCrossSection(G4double eKin, G4double logE, G4int ZZ)
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->LogVectorValue(ekin, logEkin)
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::barn << G4endl;
153 }
154 #endif
155 return xs;
156 }
157
158 G4double
159 G4NeutronCaptureXS::ComputeIsoCrossSection(G4double ekin, G4double loge,
160 const G4ParticleDefinition*,
161 G4int Z, G4int A,
162 const G4Isotope*, const G4Element*,
163 const G4Material*)
164 {
165 return IsoCrossSection(ekin, loge, Z, A);
166 }
167
168 G4double
169 G4NeutronCaptureXS::GetIsoCrossSection(const G4DynamicParticle* aParticle,
170 G4int Z, G4int A,
171 const G4Isotope*, const G4Element*,
172 const G4Material*)
173 {
174 return IsoCrossSection(aParticle->GetKineticEnergy(),
175 aParticle->GetLogKineticEnergy(),
176 Z, A);
177 }
178
179 G4double G4NeutronCaptureXS::IsoCrossSection(G4double eKin, G4double logE,
180 G4int ZZ, G4int A)
181 {
182 G4double xs = 0.0;
183 if (eKin > emax) { return xs; }
184
185 G4int Z = std::min(ZZ, MAXZCAPTURE);
186 G4double ekin = eKin;
187 G4double logEkin = logE;
188 if (ekin < elimit) {
189 ekin = elimit;
190 logEkin = logElimit;
191 }
192
193 auto pv = GetPhysicsVector(Z);
194 if (pv == nullptr) { return xs; }
195
196 // use isotope x-section if possible
197 if (data->GetNumberOfComponents(Z) > 0) {
198 G4PhysicsVector* pviso = data->GetComponentDataByID(Z, A);
199 if(pviso != nullptr) {
200 const G4double e0 = pviso->Energy(0);
201 xs = (ekin >= e0) ? pviso->LogVectorValue(ekin, logEkin)
202 : (*pviso)[0]*std::sqrt(e0/ekin);
203 #ifdef G4VERBOSE
204 if(verboseLevel > 0) {
205 G4cout << "G4NeutronCaptureXS::IsoXS: Ekin(MeV)= " << ekin/MeV
206 << " xs(b)= " << xs/barn
207 << " Z= " << Z << " A= " << A << G4endl;
208 }
209 #endif
210 return xs;
211 }
212 }
213 // isotope data are not available or applicable
214 const G4double e0 = pv->Energy(0);
215 xs = (ekin >= e0) ? pv->LogVectorValue(ekin, logEkin)
216 : (*pv)[0]*std::sqrt(e0/ekin);
217 #ifdef G4VERBOSE
218 if (verboseLevel > 0) {
219 G4cout << "G4NeutronCaptureXS::IsoXS: Ekin(MeV)= " << ekin/MeV
220 << " xs(b)= " << xs/barn
221 << " Z= " << Z << " A= " << A << " no iso XS" << G4endl;
222 }
223 #endif
224 return xs;
225 }
226
227 const G4Isotope*
228 G4NeutronCaptureXS::SelectIsotope(const G4Element* anElement,
229 G4double kinEnergy, G4double logE)
230 {
231 G4int nIso = (G4int)anElement->GetNumberOfIsotopes();
232 const G4Isotope* iso = anElement->GetIsotope(0);
233
234 //G4cout << "SelectIsotope NIso= " << nIso << G4endl;
235 if(1 == nIso) { return iso; }
236
237 // more than 1 isotope
238 G4int Z = anElement->GetZasInt();
239 if (nullptr == data->GetElementData(Z)) { InitialiseOnFly(Z); }
240
241 const G4double* abundVector = anElement->GetRelativeAbundanceVector();
242 G4double q = G4UniformRand();
243 G4double sum = 0.0;
244
245 // is there isotope wise cross section?
246 G4int j;
247 if (Z > MAXZCAPTURE || 0 == data->GetNumberOfComponents(Z)) {
248 for (j = 0; j<nIso; ++j) {
249 sum += abundVector[j];
250 if(q <= sum) {
251 iso = anElement->GetIsotope(j);
252 break;
253 }
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(kinEnergy, logE, Z,
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 }
271 }
272 return iso;
273 }
274
275 void
276 G4NeutronCaptureXS::BuildPhysicsTable(const G4ParticleDefinition& p)
277 {
278 if (verboseLevel > 0){
279 G4cout << "G4NeutronCaptureXS::BuildPhysicsTable for "
280 << p.GetParticleName() << G4endl;
281 }
282 if (p.GetParticleName() != "neutron") {
283 G4ExceptionDescription ed;
284 ed << p.GetParticleName() << " is a wrong particle type -"
285 << " only neutron is allowed";
286 G4Exception("G4NeutronCaptureXS::BuildPhysicsTable(..)","had012",
287 FatalException, ed, "");
288 return;
289 }
290
291 // it is possible re-initialisation for the second run
292 const G4ElementTable* table = G4Element::GetElementTable();
293
294 // initialise static tables only once
295 std::call_once(applyOnce, [this]() { isInitializer = true; });
296
297 if (isInitializer) {
298 G4AutoLock l(&neutronCaptureXSMutex);
299 // Access to elements
300 for ( auto const & elm : *table ) {
301 G4int Z = std::max( 1, std::min( elm->GetZasInt(), MAXZCAPTURE) );
302 if ( nullptr == data->GetElementData(Z) ) { Initialise(Z); }
303 }
304 l.unlock();
305 }
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::FindDirectoryPath()
317 {
318 // build the complete string identifying the file with the data set
319 if(gDataDirectory.empty()) {
320 std::ostringstream ost;
321 ost << G4HadronicParameters::Instance()->GetDirPARTICLEXS() << "/neutron/cap";
322 gDataDirectory = ost.str();
323 }
324 return gDataDirectory;
325 }
326
327 void G4NeutronCaptureXS::InitialiseOnFly(G4int Z)
328 {
329 G4AutoLock l(&neutronCaptureXSMutex);
330 Initialise(Z);
331 l.unlock();
332 }
333
334 void G4NeutronCaptureXS::Initialise(G4int Z)
335 {
336 if (nullptr != data->GetElementData(Z)) { return; }
337
338 // upload element data
339 std::ostringstream ost;
340 ost << FindDirectoryPath() << Z ;
341 G4PhysicsVector* v = RetrieveVector(ost, true);
342 data->InitialiseForElement(Z, v);
343
344 // upload isotope data
345 G4bool noComp = true;
346 if (amin[Z] < amax[Z]) {
347 for(G4int A=amin[Z]; A<=amax[Z]; ++A) {
348 std::ostringstream ost1;
349 ost1 << gDataDirectory << Z << "_" << A;
350 G4PhysicsVector* v1 = RetrieveVector(ost1, false);
351 if (nullptr != v1) {
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 }
360 }
361 // no components case
362 if (noComp) { data->InitialiseForComponent(Z, 0); }
363 }
364
365 G4PhysicsVector*
366 G4NeutronCaptureXS::RetrieveVector(std::ostringstream& ost, G4bool warn)
367 {
368 G4PhysicsLogVector* v = nullptr;
369 std::ifstream filein(ost.str().c_str());
370 if (!filein.is_open()) {
371 if (warn) {
372 G4ExceptionDescription ed;
373 ed << "Data file <" << ost.str().c_str()
374 << "> is not opened!";
375 G4Exception("G4NeutronCaptureXS::RetrieveVector(..)","had014",
376 FatalException, ed, "Check G4PARTICLEXSDATA");
377 }
378 } else {
379 if (verboseLevel > 1) {
380 G4cout << "File " << ost.str()
381 << " is opened by G4NeutronCaptureXS" << G4endl;
382 }
383 // retrieve data from DB
384 v = new G4PhysicsLogVector();
385 if (!v->Retrieve(filein, true)) {
386 G4ExceptionDescription ed;
387 ed << "Data file <" << ost.str().c_str()
388 << "> is not retrieved!";
389 G4Exception("G4NeutronCaptureXS::RetrieveVector(..)","had015",
390 FatalException, ed, "Check G4PARTICLEXSDATA");
391 }
392 }
393 return v;
394 }
395