| 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 // particle_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 //
39 #include "G4ParticleHPInelasticData.hh"
40
41 #include "G4ElementTable.hh"
42 #include "G4HadronicParameters.hh"
43 #include "G4Neutron.hh"
44 #include "G4NucleiProperties.hh"
45 #include "G4ParticleHPData.hh"
46 #include "G4ParticleHPManager.hh"
47 #include "G4Pow.hh"
48
49 G4ParticleHPInelasticData::G4ParticleHPInelasticData(G4ParticleDefinition* projectile)
50 : G4VCrossSectionDataSet("")
51 {
52 const char* dataDirVariable;
53 G4String particleName;
54 if (projectile == G4Neutron::Neutron()) {
55 dataDirVariable = "G4NEUTRONHPDATA";
56 }
57 else if (projectile == G4Proton::Proton()) {
58 dataDirVariable = "G4PROTONHPDATA";
59 particleName = "Proton";
60 }
61 else if (projectile == G4Deuteron::Deuteron()) {
62 dataDirVariable = "G4DEUTERONHPDATA";
63 particleName = "Deuteron";
64 }
65 else if (projectile == G4Triton::Triton()) {
66 dataDirVariable = "G4TRITONHPDATA";
67 particleName = "Triton";
68 }
69 else if (projectile == G4He3::He3()) {
70 dataDirVariable = "G4HE3HPDATA";
71 particleName = "He3";
72 }
73 else if (projectile == G4Alpha::Alpha()) {
74 dataDirVariable = "G4ALPHAHPDATA";
75 particleName = "Alpha";
76 }
77 else {
78 G4String message(
79 "G4ParticleHPInelasticData may only be called for neutron, proton, deuteron, triton, He3 or "
80 "alpha, while it is called for "
81 + projectile->GetParticleName());
82 throw G4HadronicException(__FILE__, __LINE__, message.c_str());
83 }
84 G4String dataName = projectile->GetParticleName() + "HPInelasticXS";
85 dataName.at(0) = (char)std::toupper(dataName.at(0));
86 SetName(dataName);
87
88 if ((G4FindDataDir(dataDirVariable) == nullptr) && (G4FindDataDir("G4PARTICLEHPDATA") == nullptr))
89 {
90 G4String message("Please setenv G4PARTICLEHPDATA (recommended) or, at least setenv "
91 + G4String(dataDirVariable) + " to point to the "
92 + projectile->GetParticleName() + " cross-section files.");
93 throw G4HadronicException(__FILE__, __LINE__, message.c_str());
94 }
95
96 G4String dirName;
97 if (G4FindDataDir(dataDirVariable) != nullptr) {
98 dirName = G4FindDataDir(dataDirVariable);
99 }
100 else {
101 G4String baseName = G4FindDataDir("G4PARTICLEHPDATA");
102 dirName = baseName + "/" + particleName;
103 }
104 #ifdef G4VERBOSE
105 if (G4HadronicParameters::Instance()->GetVerboseLevel() > 0) {
106 G4cout << "@@@ G4ParticleHPInelasticData instantiated for particle "
107 << projectile->GetParticleName() << " data directory variable is " << dataDirVariable
108 << " pointing to " << dirName << G4endl;
109 }
110 #endif
111
112 SetMinKinEnergy(0 * CLHEP::MeV);
113 SetMaxKinEnergy(20 * CLHEP::MeV);
114
115 theCrossSections = nullptr;
116 theProjectile = projectile;
117
118 theHPData = nullptr;
119 instanceOfWorker = false;
120 if (G4Threading::IsMasterThread()) {
121 theHPData = new G4ParticleHPData(theProjectile);
122 }
123 else {
124 instanceOfWorker = true;
125 }
126 element_cache = nullptr;
127 material_cache = nullptr;
128 ke_cache = 0.0;
129 xs_cache = 0.0;
130 }
131
132 G4ParticleHPInelasticData::~G4ParticleHPInelasticData()
133 {
134 if (theCrossSections != nullptr && !instanceOfWorker) {
135 theCrossSections->clearAndDestroy();
136 delete theCrossSections;
137 theCrossSections = nullptr;
138 }
139 if (theHPData != nullptr && !instanceOfWorker) {
140 delete theHPData;
141 theHPData = nullptr;
142 }
143 }
144
145 G4bool G4ParticleHPInelasticData::IsIsoApplicable(const G4DynamicParticle* dp, G4int /*Z*/,
146 G4int /*A*/, const G4Element* /*elm*/,
147 const G4Material* /*mat*/)
148 {
149 G4double eKin = dp->GetKineticEnergy();
150 return eKin <= GetMaxKinEnergy() && eKin >= GetMinKinEnergy()
151 && dp->GetDefinition() == theProjectile;
152 }
153
154 G4double G4ParticleHPInelasticData::GetIsoCrossSection(const G4DynamicParticle* dp, G4int /*Z*/,
155 G4int /*A*/, const G4Isotope* /*iso*/,
156 const G4Element* element,
157 const G4Material* material)
158 {
159 if (dp->GetKineticEnergy() == ke_cache && element == element_cache && material == material_cache)
160 return xs_cache;
161
162 ke_cache = dp->GetKineticEnergy();
163 element_cache = element;
164 material_cache = material;
165 G4double xs = GetCrossSection(dp, element, material->GetTemperature());
166 xs_cache = xs;
167 return xs;
168 }
169
170 void G4ParticleHPInelasticData::BuildPhysicsTable(const G4ParticleDefinition& projectile)
171 {
172 if (G4Threading::IsWorkerThread()) {
173 theCrossSections = G4ParticleHPManager::GetInstance()->GetInelasticCrossSections(&projectile);
174 return;
175 }
176 if (theHPData == nullptr)
177 theHPData = G4ParticleHPData::Instance(const_cast<G4ParticleDefinition*>(&projectile));
178
179 std::size_t numberOfElements = G4Element::GetNumberOfElements();
180 if (theCrossSections == nullptr)
181 theCrossSections = new G4PhysicsTable(numberOfElements);
182 else
183 theCrossSections->clearAndDestroy();
184
185 // make a PhysicsVector for each element
186
187 auto theElementTable = G4Element::GetElementTable();
188 for (std::size_t i = 0; i < numberOfElements; ++i) {
189 G4PhysicsVector* physVec = theHPData->MakePhysicsVector((*theElementTable)[i], this);
190 theCrossSections->push_back(physVec);
191 }
192 G4ParticleHPManager::GetInstance()->RegisterInelasticCrossSections(&projectile, theCrossSections);
193 }
194
195 void G4ParticleHPInelasticData::DumpPhysicsTable(const G4ParticleDefinition&)
196 {
197 #ifdef G4VERBOSE
198 if (G4HadronicParameters::Instance()->GetVerboseLevel() == 0) return;
199
200 // Dump element based cross section
201 // range 10e-5 eV to 20 MeV
202 // 10 point per decade
203 // in barn
204
205 G4cout << G4endl;
206 G4cout << G4endl;
207 G4cout << "Inelastic Cross Section of Neutron HP" << G4endl;
208 G4cout << "(Pointwise cross-section at 0 Kelvin.)" << G4endl;
209 G4cout << G4endl;
210 G4cout << "Name of Element" << G4endl;
211 G4cout << "Energy[eV] XS[barn]" << G4endl;
212 G4cout << G4endl;
213
214 std::size_t numberOfElements = G4Element::GetNumberOfElements();
215 auto theElementTable = G4Element::GetElementTable();
216
217 for (std::size_t i = 0; i < numberOfElements; ++i) {
218 G4cout << (*theElementTable)[i]->GetName() << G4endl;
219
220 G4int ie = 0;
221
222 for (ie = 0; ie < 130; ie++) {
223 G4double eKinetic = 1.0e-5 * G4Pow::GetInstance()->powA(10.0, ie / 10.0) * CLHEP::eV;
224 G4bool outOfRange = false;
225
226 if (eKinetic < 20 * CLHEP::MeV) {
227 G4cout << eKinetic / CLHEP::eV << " "
228 << (*((*theCrossSections)(i))).GetValue(eKinetic, outOfRange) / CLHEP::barn
229 << G4endl;
230 }
231 }
232 G4cout << G4endl;
233 }
234 #endif
235 }
236
237 G4double G4ParticleHPInelasticData::GetCrossSection(const G4DynamicParticle* projectile,
238 const G4Element* anE, G4double aT)
239 {
240 G4double result = 0;
241 G4bool outOfRange;
242 std::size_t index = anE->GetIndex();
243
244 // prepare neutron
245 G4double eKinetic = projectile->GetKineticEnergy();
246
247 if (G4ParticleHPManager::GetInstance()->GetNeglectDoppler()) {
248 // NEGLECT_DOPPLER
249 G4double factor = 1.0;
250 if (eKinetic < aT * CLHEP::k_Boltzmann) {
251 // below 0.1 eV neutrons
252 // Have to do some, but now just igonre.
253 // Will take care after performance check.
254 // factor = factor * targetV;
255 }
256 return ((*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange)) * factor;
257 }
258
259 G4ReactionProduct theNeutron(projectile->GetDefinition());
260 theNeutron.SetMomentum(projectile->GetMomentum());
261 theNeutron.SetKineticEnergy(eKinetic);
262
263 // prepare thermal nucleus
264 G4Nucleus aNuc;
265 G4double eps = 0.0001;
266 G4double theA = anE->GetN();
267 G4double theZ = anE->GetZ();
268 G4double eleMass;
269 eleMass = G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theA + eps),
270 static_cast<G4int>(theZ + eps));
271
272 G4ReactionProduct boosted;
273 G4double aXsection;
274
275 // MC integration loop
276 G4int counter = 0;
277 G4int failCount = 0;
278 G4double buffer = 0;
279 G4int size = G4int(std::max(10., aT / 60 * CLHEP::kelvin));
280 G4ThreeVector neutronVelocity = 1. / theProjectile->GetPDGMass() * theNeutron.GetMomentum();
281 G4double neutronVMag = neutronVelocity.mag();
282
283 #ifndef G4PHP_DOPPLER_LOOP_ONCE
284 while (counter == 0
285 || std::abs(buffer - result / std::max(1, counter))
286 > 0.01 * buffer) // Loop checking, 11.05.2015, T. Koi
287 {
288 if (counter != 0) buffer = result / counter;
289 while (counter < size) // Loop checking, 11.05.2015, T. Koi
290 {
291 ++counter;
292 #endif
293 G4ReactionProduct aThermalNuc =
294 aNuc.GetThermalNucleus(eleMass / G4Neutron::Neutron()->GetPDGMass(), aT);
295 boosted.Lorentz(theNeutron, aThermalNuc);
296 G4double theEkin = boosted.GetKineticEnergy();
297 aXsection = (*((*theCrossSections)(index))).GetValue(theEkin, outOfRange);
298 if (aXsection < 0) {
299 if (failCount < 1000) {
300 ++failCount;
301 #ifndef G4PHP_DOPPLER_LOOP_ONCE
302 --counter;
303 continue;
304 #endif
305 }
306 else {
307 aXsection = 0;
308 }
309 }
310 // velocity correction.
311 G4ThreeVector targetVelocity = 1. / aThermalNuc.GetMass() * aThermalNuc.GetMomentum();
312 aXsection *= (targetVelocity - neutronVelocity).mag() / neutronVMag;
313 result += aXsection;
314 }
315 #ifndef G4PHP_DOPPLER_LOOP_ONCE
316 size += size;
317 }
318 result /= counter;
319 #endif
320
321 return result;
322 }
323
324 G4int G4ParticleHPInelasticData::GetVerboseLevel() const
325 {
326 return G4ParticleHPManager::GetInstance()->GetVerboseLevel();
327 }
328
329 void G4ParticleHPInelasticData::SetVerboseLevel(G4int newValue)
330 {
331 G4ParticleHPManager::GetInstance()->SetVerboseLevel(newValue);
332 }
333
334 void G4ParticleHPInelasticData::CrossSectionDescription(std::ostream& outFile) const
335 {
336 outFile << "Extension of High Precision cross section for inelastic reaction of proton, "
337 "deuteron, triton, He3 and alpha below 20MeV\n";
338 }
339