| 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 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron transport model.
29 //
30 // 070523 add neglecting doppler broadening on the fly. T. Koi
31 // 070613 fix memory leaking by T. Koi
32 // 071002 enable cross section dump by T. Koi
33 // 080428 change checking point of "neglecting doppler broadening" flag
34 // from GetCrossSection to BuildPhysicsTable by T. Koi
35 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
36 //
37 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
38 // V. Ivanchenko July-2023 converted back
39 //
40 #include "G4NeutronHPCaptureData.hh"
41
42 #include "G4ElementTable.hh"
43 #include "G4HadronicParameters.hh"
44 #include "G4Neutron.hh"
45 #include "G4NucleiProperties.hh"
46 #include "G4ParticleHPData.hh"
47 #include "G4ParticleHPManager.hh"
48 #include "G4Element.hh"
49 #include "G4Material.hh"
50 #include "G4ParticleDefinition.hh"
51 #include "G4PhysicsTable.hh"
52 #include "G4PhysicalConstants.hh"
53 #include "G4Pow.hh"
54 #include "G4SystemOfUnits.hh"
55 #include "G4AutoLock.hh"
56
57 G4bool G4NeutronHPCaptureData::fLock = true;
58 G4PhysicsTable* G4NeutronHPCaptureData::theCrossSections = nullptr;
59
60 namespace
61 {
62 G4Mutex theHPCaptureData = G4MUTEX_INITIALIZER;
63 }
64
65 G4NeutronHPCaptureData::G4NeutronHPCaptureData()
66 : G4VCrossSectionDataSet("NeutronHPCaptureXS")
67 {
68 emax = 20.*CLHEP::MeV;
69 fManager = G4ParticleHPManager::GetInstance();
70 }
71
72 G4NeutronHPCaptureData::~G4NeutronHPCaptureData()
73 {
74 if (isFirst) {
75 if (nullptr != theCrossSections)
76 theCrossSections->clearAndDestroy();
77 delete theCrossSections;
78 theCrossSections = nullptr;
79 }
80 }
81
82 G4bool G4NeutronHPCaptureData::IsIsoApplicable(const G4DynamicParticle*,
83 G4int, G4int,
84 const G4Element*,
85 const G4Material*)
86 {
87 return true;
88 }
89
90 G4double G4NeutronHPCaptureData::GetIsoCrossSection(const G4DynamicParticle* dp,
91 G4int /*Z*/, G4int /*A*/,
92 const G4Isotope* /*iso*/,
93 const G4Element* element,
94 const G4Material* material)
95 {
96 if (dp->GetKineticEnergy() == ke_cache && element == element_cache &&
97 material == material_cache)
98 return xs_cache;
99
100 ke_cache = dp->GetKineticEnergy();
101 element_cache = element;
102 material_cache = material;
103 G4double xs = GetCrossSection(dp, element, material->GetTemperature());
104 xs_cache = xs;
105 return xs;
106 }
107
108 void G4NeutronHPCaptureData::BuildPhysicsTable(const G4ParticleDefinition& p)
109 {
110 // the choice of the first instance
111 if (fLock) {
112 G4AutoLock l(&theHPCaptureData);
113 if (fLock) {
114 isFirst = true;
115 fLock = false;
116 }
117 l.unlock();
118 }
119 if (!isFirst) { return; }
120 if (p.GetParticleName() != "neutron") {
121 G4ExceptionDescription ed;
122 ed << p.GetParticleName() << " is a wrong particle type -"
123 << " only neutron is allowed";
124 G4Exception("G4NeutronHPCaptureData::BuildPhysicsTable(..)","had012",
125 FatalException, ed, "");
126 return;
127 }
128
129 // initialisation for the first instance, others are locked
130 G4AutoLock l(&theHPCaptureData);
131 if (theCrossSections != nullptr) {
132 theCrossSections->clearAndDestroy();
133 delete theCrossSections;
134 }
135 std::size_t numberOfElements = G4Element::GetNumberOfElements();
136 theCrossSections = new G4PhysicsTable(numberOfElements);
137
138 // make a PhysicsVector for each element
139 auto theElementTable = G4Element::GetElementTable();
140 for (std::size_t i = 0; i < numberOfElements; ++i) {
141 auto elm = (*theElementTable)[i];
142 #ifdef G4VERBOSE
143 if (fManager->GetDEBUG()) {
144 G4cout << "ElementIndex " << elm->GetIndex() << " "
145 << elm->GetName() << G4endl;
146 }
147 #endif
148 G4PhysicsVector* physVec =
149 G4ParticleHPData::Instance(G4Neutron::Neutron())->MakePhysicsVector(elm, this);
150 theCrossSections->push_back(physVec);
151 }
152 fManager->RegisterCaptureCrossSections(theCrossSections);
153 l.unlock();
154 }
155
156 void G4NeutronHPCaptureData::DumpPhysicsTable(const G4ParticleDefinition&)
157 {
158 #ifdef G4VERBOSE
159 if (fManager->GetVerboseLevel() == 0) return;
160
161 //
162 // Dump element based cross section
163 // range 10e-5 eV to 20 MeV
164 // 10 point per decade
165 // in barn
166 //
167
168 G4cout << G4endl;
169 G4cout << G4endl;
170 G4cout << "Capture Cross Section of Neutron HP" << G4endl;
171 G4cout << "(Pointwise cross-section at 0 Kelvin.)" << G4endl;
172 G4cout << G4endl;
173 G4cout << "Name of Element" << G4endl;
174 G4cout << "Energy[eV] XS[barn]" << G4endl;
175 G4cout << G4endl;
176
177 std::size_t numberOfElements = G4Element::GetNumberOfElements();
178 auto theElementTable = G4Element::GetElementTable();
179
180 for (std::size_t i = 0; i < numberOfElements; ++i) {
181 G4cout << (*theElementTable)[i]->GetName() << G4endl;
182 G4cout << *((*theCrossSections)(i)) << G4endl;
183 }
184 #endif
185 }
186
187 G4double G4NeutronHPCaptureData::GetCrossSection(const G4DynamicParticle* aP,
188 const G4Element* anE,
189 G4double aT)
190 {
191 G4double result = 0;
192 auto idx = (G4int)anE->GetIndex();
193
194 // prepare neutron
195 G4double eKinetic = aP->GetKineticEnergy();
196 if (eKinetic >= emax) { return 0.0; }
197
198 // NEGLECT_DOPPLER
199 if (fManager->GetNeglectDoppler()) {
200 return (*((*theCrossSections)(idx))).Value(eKinetic);
201 }
202
203 G4ReactionProduct theNeutron(aP->GetDefinition());
204 theNeutron.SetMomentum(aP->GetMomentum());
205 theNeutron.SetKineticEnergy(eKinetic);
206
207 // prepare thermal nucleus
208 G4Nucleus aNuc;
209 G4int theA = anE->GetN();
210 G4int theZ = anE->GetZasInt();
211 G4double eleMass = G4NucleiProperties::GetNuclearMass(theA, theZ)
212 / CLHEP::neutron_mass_c2;
213
214 G4ReactionProduct boosted;
215 G4double aXsection;
216
217 // MC integration loop
218 G4int counter = 0;
219 G4double buffer = 0;
220 G4int size = G4int(std::max(10., aT / 60 * kelvin));
221 G4ThreeVector neutronVelocity =
222 1. / G4Neutron::Neutron()->GetPDGMass() * theNeutron.GetMomentum();
223 G4double neutronVMag = neutronVelocity.mag();
224
225 while (counter == 0 ||
226 std::abs(buffer - result / std::max(1, counter)) > 0.03 * buffer)
227 // Loop checking, 11.05.2015, T. Koi
228 {
229 if (counter != 0) buffer = result / counter;
230 while (counter < size) // Loop checking, 11.05.2015, T. Koi
231 {
232 ++counter;
233 G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus(eleMass, aT);
234 boosted.Lorentz(theNeutron, aThermalNuc);
235 G4double theEkin = boosted.GetKineticEnergy();
236 aXsection = (*((*theCrossSections)(idx))).Value(theEkin);
237 // velocity correction, or luminosity factor...
238 G4ThreeVector targetVelocity = 1. / aThermalNuc.GetMass() * aThermalNuc.GetMomentum();
239 aXsection *= (targetVelocity - neutronVelocity).mag() / neutronVMag;
240 result += aXsection;
241 }
242 size += size;
243 }
244 result /= counter;
245 /*
246 // Checking impact of G4NEUTRONHP_NEGLECT_DOPPLER
247 G4cout << " result " << result << " "
248 << (*((*theCrossSections)(index))).Value(eKinetic) << " "
249 << (*((*theCrossSections)(index))).Value(eKinetic) /result << G4endl;
250 */
251 return result;
252 }
253
254 void G4NeutronHPCaptureData::CrossSectionDescription(std::ostream& outF) const
255 {
256 outF << "High Precision cross data based on Evaluated Nuclear Data Files"
257 << " (ENDF) for radiative capture reaction of neutrons below 20 MeV";
258 }
259