| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4ParticleHPThermalScatteringData
27 //
28 // 15-Nov-06 First implementation is done by T
29 // 070625 implement clearCurrentXSData to fix
30 // P. Arce, June-2014 Conversion neutron_hp to
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::G4ParticleH
45 : G4VCrossSectionDataSet("NeutronHPThermalSc
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 G4ParticleHPThermalScatteringNam
60 }
61
62 G4ParticleHPThermalScatteringData::~G4Particle
63 {
64 clearCurrentXSData();
65
66 delete names;
67 }
68
69 G4bool G4ParticleHPThermalScatteringData::IsIs
70
71
72 {
73 G4double eKin = dp->GetKineticEnergy();
74 if (eKin > 4.0 * eV // GetMaxKinEnergy()
75 || eKin < 0 // GetMinKinEnergy()
76 || dp->GetDefinition() != G4Neutron::Neu
77 return false;
78
79 if (dic.find(std::pair<const G4Material*, co
80 != dic.end()
81 || dic.find(std::pair<const G4Material*,
82 return true;
83
84 return false;
85 }
86
87 G4double G4ParticleHPThermalScatteringData::Ge
88
89
90
91
92 {
93 ke_cache = dp->GetKineticEnergy();
94 element_cache = element;
95 material_cache = material;
96 G4double xs = GetCrossSection(dp, element, m
97 xs_cache = xs;
98 return xs;
99 }
100
101 void G4ParticleHPThermalScatteringData::clearC
102 {
103 if (coherent != nullptr) {
104 for (auto it = coherent->cbegin(); it != c
105 if (it->second != nullptr) {
106 for (auto itt = it->second->cbegin();
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 !=
117 if (it->second != nullptr) {
118 for (auto itt = it->second->cbegin();
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 !=
129 if (it->second != nullptr) {
130 for (auto itt = it->second->cbegin();
131 delete itt->second;
132 }
133 }
134 delete it->second;
135 }
136 inelastic->clear();
137 }
138 }
139
140 G4bool G4ParticleHPThermalScatteringData::IsAp
141
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::Neut
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::BuildP
162 {
163 if (&aP != G4Neutron::Neutron())
164 throw G4HadronicException(__FILE__, __LINE
165 "Attempt to use
166
167 // std::map < std::pair < G4Material* , cons
168 //
169 dic.clear();
170 if (G4Threading::IsMasterThread()) clearCurr
171
172 std::map<G4String, G4int> co_dic;
173
174 // Searching Nist Materials
175 static G4ThreadLocal G4MaterialTable* theMat
176 if (theMaterialTable == nullptr) theMaterial
177 std::size_t numberOfMaterials = G4Material::
178 for (std::size_t i = 0; i < numberOfMaterial
179 G4Material* material = (*theMaterialTable)
180 auto numberOfElements = (G4int)material->G
181 for (G4int j = 0; j < numberOfElements; ++
182 const G4Element* element = material->Get
183 if (names->IsThisThermalElement(material
184 G4int ts_ID_of_this_geometry;
185 G4String ts_ndl_name = names->GetTS_ND
186 if (co_dic.find(ts_ndl_name) != co_dic
187 ts_ID_of_this_geometry = co_dic.find
188 }
189 else {
190 ts_ID_of_this_geometry = (G4int)co_d
191 co_dic.insert(std::pair<G4String, G4
192 }
193
194 dic.insert(std::pair<std::pair<G4Mater
195 std::pair<G4Material*, const G4Eleme
196 }
197 }
198 }
199
200 // Searching TS Elements
201 auto theElementTable = G4Element::GetElement
202 std::size_t numberOfElements = G4Element::Ge
203
204 for (std::size_t i = 0; i < numberOfElements
205 const G4Element* element = (*theElementTab
206 if (names->IsThisThermalElement(element->G
207 if (names->IsThisThermalElement(element-
208 G4int ts_ID_of_this_geometry;
209 G4String ts_ndl_name = names->GetTS_ND
210 if (co_dic.find(ts_ndl_name) != co_dic
211 ts_ID_of_this_geometry = co_dic.find
212 }
213 else {
214 ts_ID_of_this_geometry = (G4int)co_d
215 co_dic.insert(std::pair<G4String, G4
216 }
217
218 dic.insert(std::pair<std::pair<const G
219 std::pair<const G4Material*, const G
220 ts_ID_of_this_geometry));
221 }
222 }
223 }
224
225 G4cout << G4endl;
226 G4cout << "Neutron HP Thermal Scattering Dat
227 "are registered."
228 << G4endl;
229 for (const auto& it : dic) {
230 if (it.first.first != nullptr) {
231 G4cout << "Material " << it.first.first-
232 << it.first.second->GetName() <<
233 << G4endl;
234 }
235 else {
236 G4cout << "Element " << it.first.second-
237 << it.second << G4endl;
238 }
239 }
240 G4cout << G4endl;
241
242 G4ParticleHPManager* hpmanager = G4ParticleH
243
244 coherent = hpmanager->GetThermalScatteringCo
245 incoherent = hpmanager->GetThermalScattering
246 inelastic = hpmanager->GetThermalScatteringI
247
248 if (G4Threading::IsMasterThread()) {
249 if (coherent == nullptr)
250 coherent = new std::map<G4int, std::map<
251 if (incoherent == nullptr)
252 incoherent = new std::map<G4int, std::ma
253 if (inelastic == nullptr)
254 inelastic = new std::map<G4int, std::map
255
256 // Read Cross Section Data files
257
258 G4String dirName;
259 if (G4FindDataDir("G4NEUTRONHPDATA") == nu
260 throw G4HadronicException(
261 __FILE__, __LINE__,
262 "Please setenv G4NEUTRONHPDATA to poin
263 G4String baseName = G4FindDataDir("G4NEUTR
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/CrossSe
276 auto coh_amapTemp_EnergyCross = readData
277 coherent->insert(std::pair<G4int, std::m
278 ts_ID, coh_amapTemp_EnergyCross));
279
280 // Incoherent
281 full_name = dirName + "/Incoherent/Cross
282 auto incoh_amapTemp_EnergyCross = readDa
283 incoherent->insert(std::pair<G4int, std:
284 ts_ID, incoh_amapTemp_EnergyCross));
285
286 // Inelastic
287 full_name = dirName + "/Inelastic/CrossS
288 auto inela_amapTemp_EnergyCross = readDa
289 inelastic->insert(std::pair<G4int, std::
290 ts_ID, inela_amapTemp_EnergyCross));
291 }
292 hpmanager->RegisterThermalScatteringCohere
293 hpmanager->RegisterThermalScatteringIncohe
294 hpmanager->RegisterThermalScatteringInelas
295 }
296 }
297
298 std::map<G4double, G4ParticleHPVector*>*
299 G4ParticleHPThermalScatteringData::readData(co
300 {
301 auto aData = new std::map<G4double, G4Partic
302
303 std::istringstream theChannel;
304 G4ParticleHPManager::GetInstance()->GetDataS
305
306 G4int dummy;
307 while (theChannel >> dummy) // MF // Loop c
308 {
309 theChannel >> dummy; // MT
310 G4double temp;
311 theChannel >> temp;
312 auto anEnergyCross = new G4ParticleHPVecto
313 G4int nData;
314 theChannel >> nData;
315 anEnergyCross->Init(theChannel, nData, eV,
316 aData->insert(std::pair<G4double, G4Partic
317 }
318
319 return aData;
320 }
321
322 void G4ParticleHPThermalScatteringData::DumpPh
323 {
324 if (&aP != G4Neutron::Neutron())
325 throw G4HadronicException(__FILE__, __LINE
326 "Attempt to use
327 }
328
329 G4double G4ParticleHPThermalScatteringData::Ge
330
331
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(
342 G4double Xincoh = GetX(aP, aT, incoherent->f
343 G4double Xinela = GetX(aP, aT, inelastic->fi
344
345 result = Xcoh + Xincoh + Xinela;
346
347 return result;
348 }
349
350 G4double G4ParticleHPThermalScatteringData::Ge
351
352
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)
358 return result;
359 }
360
361 G4double G4ParticleHPThermalScatteringData::Ge
362
363
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)-
369 return result;
370 }
371
372 G4double G4ParticleHPThermalScatteringData::Ge
373
374
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
380 return result;
381 }
382
383 G4int G4ParticleHPThermalScatteringData::getTS
384
385 {
386 G4int result = -1;
387 if (dic.find(std::pair<const G4Material*, co
388 != dic.end())
389 return dic.find(std::pair<const G4Material
390 ->second;
391 if (dic.find(std::pair<const G4Material*, co
392 return dic.find(std::pair<const G4Material
393 return result;
394 }
395
396 G4double G4ParticleHPThermalScatteringData::
397 GetX(const G4DynamicParticle* aP, G4double aT,
398 std::map<G4double, G4ParticleHPVector*>*
399 {
400 G4double result = 0;
401 if (amapTemp_EnergyCross->empty()) return re
402
403 G4double eKinetic = aP->GetKineticEnergy();
404
405 if (amapTemp_EnergyCross->size() == 1) {
406 if (std::fabs(aT - amapTemp_EnergyCross->c
407 > 0.1)
408 {
409 G4cout
410 << "G4ParticleHPThermalScatteringData:
411 << "K) is different more than 10% from
412 << amapTemp_EnergyCross->begin()->firs
413 }
414 result = amapTemp_EnergyCross->begin()->se
415 return result;
416 }
417
418 auto it = amapTemp_EnergyCross->cbegin();
419 for (it = amapTemp_EnergyCross->cbegin(); 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()) --
433 G4double TL = it->first;
434 G4double XL = it->second->GetXsec(eKinetic);
435
436 if (TH == TL) throw G4HadronicException(__FI
437
438 G4double T = aT;
439 G4double X = (XH - XL) / (TH - TL) * (T - TL
440 result = X;
441
442 return result;
443 }
444
445 void G4ParticleHPThermalScatteringData::AddUse
446
447 {
448 names->AddThermalElement(nameG4Element, file
449 }
450
451 void G4ParticleHPThermalScatteringData::CrossS
452 {
453 outFile << "High Precision cross data based
454 "libraries for neutrons below 5eV
455 }
456