| 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 // G4ParticleHPThermalScatteringData
27 //
28 // 15-Nov-06 First implementation is done by T. Koi (SLAC/SCCS)
29 // 070625 implement clearCurrentXSData to fix memory leaking by T. Koi
30 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
31 // ---------------------------------------------------------------------
32
33 #include "G4ParticleHPThermalScatteringData.hh"
34
35 #include "G4ElementTable.hh"
36 #include "G4Neutron.hh"
37 #include "G4ParticleHPManager.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4Threading.hh"
40
41 #include <algorithm>
42 #include <list>
43
44 G4ParticleHPThermalScatteringData::G4ParticleHPThermalScatteringData()
45 : G4VCrossSectionDataSet("NeutronHPThermalScatteringData")
46 {
47 // Upper limit of neutron energy
48 emax = 4 * eV;
49 SetMinKinEnergy(0 * MeV);
50 SetMaxKinEnergy(emax);
51
52 ke_cache = 0.0;
53 xs_cache = 0.0;
54 element_cache = nullptr;
55 material_cache = nullptr;
56
57 indexOfThermalElement.clear();
58
59 names = new G4ParticleHPThermalScatteringNames();
60 }
61
62 G4ParticleHPThermalScatteringData::~G4ParticleHPThermalScatteringData()
63 {
64 clearCurrentXSData();
65
66 delete names;
67 }
68
69 G4bool G4ParticleHPThermalScatteringData::IsIsoApplicable(const G4DynamicParticle* dp, G4int /*Z*/,
70 G4int /*A*/, const G4Element* element,
71 const G4Material* material)
72 {
73 G4double eKin = dp->GetKineticEnergy();
74 if (eKin > 4.0 * eV // GetMaxKinEnergy()
75 || eKin < 0 // GetMinKinEnergy()
76 || dp->GetDefinition() != G4Neutron::Neutron())
77 return false;
78
79 if (dic.find(std::pair<const G4Material*, const G4Element*>((G4Material*)nullptr, element))
80 != dic.end()
81 || dic.find(std::pair<const G4Material*, const G4Element*>(material, element)) != dic.end())
82 return true;
83
84 return false;
85 }
86
87 G4double G4ParticleHPThermalScatteringData::GetIsoCrossSection(const G4DynamicParticle* dp,
88 G4int /*Z*/, G4int /*A*/,
89 const G4Isotope* /*iso*/,
90 const G4Element* element,
91 const G4Material* material)
92 {
93 ke_cache = dp->GetKineticEnergy();
94 element_cache = element;
95 material_cache = material;
96 G4double xs = GetCrossSection(dp, element, material);
97 xs_cache = xs;
98 return xs;
99 }
100
101 void G4ParticleHPThermalScatteringData::clearCurrentXSData()
102 {
103 if (coherent != nullptr) {
104 for (auto it = coherent->cbegin(); it != coherent->cend(); ++it) {
105 if (it->second != nullptr) {
106 for (auto itt = it->second->cbegin(); itt != it->second->cend(); ++itt) {
107 delete itt->second;
108 }
109 }
110 delete it->second;
111 }
112 coherent->clear();
113 }
114
115 if (incoherent != nullptr) {
116 for (auto it = incoherent->cbegin(); it != incoherent->cend(); ++it) {
117 if (it->second != nullptr) {
118 for (auto itt = it->second->cbegin(); itt != it->second->cend(); ++itt) {
119 delete itt->second;
120 }
121 }
122 delete it->second;
123 }
124 incoherent->clear();
125 }
126
127 if (inelastic != nullptr) {
128 for (auto it = inelastic->cbegin(); it != inelastic->cend(); ++it) {
129 if (it->second != nullptr) {
130 for (auto itt = it->second->cbegin(); itt != it->second->cend(); ++itt) {
131 delete itt->second;
132 }
133 }
134 delete it->second;
135 }
136 inelastic->clear();
137 }
138 }
139
140 G4bool G4ParticleHPThermalScatteringData::IsApplicable(const G4DynamicParticle* aP,
141 const G4Element* anEle)
142 {
143 G4bool result = false;
144
145 G4double eKin = aP->GetKineticEnergy();
146 // Check energy
147 if (eKin < emax) {
148 // Check Particle Species
149 if (aP->GetDefinition() == G4Neutron::Neutron()) {
150 // anEle is one of Thermal elements
151 auto ie = (G4int)anEle->GetIndex();
152 for (int it : indexOfThermalElement) {
153 if (ie == it) return true;
154 }
155 }
156 }
157
158 return result;
159 }
160
161 void G4ParticleHPThermalScatteringData::BuildPhysicsTable(const G4ParticleDefinition& aP)
162 {
163 if (&aP != G4Neutron::Neutron())
164 throw G4HadronicException(__FILE__, __LINE__,
165 "Attempt to use NeutronHP data for particles other than neutrons!!!");
166
167 // std::map < std::pair < G4Material* , const G4Element* > , G4int > dic;
168 //
169 dic.clear();
170 if (G4Threading::IsMasterThread()) clearCurrentXSData();
171
172 std::map<G4String, G4int> co_dic;
173
174 // Searching Nist Materials
175 static G4ThreadLocal G4MaterialTable* theMaterialTable = nullptr;
176 if (theMaterialTable == nullptr) theMaterialTable = G4Material::GetMaterialTable();
177 std::size_t numberOfMaterials = G4Material::GetNumberOfMaterials();
178 for (std::size_t i = 0; i < numberOfMaterials; ++i) {
179 G4Material* material = (*theMaterialTable)[i];
180 auto numberOfElements = (G4int)material->GetNumberOfElements();
181 for (G4int j = 0; j < numberOfElements; ++j) {
182 const G4Element* element = material->GetElement(j);
183 if (names->IsThisThermalElement(material->GetName(), element->GetName())) {
184 G4int ts_ID_of_this_geometry;
185 G4String ts_ndl_name = names->GetTS_NDL_Name(material->GetName(), element->GetName());
186 if (co_dic.find(ts_ndl_name) != co_dic.cend()) {
187 ts_ID_of_this_geometry = co_dic.find(ts_ndl_name)->second;
188 }
189 else {
190 ts_ID_of_this_geometry = (G4int)co_dic.size();
191 co_dic.insert(std::pair<G4String, G4int>(ts_ndl_name, ts_ID_of_this_geometry));
192 }
193
194 dic.insert(std::pair<std::pair<G4Material*, const G4Element*>, G4int>(
195 std::pair<G4Material*, const G4Element*>(material, element), ts_ID_of_this_geometry));
196 }
197 }
198 }
199
200 // Searching TS Elements
201 auto theElementTable = G4Element::GetElementTable();
202 std::size_t numberOfElements = G4Element::GetNumberOfElements();
203
204 for (std::size_t i = 0; i < numberOfElements; ++i) {
205 const G4Element* element = (*theElementTable)[i];
206 if (names->IsThisThermalElement(element->GetName())) {
207 if (names->IsThisThermalElement(element->GetName())) {
208 G4int ts_ID_of_this_geometry;
209 G4String ts_ndl_name = names->GetTS_NDL_Name(element->GetName());
210 if (co_dic.find(ts_ndl_name) != co_dic.cend()) {
211 ts_ID_of_this_geometry = co_dic.find(ts_ndl_name)->second;
212 }
213 else {
214 ts_ID_of_this_geometry = (G4int)co_dic.size();
215 co_dic.insert(std::pair<G4String, G4int>(ts_ndl_name, ts_ID_of_this_geometry));
216 }
217
218 dic.insert(std::pair<std::pair<const G4Material*, const G4Element*>, G4int>(
219 std::pair<const G4Material*, const G4Element*>((G4Material*)nullptr, element),
220 ts_ID_of_this_geometry));
221 }
222 }
223 }
224
225 G4cout << G4endl;
226 G4cout << "Neutron HP Thermal Scattering Data: Following material-element pairs and/or elements "
227 "are registered."
228 << G4endl;
229 for (const auto& it : dic) {
230 if (it.first.first != nullptr) {
231 G4cout << "Material " << it.first.first->GetName() << " - Element "
232 << it.first.second->GetName() << ", internal thermal scattering id " << it.second
233 << G4endl;
234 }
235 else {
236 G4cout << "Element " << it.first.second->GetName() << ", internal thermal scattering id "
237 << it.second << G4endl;
238 }
239 }
240 G4cout << G4endl;
241
242 G4ParticleHPManager* hpmanager = G4ParticleHPManager::GetInstance();
243
244 coherent = hpmanager->GetThermalScatteringCoherentCrossSections();
245 incoherent = hpmanager->GetThermalScatteringIncoherentCrossSections();
246 inelastic = hpmanager->GetThermalScatteringInelasticCrossSections();
247
248 if (G4Threading::IsMasterThread()) {
249 if (coherent == nullptr)
250 coherent = new std::map<G4int, std::map<G4double, G4ParticleHPVector*>*>;
251 if (incoherent == nullptr)
252 incoherent = new std::map<G4int, std::map<G4double, G4ParticleHPVector*>*>;
253 if (inelastic == nullptr)
254 inelastic = new std::map<G4int, std::map<G4double, G4ParticleHPVector*>*>;
255
256 // Read Cross Section Data files
257
258 G4String dirName;
259 if (G4FindDataDir("G4NEUTRONHPDATA") == nullptr)
260 throw G4HadronicException(
261 __FILE__, __LINE__,
262 "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
263 G4String baseName = G4FindDataDir("G4NEUTRONHPDATA");
264
265 dirName = baseName + "/ThermalScattering";
266
267 G4String ndl_filename;
268 G4String full_name;
269
270 for (const auto& it : co_dic) {
271 ndl_filename = it.first;
272 G4int ts_ID = it.second;
273
274 // Coherent
275 full_name = dirName + "/Coherent/CrossSection/" + ndl_filename;
276 auto coh_amapTemp_EnergyCross = readData(full_name);
277 coherent->insert(std::pair<G4int, std::map<G4double, G4ParticleHPVector*>*>(
278 ts_ID, coh_amapTemp_EnergyCross));
279
280 // Incoherent
281 full_name = dirName + "/Incoherent/CrossSection/" + ndl_filename;
282 auto incoh_amapTemp_EnergyCross = readData(full_name);
283 incoherent->insert(std::pair<G4int, std::map<G4double, G4ParticleHPVector*>*>(
284 ts_ID, incoh_amapTemp_EnergyCross));
285
286 // Inelastic
287 full_name = dirName + "/Inelastic/CrossSection/" + ndl_filename;
288 auto inela_amapTemp_EnergyCross = readData(full_name);
289 inelastic->insert(std::pair<G4int, std::map<G4double, G4ParticleHPVector*>*>(
290 ts_ID, inela_amapTemp_EnergyCross));
291 }
292 hpmanager->RegisterThermalScatteringCoherentCrossSections(coherent);
293 hpmanager->RegisterThermalScatteringIncoherentCrossSections(incoherent);
294 hpmanager->RegisterThermalScatteringInelasticCrossSections(inelastic);
295 }
296 }
297
298 std::map<G4double, G4ParticleHPVector*>*
299 G4ParticleHPThermalScatteringData::readData(const G4String& full_name)
300 {
301 auto aData = new std::map<G4double, G4ParticleHPVector*>;
302
303 std::istringstream theChannel;
304 G4ParticleHPManager::GetInstance()->GetDataStream(full_name, theChannel);
305
306 G4int dummy;
307 while (theChannel >> dummy) // MF // Loop checking, 11.05.2015, T. Koi
308 {
309 theChannel >> dummy; // MT
310 G4double temp;
311 theChannel >> temp;
312 auto anEnergyCross = new G4ParticleHPVector;
313 G4int nData;
314 theChannel >> nData;
315 anEnergyCross->Init(theChannel, nData, eV, barn);
316 aData->insert(std::pair<G4double, G4ParticleHPVector*>(temp, anEnergyCross));
317 }
318
319 return aData;
320 }
321
322 void G4ParticleHPThermalScatteringData::DumpPhysicsTable(const G4ParticleDefinition& aP)
323 {
324 if (&aP != G4Neutron::Neutron())
325 throw G4HadronicException(__FILE__, __LINE__,
326 "Attempt to use NeutronHP data for particles other than neutrons!!!");
327 }
328
329 G4double G4ParticleHPThermalScatteringData::GetCrossSection(const G4DynamicParticle* aP,
330 const G4Element* anE,
331 const G4Material* aM)
332 {
333 G4double result = 0;
334
335 G4int ts_id = getTS_ID(aM, anE);
336
337 if (ts_id == -1) return result;
338
339 G4double aT = aM->GetTemperature();
340
341 G4double Xcoh = GetX(aP, aT, coherent->find(ts_id)->second);
342 G4double Xincoh = GetX(aP, aT, incoherent->find(ts_id)->second);
343 G4double Xinela = GetX(aP, aT, inelastic->find(ts_id)->second);
344
345 result = Xcoh + Xincoh + Xinela;
346
347 return result;
348 }
349
350 G4double G4ParticleHPThermalScatteringData::GetInelasticCrossSection(const G4DynamicParticle* aP,
351 const G4Element* anE,
352 const G4Material* aM)
353 {
354 G4double result = 0;
355 G4int ts_id = getTS_ID(aM, anE);
356 G4double aT = aM->GetTemperature();
357 result = GetX(aP, aT, inelastic->find(ts_id)->second);
358 return result;
359 }
360
361 G4double G4ParticleHPThermalScatteringData::GetCoherentCrossSection(const G4DynamicParticle* aP,
362 const G4Element* anE,
363 const G4Material* aM)
364 {
365 G4double result = 0;
366 G4int ts_id = getTS_ID(aM, anE);
367 G4double aT = aM->GetTemperature();
368 result = GetX(aP, aT, coherent->find(ts_id)->second);
369 return result;
370 }
371
372 G4double G4ParticleHPThermalScatteringData::GetIncoherentCrossSection(const G4DynamicParticle* aP,
373 const G4Element* anE,
374 const G4Material* aM)
375 {
376 G4double result = 0;
377 G4int ts_id = getTS_ID(aM, anE);
378 G4double aT = aM->GetTemperature();
379 result = GetX(aP, aT, incoherent->find(ts_id)->second);
380 return result;
381 }
382
383 G4int G4ParticleHPThermalScatteringData::getTS_ID(const G4Material* material,
384 const G4Element* element)
385 {
386 G4int result = -1;
387 if (dic.find(std::pair<const G4Material*, const G4Element*>((G4Material*)nullptr, element))
388 != dic.end())
389 return dic.find(std::pair<const G4Material*, const G4Element*>((G4Material*)nullptr, element))
390 ->second;
391 if (dic.find(std::pair<const G4Material*, const G4Element*>(material, element)) != dic.end())
392 return dic.find(std::pair<const G4Material*, const G4Element*>(material, element))->second;
393 return result;
394 }
395
396 G4double G4ParticleHPThermalScatteringData::
397 GetX(const G4DynamicParticle* aP, G4double aT,
398 std::map<G4double, G4ParticleHPVector*>* amapTemp_EnergyCross)
399 {
400 G4double result = 0;
401 if (amapTemp_EnergyCross->empty()) return result;
402
403 G4double eKinetic = aP->GetKineticEnergy();
404
405 if (amapTemp_EnergyCross->size() == 1) {
406 if (std::fabs(aT - amapTemp_EnergyCross->cbegin()->first) / amapTemp_EnergyCross->begin()->first
407 > 0.1)
408 {
409 G4cout
410 << "G4ParticleHPThermalScatteringData:: The temperature of material (" << aT / kelvin
411 << "K) is different more than 10% from temperature of thermal scattering file expected ("
412 << amapTemp_EnergyCross->begin()->first << "K). Result may not be reliable." << G4endl;
413 }
414 result = amapTemp_EnergyCross->begin()->second->GetXsec(eKinetic);
415 return result;
416 }
417
418 auto it = amapTemp_EnergyCross->cbegin();
419 for (it = amapTemp_EnergyCross->cbegin(); it != amapTemp_EnergyCross->cend(); ++it) {
420 if (aT < it->first) break;
421 }
422 if (it == amapTemp_EnergyCross->cbegin()) {
423 ++it; // lower than the first
424 }
425 else if (it == amapTemp_EnergyCross->cend()) {
426 --it; // upper than the last
427 }
428
429 G4double TH = it->first;
430 G4double XH = it->second->GetXsec(eKinetic);
431
432 if (it != amapTemp_EnergyCross->cbegin()) --it;
433 G4double TL = it->first;
434 G4double XL = it->second->GetXsec(eKinetic);
435
436 if (TH == TL) throw G4HadronicException(__FILE__, __LINE__, "Thermal Scattering Data Error!");
437
438 G4double T = aT;
439 G4double X = (XH - XL) / (TH - TL) * (T - TL) + XL;
440 result = X;
441
442 return result;
443 }
444
445 void G4ParticleHPThermalScatteringData::AddUserThermalScatteringFile(const G4String& nameG4Element,
446 const G4String& filename)
447 {
448 names->AddThermalElement(nameG4Element, filename);
449 }
450
451 void G4ParticleHPThermalScatteringData::CrossSectionDescription(std::ostream& outFile) const
452 {
453 outFile << "High Precision cross data based on thermal scattering data in evaluated nuclear data "
454 "libraries for neutrons below 5eV on specific materials\n";
455 }
456