| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4ParticleHPThermalScatteringData << 26 // Thermal Neutron Scattering
>> 27 // Koi, Tatsumi (SCCS/SLAC)
27 // 28 //
>> 29 // Class Description
>> 30 // Cross Sections for a high precision (based on evaluated data
>> 31 // libraries) description of themal neutron scattering below 4 eV;
>> 32 // Based on Thermal neutron scattering files
>> 33 // from the evaluated nuclear data files ENDF/B-VI, Release2
>> 34 // To be used in your physics list in case you need this physics.
>> 35 // In this case you want to register an object of this class with
>> 36 // the corresponding process.
>> 37 // Class Description - End
>> 38
28 // 15-Nov-06 First implementation is done by T 39 // 15-Nov-06 First implementation is done by T. Koi (SLAC/SCCS)
29 // 070625 implement clearCurrentXSData to fix 40 // 070625 implement clearCurrentXSData to fix memory leaking by T. Koi
30 // P. Arce, June-2014 Conversion neutron_hp to 41 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
31 // ------------------------------------------- << 42 //
32 43
33 #include "G4ParticleHPThermalScatteringData.hh 44 #include "G4ParticleHPThermalScatteringData.hh"
34 45
35 #include "G4ElementTable.hh" 46 #include "G4ElementTable.hh"
36 #include "G4Neutron.hh" 47 #include "G4Neutron.hh"
37 #include "G4ParticleHPManager.hh" 48 #include "G4ParticleHPManager.hh"
38 #include "G4SystemOfUnits.hh" 49 #include "G4SystemOfUnits.hh"
39 #include "G4Threading.hh" 50 #include "G4Threading.hh"
40 51
41 #include <algorithm> 52 #include <algorithm>
42 #include <list> 53 #include <list>
43 54
44 G4ParticleHPThermalScatteringData::G4ParticleH 55 G4ParticleHPThermalScatteringData::G4ParticleHPThermalScatteringData()
45 : G4VCrossSectionDataSet("NeutronHPThermalSc 56 : G4VCrossSectionDataSet("NeutronHPThermalScatteringData")
46 { 57 {
47 // Upper limit of neutron energy 58 // Upper limit of neutron energy
48 emax = 4 * eV; 59 emax = 4 * eV;
49 SetMinKinEnergy(0 * MeV); 60 SetMinKinEnergy(0 * MeV);
50 SetMaxKinEnergy(emax); 61 SetMaxKinEnergy(emax);
51 62
52 ke_cache = 0.0; 63 ke_cache = 0.0;
53 xs_cache = 0.0; 64 xs_cache = 0.0;
54 element_cache = nullptr; 65 element_cache = nullptr;
55 material_cache = nullptr; 66 material_cache = nullptr;
56 67
57 indexOfThermalElement.clear(); 68 indexOfThermalElement.clear();
58 69
59 names = new G4ParticleHPThermalScatteringNam 70 names = new G4ParticleHPThermalScatteringNames();
60 } 71 }
61 72
62 G4ParticleHPThermalScatteringData::~G4Particle 73 G4ParticleHPThermalScatteringData::~G4ParticleHPThermalScatteringData()
63 { 74 {
64 clearCurrentXSData(); 75 clearCurrentXSData();
65 76
66 delete names; 77 delete names;
67 } 78 }
68 79
69 G4bool G4ParticleHPThermalScatteringData::IsIs 80 G4bool G4ParticleHPThermalScatteringData::IsIsoApplicable(const G4DynamicParticle* dp, G4int /*Z*/,
70 81 G4int /*A*/, const G4Element* element,
71 82 const G4Material* material)
72 { 83 {
73 G4double eKin = dp->GetKineticEnergy(); 84 G4double eKin = dp->GetKineticEnergy();
74 if (eKin > 4.0 * eV // GetMaxKinEnergy() 85 if (eKin > 4.0 * eV // GetMaxKinEnergy()
75 || eKin < 0 // GetMinKinEnergy() 86 || eKin < 0 // GetMinKinEnergy()
76 || dp->GetDefinition() != G4Neutron::Neu 87 || dp->GetDefinition() != G4Neutron::Neutron())
77 return false; 88 return false;
78 89
79 if (dic.find(std::pair<const G4Material*, co 90 if (dic.find(std::pair<const G4Material*, const G4Element*>((G4Material*)nullptr, element))
80 != dic.end() 91 != dic.end()
81 || dic.find(std::pair<const G4Material*, 92 || dic.find(std::pair<const G4Material*, const G4Element*>(material, element)) != dic.end())
82 return true; 93 return true;
83 94
84 return false; 95 return false;
85 } 96 }
86 97
87 G4double G4ParticleHPThermalScatteringData::Ge 98 G4double G4ParticleHPThermalScatteringData::GetIsoCrossSection(const G4DynamicParticle* dp,
88 99 G4int /*Z*/, G4int /*A*/,
89 100 const G4Isotope* /*iso*/,
90 101 const G4Element* element,
91 102 const G4Material* material)
92 { 103 {
93 ke_cache = dp->GetKineticEnergy(); 104 ke_cache = dp->GetKineticEnergy();
94 element_cache = element; 105 element_cache = element;
95 material_cache = material; 106 material_cache = material;
96 G4double xs = GetCrossSection(dp, element, m 107 G4double xs = GetCrossSection(dp, element, material);
97 xs_cache = xs; 108 xs_cache = xs;
98 return xs; 109 return xs;
99 } 110 }
100 111
101 void G4ParticleHPThermalScatteringData::clearC 112 void G4ParticleHPThermalScatteringData::clearCurrentXSData()
102 { 113 {
103 if (coherent != nullptr) { 114 if (coherent != nullptr) {
104 for (auto it = coherent->cbegin(); it != c 115 for (auto it = coherent->cbegin(); it != coherent->cend(); ++it) {
105 if (it->second != nullptr) { 116 if (it->second != nullptr) {
106 for (auto itt = it->second->cbegin(); 117 for (auto itt = it->second->cbegin(); itt != it->second->cend(); ++itt) {
107 delete itt->second; 118 delete itt->second;
108 } 119 }
109 } 120 }
110 delete it->second; 121 delete it->second;
111 } 122 }
112 coherent->clear(); 123 coherent->clear();
113 } 124 }
114 125
115 if (incoherent != nullptr) { 126 if (incoherent != nullptr) {
116 for (auto it = incoherent->cbegin(); it != 127 for (auto it = incoherent->cbegin(); it != incoherent->cend(); ++it) {
117 if (it->second != nullptr) { 128 if (it->second != nullptr) {
118 for (auto itt = it->second->cbegin(); 129 for (auto itt = it->second->cbegin(); itt != it->second->cend(); ++itt) {
119 delete itt->second; 130 delete itt->second;
120 } 131 }
121 } 132 }
122 delete it->second; 133 delete it->second;
123 } 134 }
124 incoherent->clear(); 135 incoherent->clear();
125 } 136 }
126 137
127 if (inelastic != nullptr) { 138 if (inelastic != nullptr) {
128 for (auto it = inelastic->cbegin(); it != 139 for (auto it = inelastic->cbegin(); it != inelastic->cend(); ++it) {
129 if (it->second != nullptr) { 140 if (it->second != nullptr) {
130 for (auto itt = it->second->cbegin(); 141 for (auto itt = it->second->cbegin(); itt != it->second->cend(); ++itt) {
131 delete itt->second; 142 delete itt->second;
132 } 143 }
133 } 144 }
134 delete it->second; 145 delete it->second;
135 } 146 }
136 inelastic->clear(); 147 inelastic->clear();
137 } 148 }
138 } 149 }
139 150
140 G4bool G4ParticleHPThermalScatteringData::IsAp 151 G4bool G4ParticleHPThermalScatteringData::IsApplicable(const G4DynamicParticle* aP,
141 152 const G4Element* anEle)
142 { 153 {
143 G4bool result = false; 154 G4bool result = false;
144 155
145 G4double eKin = aP->GetKineticEnergy(); 156 G4double eKin = aP->GetKineticEnergy();
146 // Check energy 157 // Check energy
147 if (eKin < emax) { 158 if (eKin < emax) {
148 // Check Particle Species 159 // Check Particle Species
149 if (aP->GetDefinition() == G4Neutron::Neut 160 if (aP->GetDefinition() == G4Neutron::Neutron()) {
150 // anEle is one of Thermal elements 161 // anEle is one of Thermal elements
151 auto ie = (G4int)anEle->GetIndex(); 162 auto ie = (G4int)anEle->GetIndex();
152 for (int it : indexOfThermalElement) { 163 for (int it : indexOfThermalElement) {
153 if (ie == it) return true; 164 if (ie == it) return true;
154 } 165 }
155 } 166 }
156 } 167 }
157 168
158 return result; 169 return result;
159 } 170 }
160 171
161 void G4ParticleHPThermalScatteringData::BuildP 172 void G4ParticleHPThermalScatteringData::BuildPhysicsTable(const G4ParticleDefinition& aP)
162 { 173 {
163 if (&aP != G4Neutron::Neutron()) 174 if (&aP != G4Neutron::Neutron())
164 throw G4HadronicException(__FILE__, __LINE 175 throw G4HadronicException(__FILE__, __LINE__,
165 "Attempt to use 176 "Attempt to use NeutronHP data for particles other than neutrons!!!");
166 177
167 // std::map < std::pair < G4Material* , cons 178 // std::map < std::pair < G4Material* , const G4Element* > , G4int > dic;
168 // 179 //
169 dic.clear(); 180 dic.clear();
170 if (G4Threading::IsMasterThread()) clearCurr 181 if (G4Threading::IsMasterThread()) clearCurrentXSData();
171 182
172 std::map<G4String, G4int> co_dic; 183 std::map<G4String, G4int> co_dic;
173 184
174 // Searching Nist Materials 185 // Searching Nist Materials
175 static G4ThreadLocal G4MaterialTable* theMat 186 static G4ThreadLocal G4MaterialTable* theMaterialTable = nullptr;
176 if (theMaterialTable == nullptr) theMaterial 187 if (theMaterialTable == nullptr) theMaterialTable = G4Material::GetMaterialTable();
177 std::size_t numberOfMaterials = G4Material:: 188 std::size_t numberOfMaterials = G4Material::GetNumberOfMaterials();
178 for (std::size_t i = 0; i < numberOfMaterial 189 for (std::size_t i = 0; i < numberOfMaterials; ++i) {
179 G4Material* material = (*theMaterialTable) 190 G4Material* material = (*theMaterialTable)[i];
180 auto numberOfElements = (G4int)material->G 191 auto numberOfElements = (G4int)material->GetNumberOfElements();
181 for (G4int j = 0; j < numberOfElements; ++ 192 for (G4int j = 0; j < numberOfElements; ++j) {
182 const G4Element* element = material->Get 193 const G4Element* element = material->GetElement(j);
183 if (names->IsThisThermalElement(material 194 if (names->IsThisThermalElement(material->GetName(), element->GetName())) {
184 G4int ts_ID_of_this_geometry; 195 G4int ts_ID_of_this_geometry;
185 G4String ts_ndl_name = names->GetTS_ND 196 G4String ts_ndl_name = names->GetTS_NDL_Name(material->GetName(), element->GetName());
186 if (co_dic.find(ts_ndl_name) != co_dic 197 if (co_dic.find(ts_ndl_name) != co_dic.cend()) {
187 ts_ID_of_this_geometry = co_dic.find 198 ts_ID_of_this_geometry = co_dic.find(ts_ndl_name)->second;
188 } 199 }
189 else { 200 else {
190 ts_ID_of_this_geometry = (G4int)co_d 201 ts_ID_of_this_geometry = (G4int)co_dic.size();
191 co_dic.insert(std::pair<G4String, G4 202 co_dic.insert(std::pair<G4String, G4int>(ts_ndl_name, ts_ID_of_this_geometry));
192 } 203 }
193 204
194 dic.insert(std::pair<std::pair<G4Mater 205 dic.insert(std::pair<std::pair<G4Material*, const G4Element*>, G4int>(
195 std::pair<G4Material*, const G4Eleme 206 std::pair<G4Material*, const G4Element*>(material, element), ts_ID_of_this_geometry));
196 } 207 }
197 } 208 }
198 } 209 }
199 210
200 // Searching TS Elements 211 // Searching TS Elements
201 auto theElementTable = G4Element::GetElement << 212 static G4ThreadLocal G4ElementTable* theElementTable = nullptr;
>> 213 if (theElementTable == nullptr) theElementTable = G4Element::GetElementTable();
202 std::size_t numberOfElements = G4Element::Ge 214 std::size_t numberOfElements = G4Element::GetNumberOfElements();
203 215
204 for (std::size_t i = 0; i < numberOfElements 216 for (std::size_t i = 0; i < numberOfElements; ++i) {
205 const G4Element* element = (*theElementTab 217 const G4Element* element = (*theElementTable)[i];
206 if (names->IsThisThermalElement(element->G 218 if (names->IsThisThermalElement(element->GetName())) {
207 if (names->IsThisThermalElement(element- 219 if (names->IsThisThermalElement(element->GetName())) {
208 G4int ts_ID_of_this_geometry; 220 G4int ts_ID_of_this_geometry;
209 G4String ts_ndl_name = names->GetTS_ND 221 G4String ts_ndl_name = names->GetTS_NDL_Name(element->GetName());
210 if (co_dic.find(ts_ndl_name) != co_dic 222 if (co_dic.find(ts_ndl_name) != co_dic.cend()) {
211 ts_ID_of_this_geometry = co_dic.find 223 ts_ID_of_this_geometry = co_dic.find(ts_ndl_name)->second;
212 } 224 }
213 else { 225 else {
214 ts_ID_of_this_geometry = (G4int)co_d 226 ts_ID_of_this_geometry = (G4int)co_dic.size();
215 co_dic.insert(std::pair<G4String, G4 227 co_dic.insert(std::pair<G4String, G4int>(ts_ndl_name, ts_ID_of_this_geometry));
216 } 228 }
217 229
218 dic.insert(std::pair<std::pair<const G 230 dic.insert(std::pair<std::pair<const G4Material*, const G4Element*>, G4int>(
219 std::pair<const G4Material*, const G 231 std::pair<const G4Material*, const G4Element*>((G4Material*)nullptr, element),
220 ts_ID_of_this_geometry)); 232 ts_ID_of_this_geometry));
221 } 233 }
222 } 234 }
223 } 235 }
224 236
225 G4cout << G4endl; 237 G4cout << G4endl;
226 G4cout << "Neutron HP Thermal Scattering Dat 238 G4cout << "Neutron HP Thermal Scattering Data: Following material-element pairs and/or elements "
227 "are registered." 239 "are registered."
228 << G4endl; 240 << G4endl;
229 for (const auto& it : dic) { 241 for (const auto& it : dic) {
230 if (it.first.first != nullptr) { 242 if (it.first.first != nullptr) {
231 G4cout << "Material " << it.first.first- 243 G4cout << "Material " << it.first.first->GetName() << " - Element "
232 << it.first.second->GetName() << 244 << it.first.second->GetName() << ", internal thermal scattering id " << it.second
233 << G4endl; 245 << G4endl;
234 } 246 }
235 else { 247 else {
236 G4cout << "Element " << it.first.second- 248 G4cout << "Element " << it.first.second->GetName() << ", internal thermal scattering id "
237 << it.second << G4endl; 249 << it.second << G4endl;
238 } 250 }
239 } 251 }
240 G4cout << G4endl; 252 G4cout << G4endl;
241 253
242 G4ParticleHPManager* hpmanager = G4ParticleH 254 G4ParticleHPManager* hpmanager = G4ParticleHPManager::GetInstance();
243 255
244 coherent = hpmanager->GetThermalScatteringCo 256 coherent = hpmanager->GetThermalScatteringCoherentCrossSections();
245 incoherent = hpmanager->GetThermalScattering 257 incoherent = hpmanager->GetThermalScatteringIncoherentCrossSections();
246 inelastic = hpmanager->GetThermalScatteringI 258 inelastic = hpmanager->GetThermalScatteringInelasticCrossSections();
247 259
248 if (G4Threading::IsMasterThread()) { 260 if (G4Threading::IsMasterThread()) {
249 if (coherent == nullptr) 261 if (coherent == nullptr)
250 coherent = new std::map<G4int, std::map< 262 coherent = new std::map<G4int, std::map<G4double, G4ParticleHPVector*>*>;
251 if (incoherent == nullptr) 263 if (incoherent == nullptr)
252 incoherent = new std::map<G4int, std::ma 264 incoherent = new std::map<G4int, std::map<G4double, G4ParticleHPVector*>*>;
253 if (inelastic == nullptr) 265 if (inelastic == nullptr)
254 inelastic = new std::map<G4int, std::map 266 inelastic = new std::map<G4int, std::map<G4double, G4ParticleHPVector*>*>;
255 267
256 // Read Cross Section Data files 268 // Read Cross Section Data files
257 269
258 G4String dirName; 270 G4String dirName;
259 if (G4FindDataDir("G4NEUTRONHPDATA") == nu 271 if (G4FindDataDir("G4NEUTRONHPDATA") == nullptr)
260 throw G4HadronicException( 272 throw G4HadronicException(
261 __FILE__, __LINE__, 273 __FILE__, __LINE__,
262 "Please setenv G4NEUTRONHPDATA to poin 274 "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
263 G4String baseName = G4FindDataDir("G4NEUTR 275 G4String baseName = G4FindDataDir("G4NEUTRONHPDATA");
264 276
265 dirName = baseName + "/ThermalScattering"; 277 dirName = baseName + "/ThermalScattering";
266 278
267 G4String ndl_filename; 279 G4String ndl_filename;
268 G4String full_name; 280 G4String full_name;
269 281
270 for (const auto& it : co_dic) { 282 for (const auto& it : co_dic) {
271 ndl_filename = it.first; 283 ndl_filename = it.first;
272 G4int ts_ID = it.second; 284 G4int ts_ID = it.second;
273 285
274 // Coherent 286 // Coherent
275 full_name = dirName + "/Coherent/CrossSe 287 full_name = dirName + "/Coherent/CrossSection/" + ndl_filename;
276 auto coh_amapTemp_EnergyCross = readData 288 auto coh_amapTemp_EnergyCross = readData(full_name);
277 coherent->insert(std::pair<G4int, std::m 289 coherent->insert(std::pair<G4int, std::map<G4double, G4ParticleHPVector*>*>(
278 ts_ID, coh_amapTemp_EnergyCross)); 290 ts_ID, coh_amapTemp_EnergyCross));
279 291
280 // Incoherent 292 // Incoherent
281 full_name = dirName + "/Incoherent/Cross 293 full_name = dirName + "/Incoherent/CrossSection/" + ndl_filename;
282 auto incoh_amapTemp_EnergyCross = readDa 294 auto incoh_amapTemp_EnergyCross = readData(full_name);
283 incoherent->insert(std::pair<G4int, std: 295 incoherent->insert(std::pair<G4int, std::map<G4double, G4ParticleHPVector*>*>(
284 ts_ID, incoh_amapTemp_EnergyCross)); 296 ts_ID, incoh_amapTemp_EnergyCross));
285 297
286 // Inelastic 298 // Inelastic
287 full_name = dirName + "/Inelastic/CrossS 299 full_name = dirName + "/Inelastic/CrossSection/" + ndl_filename;
288 auto inela_amapTemp_EnergyCross = readDa 300 auto inela_amapTemp_EnergyCross = readData(full_name);
289 inelastic->insert(std::pair<G4int, std:: 301 inelastic->insert(std::pair<G4int, std::map<G4double, G4ParticleHPVector*>*>(
290 ts_ID, inela_amapTemp_EnergyCross)); 302 ts_ID, inela_amapTemp_EnergyCross));
291 } 303 }
292 hpmanager->RegisterThermalScatteringCohere 304 hpmanager->RegisterThermalScatteringCoherentCrossSections(coherent);
293 hpmanager->RegisterThermalScatteringIncohe 305 hpmanager->RegisterThermalScatteringIncoherentCrossSections(incoherent);
294 hpmanager->RegisterThermalScatteringInelas 306 hpmanager->RegisterThermalScatteringInelasticCrossSections(inelastic);
295 } 307 }
296 } 308 }
297 309
298 std::map<G4double, G4ParticleHPVector*>* 310 std::map<G4double, G4ParticleHPVector*>*
299 G4ParticleHPThermalScatteringData::readData(co << 311 G4ParticleHPThermalScatteringData::readData(G4String full_name)
300 { 312 {
301 auto aData = new std::map<G4double, G4Partic 313 auto aData = new std::map<G4double, G4ParticleHPVector*>;
302 314
303 std::istringstream theChannel; 315 std::istringstream theChannel;
304 G4ParticleHPManager::GetInstance()->GetDataS 316 G4ParticleHPManager::GetInstance()->GetDataStream(full_name, theChannel);
305 317
306 G4int dummy; 318 G4int dummy;
307 while (theChannel >> dummy) // MF // Loop c 319 while (theChannel >> dummy) // MF // Loop checking, 11.05.2015, T. Koi
308 { 320 {
309 theChannel >> dummy; // MT 321 theChannel >> dummy; // MT
310 G4double temp; 322 G4double temp;
311 theChannel >> temp; 323 theChannel >> temp;
312 auto anEnergyCross = new G4ParticleHPVecto 324 auto anEnergyCross = new G4ParticleHPVector;
313 G4int nData; 325 G4int nData;
314 theChannel >> nData; 326 theChannel >> nData;
315 anEnergyCross->Init(theChannel, nData, eV, 327 anEnergyCross->Init(theChannel, nData, eV, barn);
316 aData->insert(std::pair<G4double, G4Partic 328 aData->insert(std::pair<G4double, G4ParticleHPVector*>(temp, anEnergyCross));
317 } 329 }
318 330
319 return aData; 331 return aData;
320 } 332 }
321 333
322 void G4ParticleHPThermalScatteringData::DumpPh 334 void G4ParticleHPThermalScatteringData::DumpPhysicsTable(const G4ParticleDefinition& aP)
323 { 335 {
324 if (&aP != G4Neutron::Neutron()) 336 if (&aP != G4Neutron::Neutron())
325 throw G4HadronicException(__FILE__, __LINE 337 throw G4HadronicException(__FILE__, __LINE__,
326 "Attempt to use 338 "Attempt to use NeutronHP data for particles other than neutrons!!!");
327 } 339 }
328 340
329 G4double G4ParticleHPThermalScatteringData::Ge 341 G4double G4ParticleHPThermalScatteringData::GetCrossSection(const G4DynamicParticle* aP,
330 342 const G4Element* anE,
331 343 const G4Material* aM)
332 { 344 {
333 G4double result = 0; 345 G4double result = 0;
334 346
335 G4int ts_id = getTS_ID(aM, anE); 347 G4int ts_id = getTS_ID(aM, anE);
336 348
337 if (ts_id == -1) return result; 349 if (ts_id == -1) return result;
338 350
339 G4double aT = aM->GetTemperature(); 351 G4double aT = aM->GetTemperature();
340 352
341 G4double Xcoh = GetX(aP, aT, coherent->find( 353 G4double Xcoh = GetX(aP, aT, coherent->find(ts_id)->second);
342 G4double Xincoh = GetX(aP, aT, incoherent->f 354 G4double Xincoh = GetX(aP, aT, incoherent->find(ts_id)->second);
343 G4double Xinela = GetX(aP, aT, inelastic->fi 355 G4double Xinela = GetX(aP, aT, inelastic->find(ts_id)->second);
344 356
345 result = Xcoh + Xincoh + Xinela; 357 result = Xcoh + Xincoh + Xinela;
346 358
347 return result; 359 return result;
348 } 360 }
349 361
350 G4double G4ParticleHPThermalScatteringData::Ge 362 G4double G4ParticleHPThermalScatteringData::GetInelasticCrossSection(const G4DynamicParticle* aP,
351 363 const G4Element* anE,
352 364 const G4Material* aM)
353 { 365 {
354 G4double result = 0; 366 G4double result = 0;
355 G4int ts_id = getTS_ID(aM, anE); 367 G4int ts_id = getTS_ID(aM, anE);
356 G4double aT = aM->GetTemperature(); 368 G4double aT = aM->GetTemperature();
357 result = GetX(aP, aT, inelastic->find(ts_id) 369 result = GetX(aP, aT, inelastic->find(ts_id)->second);
358 return result; 370 return result;
359 } 371 }
360 372
361 G4double G4ParticleHPThermalScatteringData::Ge 373 G4double G4ParticleHPThermalScatteringData::GetCoherentCrossSection(const G4DynamicParticle* aP,
362 374 const G4Element* anE,
363 375 const G4Material* aM)
364 { 376 {
365 G4double result = 0; 377 G4double result = 0;
366 G4int ts_id = getTS_ID(aM, anE); 378 G4int ts_id = getTS_ID(aM, anE);
367 G4double aT = aM->GetTemperature(); 379 G4double aT = aM->GetTemperature();
368 result = GetX(aP, aT, coherent->find(ts_id)- 380 result = GetX(aP, aT, coherent->find(ts_id)->second);
369 return result; 381 return result;
370 } 382 }
371 383
372 G4double G4ParticleHPThermalScatteringData::Ge 384 G4double G4ParticleHPThermalScatteringData::GetIncoherentCrossSection(const G4DynamicParticle* aP,
373 385 const G4Element* anE,
374 386 const G4Material* aM)
375 { 387 {
376 G4double result = 0; 388 G4double result = 0;
377 G4int ts_id = getTS_ID(aM, anE); 389 G4int ts_id = getTS_ID(aM, anE);
378 G4double aT = aM->GetTemperature(); 390 G4double aT = aM->GetTemperature();
379 result = GetX(aP, aT, incoherent->find(ts_id 391 result = GetX(aP, aT, incoherent->find(ts_id)->second);
380 return result; 392 return result;
381 } 393 }
382 394
383 G4int G4ParticleHPThermalScatteringData::getTS 395 G4int G4ParticleHPThermalScatteringData::getTS_ID(const G4Material* material,
384 396 const G4Element* element)
385 { 397 {
386 G4int result = -1; 398 G4int result = -1;
387 if (dic.find(std::pair<const G4Material*, co 399 if (dic.find(std::pair<const G4Material*, const G4Element*>((G4Material*)nullptr, element))
388 != dic.end()) 400 != dic.end())
389 return dic.find(std::pair<const G4Material 401 return dic.find(std::pair<const G4Material*, const G4Element*>((G4Material*)nullptr, element))
390 ->second; 402 ->second;
391 if (dic.find(std::pair<const G4Material*, co 403 if (dic.find(std::pair<const G4Material*, const G4Element*>(material, element)) != dic.end())
392 return dic.find(std::pair<const G4Material 404 return dic.find(std::pair<const G4Material*, const G4Element*>(material, element))->second;
393 return result; 405 return result;
394 } 406 }
395 407
396 G4double G4ParticleHPThermalScatteringData:: << 408 G4double G4ParticleHPThermalScatteringData::GetX(
397 GetX(const G4DynamicParticle* aP, G4double aT, << 409 const G4DynamicParticle* aP, G4double aT,
398 std::map<G4double, G4ParticleHPVector*>* << 410 std::map<G4double, G4ParticleHPVector*>* amapTemp_EnergyCross)
399 { 411 {
400 G4double result = 0; 412 G4double result = 0;
401 if (amapTemp_EnergyCross->empty()) return re 413 if (amapTemp_EnergyCross->empty()) return result;
402 414
403 G4double eKinetic = aP->GetKineticEnergy(); 415 G4double eKinetic = aP->GetKineticEnergy();
404 416
405 if (amapTemp_EnergyCross->size() == 1) { 417 if (amapTemp_EnergyCross->size() == 1) {
406 if (std::fabs(aT - amapTemp_EnergyCross->c 418 if (std::fabs(aT - amapTemp_EnergyCross->cbegin()->first) / amapTemp_EnergyCross->begin()->first
407 > 0.1) 419 > 0.1)
408 { 420 {
409 G4cout 421 G4cout
410 << "G4ParticleHPThermalScatteringData: 422 << "G4ParticleHPThermalScatteringData:: The temperature of material (" << aT / kelvin
411 << "K) is different more than 10% from 423 << "K) is different more than 10% from temperature of thermal scattering file expected ("
412 << amapTemp_EnergyCross->begin()->firs 424 << amapTemp_EnergyCross->begin()->first << "K). Result may not be reliable." << G4endl;
413 } 425 }
414 result = amapTemp_EnergyCross->begin()->se 426 result = amapTemp_EnergyCross->begin()->second->GetXsec(eKinetic);
415 return result; 427 return result;
416 } 428 }
417 429
418 auto it = amapTemp_EnergyCross->cbegin(); 430 auto it = amapTemp_EnergyCross->cbegin();
419 for (it = amapTemp_EnergyCross->cbegin(); it 431 for (it = amapTemp_EnergyCross->cbegin(); it != amapTemp_EnergyCross->cend(); ++it) {
420 if (aT < it->first) break; 432 if (aT < it->first) break;
421 } 433 }
422 if (it == amapTemp_EnergyCross->cbegin()) { 434 if (it == amapTemp_EnergyCross->cbegin()) {
423 ++it; // lower than the first 435 ++it; // lower than the first
424 } 436 }
425 else if (it == amapTemp_EnergyCross->cend()) 437 else if (it == amapTemp_EnergyCross->cend()) {
426 --it; // upper than the last 438 --it; // upper than the last
427 } 439 }
428 440
429 G4double TH = it->first; 441 G4double TH = it->first;
430 G4double XH = it->second->GetXsec(eKinetic); 442 G4double XH = it->second->GetXsec(eKinetic);
431 443
432 if (it != amapTemp_EnergyCross->cbegin()) -- 444 if (it != amapTemp_EnergyCross->cbegin()) --it;
433 G4double TL = it->first; 445 G4double TL = it->first;
434 G4double XL = it->second->GetXsec(eKinetic); 446 G4double XL = it->second->GetXsec(eKinetic);
435 447
436 if (TH == TL) throw G4HadronicException(__FI 448 if (TH == TL) throw G4HadronicException(__FILE__, __LINE__, "Thermal Scattering Data Error!");
437 449
438 G4double T = aT; 450 G4double T = aT;
439 G4double X = (XH - XL) / (TH - TL) * (T - TL 451 G4double X = (XH - XL) / (TH - TL) * (T - TL) + XL;
440 result = X; 452 result = X;
441 453
442 return result; 454 return result;
443 } 455 }
444 456
445 void G4ParticleHPThermalScatteringData::AddUse << 457 void G4ParticleHPThermalScatteringData::AddUserThermalScatteringFile(G4String nameG4Element,
446 << 458 G4String filename)
447 { 459 {
448 names->AddThermalElement(nameG4Element, file 460 names->AddThermalElement(nameG4Element, filename);
449 } 461 }
450 462
451 void G4ParticleHPThermalScatteringData::CrossS 463 void G4ParticleHPThermalScatteringData::CrossSectionDescription(std::ostream& outFile) const
452 { 464 {
453 outFile << "High Precision cross data based 465 outFile << "High Precision cross data based on thermal scattering data in evaluated nuclear data "
454 "libraries for neutrons below 5eV 466 "libraries for neutrons below 5eV on specific materials\n";
455 } 467 }
456 468