| 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 // Thermal Neutron Scattering
27 // Koi, Tatsumi (SLAC/SCCS)
28 //
29 // Class Description:
30 //
31 // Final State Generators for a high precision (based on evaluated data
32 // libraries) description of themal neutron scattering below 4 eV;
33 // Based on Thermal neutron scattering files
34 // from the evaluated nuclear data files ENDF/B-VI, Release2
35 // To be used in your physics list in case you need this physics.
36 // In this case you want to register an object of this class with
37 // the corresponding process.
38
39 // 070625 Fix memory leaking at destructor by T. Koi
40 // 081201 Fix memory leaking at destructor by T. Koi
41 // 100729 Add model name in constructor Problem #1116
42 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
43 //
44 #include "G4ParticleHPThermalScattering.hh"
45
46 #include "G4ElementTable.hh"
47 #include "G4MaterialTable.hh"
48 #include "G4Neutron.hh"
49 #include "G4ParticleHPElastic.hh"
50 #include "G4ParticleHPManager.hh"
51 #include "G4ParticleHPThermalScatteringData.hh"
52 #include "G4ParticleHPThermalScatteringNames.hh"
53 #include "G4SystemOfUnits.hh"
54 #include "G4Threading.hh"
55
56 G4ParticleHPThermalScattering::G4ParticleHPThermalScattering()
57 : G4HadronicInteraction("NeutronHPThermalScattering")
58 {
59 theHPElastic = new G4ParticleHPElastic();
60
61 SetMinEnergy(0. * eV);
62 SetMaxEnergy(4 * eV);
63 theXSection = new G4ParticleHPThermalScatteringData();
64
65 nMaterial = 0;
66 nElement = 0;
67 }
68
69 G4ParticleHPThermalScattering::~G4ParticleHPThermalScattering()
70 {
71 delete theHPElastic;
72 }
73
74 void G4ParticleHPThermalScattering::clearCurrentFSData()
75 {
76 if (incoherentFSs != nullptr) {
77 for (auto it = incoherentFSs->cbegin(); it != incoherentFSs->cend(); ++it) {
78 for (auto itt = it->second->cbegin(); itt != it->second->cend(); ++itt) {
79 for (auto ittt = itt->second->cbegin(); ittt != itt->second->cend(); ++ittt) {
80 delete *ittt;
81 }
82 delete itt->second;
83 }
84 delete it->second;
85 }
86 }
87
88 if (coherentFSs != nullptr) {
89 for (auto it = coherentFSs->cbegin(); it != coherentFSs->cend(); ++it) {
90 for (auto itt = it->second->cbegin(); itt != it->second->cend(); ++itt) {
91 for (auto ittt = itt->second->cbegin(); ittt != itt->second->cend(); ++ittt) {
92 delete *ittt;
93 }
94 delete itt->second;
95 }
96 delete it->second;
97 }
98 }
99
100 if (inelasticFSs != nullptr) {
101 for (auto it = inelasticFSs->cbegin(); it != inelasticFSs->cend(); ++it) {
102 for (auto itt = it->second->cbegin(); itt != it->second->cend(); ++itt) {
103 for (auto ittt = itt->second->cbegin(); ittt != itt->second->cend(); ++ittt) {
104 for (auto it4 = (*ittt)->vE_isoAngle.cbegin(); it4 != (*ittt)->vE_isoAngle.cend(); ++it4)
105 {
106 delete *it4;
107 }
108 delete *ittt;
109 }
110 delete itt->second;
111 }
112 delete it->second;
113 }
114 }
115
116 incoherentFSs = nullptr;
117 coherentFSs = nullptr;
118 inelasticFSs = nullptr;
119 }
120
121 void G4ParticleHPThermalScattering::BuildPhysicsTable(const G4ParticleDefinition& particle)
122 {
123 buildPhysicsTable();
124 theHPElastic->BuildPhysicsTable(particle);
125 }
126
127 std::map<G4double, std::vector<std::pair<G4double, G4double>*>*>*
128 G4ParticleHPThermalScattering::readACoherentFSDATA(const G4String& name)
129 {
130 auto aCoherentFSDATA = new std::map<G4double, std::vector<std::pair<G4double, G4double>*>*>;
131
132 std::istringstream theChannel(std::ios::in);
133 G4ParticleHPManager::GetInstance()->GetDataStream(name, theChannel);
134
135 std::vector<G4double> vBraggE;
136
137 G4int dummy;
138 while (theChannel >> dummy) // MF // Loop checking, 11.05.2015, T. Koi
139 {
140 theChannel >> dummy; // MT
141 G4double temp;
142 theChannel >> temp;
143 auto anBragE_P = new std::vector<std::pair<G4double, G4double>*>;
144
145 G4int n;
146 theChannel >> n;
147 for (G4int i = 0; i < n; ++i) {
148 G4double Ei;
149 G4double Pi;
150 if (aCoherentFSDATA->empty()) {
151 theChannel >> Ei;
152 vBraggE.push_back(Ei);
153 }
154 else {
155 Ei = vBraggE[i];
156 }
157 theChannel >> Pi;
158 anBragE_P->push_back(new std::pair<G4double, G4double>(Ei, Pi));
159 }
160 aCoherentFSDATA->insert(
161 std::pair<G4double, std::vector<std::pair<G4double, G4double>*>*>(temp, anBragE_P));
162 }
163 return aCoherentFSDATA;
164 }
165
166 std::map<G4double, std::vector<E_P_E_isoAng*>*>*
167 G4ParticleHPThermalScattering::readAnInelasticFSDATA(const G4String& name)
168 {
169 auto anT_E_P_E_isoAng = new std::map<G4double, std::vector<E_P_E_isoAng*>*>;
170
171 std::istringstream theChannel(std::ios::in);
172 G4ParticleHPManager::GetInstance()->GetDataStream(name, theChannel);
173
174 G4int dummy;
175 while (theChannel >> dummy) // MF // Loop checking, 11.05.2015, T. Koi
176 {
177 theChannel >> dummy; // MT
178 G4double temp;
179 theChannel >> temp;
180 auto vE_P_E_isoAng = new std::vector<E_P_E_isoAng*>;
181 G4int n;
182 theChannel >> n;
183 for (G4int i = 0; i < n; ++i) {
184 vE_P_E_isoAng->push_back(readAnE_P_E_isoAng(&theChannel));
185 }
186 anT_E_P_E_isoAng->insert(std::pair<G4double, std::vector<E_P_E_isoAng*>*>(temp, vE_P_E_isoAng));
187 }
188
189 return anT_E_P_E_isoAng;
190 }
191
192 E_P_E_isoAng*
193 G4ParticleHPThermalScattering::readAnE_P_E_isoAng(std::istream* file) // for inelastic
194 {
195 auto aData = new E_P_E_isoAng;
196
197 G4double dummy;
198 G4double energy;
199 G4int nep, nl;
200 *file >> dummy;
201 *file >> energy;
202 aData->energy = energy * eV;
203 *file >> dummy;
204 *file >> dummy;
205 *file >> nep;
206 *file >> nl;
207 aData->n = nep / nl;
208 for (G4int i = 0; i < aData->n; ++i) {
209 G4double prob;
210 auto anE_isoAng = new E_isoAng;
211 aData->vE_isoAngle.push_back(anE_isoAng);
212 *file >> energy;
213 anE_isoAng->energy = energy * eV;
214 anE_isoAng->n = nl - 2;
215 anE_isoAng->isoAngle.resize(anE_isoAng->n);
216 *file >> prob;
217 aData->prob.push_back(prob);
218 // G4cout << "G4ParticleHPThermalScattering inelastic " << energy/eV << " " << i << " " << prob
219 // << " " << aData->prob[ i ] << G4endl;
220 for (G4int j = 0; j < anE_isoAng->n; ++j) {
221 G4double x;
222 *file >> x;
223 anE_isoAng->isoAngle[j] = x;
224 }
225 }
226
227 // Calcuate sum_of_provXdEs
228 G4double total = 0;
229 aData->secondary_energy_cdf.push_back(0.);
230 for (G4int i = 0; i < aData->n - 1; ++i) {
231 G4double E_L = aData->vE_isoAngle[i]->energy / eV;
232 G4double E_H = aData->vE_isoAngle[i + 1]->energy / eV;
233 G4double dE = E_H - E_L;
234 G4double pdf = (aData->prob[i] + aData->prob[i + 1]) / 2. * dE;
235 total += (pdf);
236 aData->secondary_energy_cdf.push_back(total);
237 aData->secondary_energy_pdf.push_back(pdf);
238 aData->secondary_energy_value.push_back(E_L);
239 }
240
241 aData->sum_of_probXdEs = total;
242
243 // Normalize CDF
244 aData->secondary_energy_cdf_size = (G4int)aData->secondary_energy_cdf.size();
245 for (G4int i = 0; i < aData->secondary_energy_cdf_size; ++i) {
246 aData->secondary_energy_cdf[i] /= total;
247 }
248
249 return aData;
250 }
251
252 std::map<G4double, std::vector<E_isoAng*>*>*
253 G4ParticleHPThermalScattering::readAnIncoherentFSDATA(const G4String& name)
254 {
255 auto T_E = new std::map<G4double, std::vector<E_isoAng*>*>;
256
257 // std::ifstream theChannel( name.c_str() );
258 std::istringstream theChannel(std::ios::in);
259 G4ParticleHPManager::GetInstance()->GetDataStream(name, theChannel);
260
261 G4int dummy;
262 while (theChannel >> dummy) // MF // Loop checking, 11.05.2015, T. Koi
263 {
264 theChannel >> dummy; // MT
265 G4double temp;
266 theChannel >> temp;
267 auto vE_isoAng = new std::vector<E_isoAng*>;
268 G4int n;
269 theChannel >> n;
270 for (G4int i = 0; i < n; i++)
271 vE_isoAng->push_back(readAnE_isoAng(&theChannel));
272 T_E->insert(std::pair<G4double, std::vector<E_isoAng*>*>(temp, vE_isoAng));
273 }
274 // theChannel.close();
275
276 return T_E;
277 }
278
279 E_isoAng* G4ParticleHPThermalScattering::readAnE_isoAng(std::istream* file)
280 {
281 auto aData = new E_isoAng;
282
283 G4double dummy;
284 G4double energy;
285 G4int n;
286 *file >> dummy;
287 *file >> energy;
288 *file >> dummy;
289 *file >> dummy;
290 *file >> n;
291 *file >> dummy;
292 aData->energy = energy * eV;
293 aData->n = n - 2;
294 aData->isoAngle.resize(n);
295
296 *file >> dummy;
297 *file >> dummy;
298 for (G4int i = 0; i < aData->n; i++)
299 *file >> aData->isoAngle[i];
300
301 return aData;
302 }
303
304 G4HadFinalState* G4ParticleHPThermalScattering::ApplyYourself(const G4HadProjectile& aTrack,
305 G4Nucleus& aNucleus)
306 {
307 // Select Element > Reaction >
308
309 const G4Material* theMaterial = aTrack.GetMaterial();
310 G4double aTemp = theMaterial->GetTemperature();
311 auto n = (G4int)theMaterial->GetNumberOfElements();
312
313 G4bool findThermalElement = false;
314 G4int ielement;
315 const G4Element* theElement = nullptr;
316 for (G4int i = 0; i < n; ++i) {
317 theElement = theMaterial->GetElement(i);
318 // Select target element
319 if (aNucleus.GetZ_asInt() == (G4int)(theElement->GetZ() + 0.5)) {
320 // Check Applicability of Thermal Scattering
321 if (getTS_ID(nullptr, theElement) != -1) {
322 ielement = getTS_ID(nullptr, theElement);
323 findThermalElement = true;
324 break;
325 }
326 if (getTS_ID(theMaterial, theElement) != -1) {
327 ielement = getTS_ID(theMaterial, theElement);
328 findThermalElement = true;
329 break;
330 }
331 }
332 }
333
334 if (findThermalElement) {
335 // Select Reaction (Inelastic, coherent, incoherent)
336 const G4ParticleDefinition* pd = aTrack.GetDefinition();
337 auto dp = new G4DynamicParticle(pd, aTrack.Get4Momentum());
338 G4double total = theXSection->GetCrossSection(dp, theElement, theMaterial);
339 G4double inelastic = theXSection->GetInelasticCrossSection(dp, theElement, theMaterial);
340
341 G4double random = G4UniformRand();
342 if (random <= inelastic / total) {
343 // Inelastic
344
345 std::vector<G4double> v_temp;
346 v_temp.clear();
347 for (auto it = inelasticFSs->find(ielement)->second->cbegin();
348 it != inelasticFSs->find(ielement)->second->cend(); ++it)
349 {
350 v_temp.push_back(it->first);
351 }
352
353 std::pair<G4double, G4double> tempLH = find_LH(aTemp, &v_temp);
354 //
355 // For T_L aNEP_EPM_TL and T_H aNEP_EPM_TH
356 //
357 std::vector<E_P_E_isoAng*>* vNEP_EPM_TL = nullptr;
358 std::vector<E_P_E_isoAng*>* vNEP_EPM_TH = nullptr;
359
360 if (tempLH.first != 0.0 && tempLH.second != 0.0) {
361 vNEP_EPM_TL = inelasticFSs->find(ielement)->second->find(tempLH.first / kelvin)->second;
362 vNEP_EPM_TH = inelasticFSs->find(ielement)->second->find(tempLH.second / kelvin)->second;
363 }
364 else if (tempLH.first == 0.0) {
365 auto itm = inelasticFSs->find(ielement)->second->cbegin();
366 vNEP_EPM_TL = itm->second;
367 ++itm;
368 vNEP_EPM_TH = itm->second;
369 tempLH.first = tempLH.second;
370 tempLH.second = itm->first;
371 }
372 else if (tempLH.second == 0.0) {
373 auto itm = inelasticFSs->find(ielement)->second->cend();
374 --itm;
375 vNEP_EPM_TH = itm->second;
376 --itm;
377 vNEP_EPM_TL = itm->second;
378 tempLH.second = tempLH.first;
379 tempLH.first = itm->first;
380 }
381
382 G4double sE = 0., mu = 1.0;
383
384 // New Geant4 method - Stochastic temperature interpolation of the final state
385 // (continuous temperature interpolation was used previously)
386 std::pair<G4double, G4double> secondaryParam;
387 G4double rand_temp = G4UniformRand();
388 if (rand_temp < (aTemp - tempLH.first) / (tempLH.second - tempLH.first))
389 secondaryParam = sample_inelastic_E_mu(aTrack.GetKineticEnergy(), vNEP_EPM_TH);
390 else
391 secondaryParam = sample_inelastic_E_mu(aTrack.GetKineticEnergy(), vNEP_EPM_TL);
392
393 sE = secondaryParam.first;
394 mu = secondaryParam.second;
395
396 // set
397 theParticleChange.SetEnergyChange(sE);
398 G4double phi = CLHEP::twopi * G4UniformRand();
399 G4double sint = std::sqrt(1 - mu * mu);
400 theParticleChange.SetMomentumChange(sint * std::cos(phi), sint * std::sin(phi), mu);
401 }
402 else if (random
403 <= (inelastic + theXSection->GetCoherentCrossSection(dp, theElement, theMaterial))
404 / total)
405 {
406 // Coherent Elastic
407
408 G4double E = aTrack.GetKineticEnergy();
409
410 // T_L and T_H
411 std::vector<G4double> v_temp;
412 v_temp.clear();
413 for (auto it = coherentFSs->find(ielement)->second->cbegin();
414 it != coherentFSs->find(ielement)->second->cend(); ++it)
415 {
416 v_temp.push_back(it->first);
417 }
418
419 // T_L T_H
420 std::pair<G4double, G4double> tempLH = find_LH(aTemp, &v_temp);
421 //
422 //
423 // For T_L anEPM_TL and T_H anEPM_TH
424 //
425 std::vector<std::pair<G4double, G4double>*>* pvE_p_TL = nullptr;
426 std::vector<std::pair<G4double, G4double>*>* pvE_p_TH = nullptr;
427
428 if (tempLH.first != 0.0 && tempLH.second != 0.0) {
429 pvE_p_TL = coherentFSs->find(ielement)->second->find(tempLH.first / kelvin)->second;
430 pvE_p_TH = coherentFSs->find(ielement)->second->find(tempLH.first / kelvin)->second;
431 }
432 else if (tempLH.first == 0.0) {
433 pvE_p_TL = coherentFSs->find(ielement)->second->find(v_temp[0])->second;
434 pvE_p_TH = coherentFSs->find(ielement)->second->find(v_temp[1])->second;
435 tempLH.first = tempLH.second;
436 tempLH.second = v_temp[1];
437 }
438 else if (tempLH.second == 0.0) {
439 pvE_p_TH = coherentFSs->find(ielement)->second->find(v_temp.back())->second;
440 auto itv = v_temp.cend();
441 --itv;
442 --itv;
443 pvE_p_TL = coherentFSs->find(ielement)->second->find(*itv)->second;
444 tempLH.second = tempLH.first;
445 tempLH.first = *itv;
446 }
447 else {
448 // tempLH.first == 0.0 && tempLH.second
449 throw G4HadronicException(
450 __FILE__, __LINE__,
451 "A problem is found in Thermal Scattering Data! Unexpected temperature values in data");
452 }
453
454 std::vector<G4double> vE_T;
455 std::vector<G4double> vp_T;
456
457 auto n1 = (G4int)pvE_p_TL->size();
458
459 // New Geant4 method - Stochastic interpolation of the final state
460 std::vector<std::pair<G4double, G4double>*>* pvE_p_T_sampled;
461 G4double rand_temp = G4UniformRand();
462 if (rand_temp < (aTemp - tempLH.first) / (tempLH.second - tempLH.first))
463 pvE_p_T_sampled = pvE_p_TH;
464 else
465 pvE_p_T_sampled = pvE_p_TL;
466
467 // 171005 fix bug, contribution from H.N. TRAN@CEA
468 for (G4int i = 0; i < n1; ++i) {
469 vE_T.push_back((*pvE_p_T_sampled)[i]->first);
470 vp_T.push_back((*pvE_p_T_sampled)[i]->second);
471 }
472
473 G4int j = 0;
474 for (G4int i = 1; i < n1; ++i) {
475 if (E / eV < vE_T[i]) {
476 j = i - 1;
477 break;
478 }
479 }
480
481 G4double rand_for_mu = G4UniformRand();
482
483 G4int k = 0;
484 for (G4int i = 0; i <= j; ++i) {
485 G4double Pi = vp_T[i] / vp_T[j];
486 if (rand_for_mu < Pi) {
487 k = i;
488 break;
489 }
490 }
491
492 G4double Ei = vE_T[k];
493
494 G4double mu = 1 - 2 * Ei / (E / eV);
495
496 if (mu < -1.0) mu = -1.0;
497
498 theParticleChange.SetEnergyChange(E);
499 G4double phi = CLHEP::twopi * G4UniformRand();
500 G4double sint = std::sqrt(1 - mu * mu);
501 theParticleChange.SetMomentumChange(sint * std::cos(phi), sint * std::sin(phi), mu);
502 }
503 else {
504 // InCoherent Elastic
505
506 // T_L and T_H
507 std::vector<G4double> v_temp;
508 v_temp.clear();
509 for (auto it = incoherentFSs->find(ielement)->second->cbegin();
510 it != incoherentFSs->find(ielement)->second->cend(); ++it)
511 {
512 v_temp.push_back(it->first);
513 }
514
515 // T_L T_H
516 std::pair<G4double, G4double> tempLH = find_LH(aTemp, &v_temp);
517
518 //
519 // For T_L anEPM_TL and T_H anEPM_TH
520 //
521
522 E_isoAng anEPM_TL_E;
523 E_isoAng anEPM_TH_E;
524
525 if (tempLH.first != 0.0 && tempLH.second != 0.0) {
526 // Interpolate TL and TH
527 anEPM_TL_E = create_E_isoAng_from_energy(
528 aTrack.GetKineticEnergy(),
529 incoherentFSs->find(ielement)->second->find(tempLH.first / kelvin)->second);
530 anEPM_TH_E = create_E_isoAng_from_energy(
531 aTrack.GetKineticEnergy(),
532 incoherentFSs->find(ielement)->second->find(tempLH.second / kelvin)->second);
533 }
534 else if (tempLH.first == 0.0) {
535 // Extrapolate T0 and T1
536 anEPM_TL_E = create_E_isoAng_from_energy(
537 aTrack.GetKineticEnergy(),
538 incoherentFSs->find(ielement)->second->find(v_temp[0])->second);
539 anEPM_TH_E = create_E_isoAng_from_energy(
540 aTrack.GetKineticEnergy(),
541 incoherentFSs->find(ielement)->second->find(v_temp[1])->second);
542 tempLH.first = tempLH.second;
543 tempLH.second = v_temp[1];
544 }
545 else if (tempLH.second == 0.0) {
546 // Extrapolate Tmax-1 and Tmax
547 anEPM_TH_E = create_E_isoAng_from_energy(
548 aTrack.GetKineticEnergy(),
549 incoherentFSs->find(ielement)->second->find(v_temp.back())->second);
550 auto itv = v_temp.cend();
551 --itv;
552 --itv;
553 anEPM_TL_E = create_E_isoAng_from_energy(
554 aTrack.GetKineticEnergy(), incoherentFSs->find(ielement)->second->find(*itv)->second);
555 tempLH.second = tempLH.first;
556 tempLH.first = *itv;
557 }
558
559 // E_isoAng for aTemp and aTrack.GetKineticEnergy()
560 G4double mu = 1.0;
561
562 // New Geant4 method - Stochastic interpolation of the final state
563 E_isoAng anEPM_T_E_sampled;
564 G4double rand_temp = G4UniformRand();
565 if (rand_temp < (aTemp - tempLH.first) / (tempLH.second - tempLH.first))
566 anEPM_T_E_sampled = std::move(anEPM_TH_E);
567 else
568 anEPM_T_E_sampled = std::move(anEPM_TL_E);
569
570 mu = getMu(&anEPM_T_E_sampled);
571
572 // Set Final State
573 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy()); // No energy change in Elastic
574 G4double phi = CLHEP::twopi * G4UniformRand();
575 G4double sint = std::sqrt(1 - mu * mu);
576 theParticleChange.SetMomentumChange(sint * std::cos(phi), sint * std::sin(phi), mu);
577 }
578 delete dp;
579
580 return &theParticleChange;
581 }
582 // Not thermal element
583 // Neutron HP will handle
584 return theHPElastic->ApplyYourself(aTrack, aNucleus,
585 true); // L. Thulliez 2021/05/04 (CEA-Saclay)
586 }
587
588 //**********************************************************
589 // Geant4 new algorithm
590 //**********************************************************
591
592 //--------------------------------------------------
593 // New method added by L. Thulliez 2021 (CEA-Saclay)
594 //--------------------------------------------------
595 std::pair<G4double, G4int>
596 G4ParticleHPThermalScattering::sample_inelastic_E(G4double rndm1, G4double rndm2,
597 E_P_E_isoAng* anE_P_E_isoAng)
598 {
599 G4int i = 0;
600 G4double sE_value = 0;
601
602 for (; i < anE_P_E_isoAng->secondary_energy_cdf_size - 1; ++i) {
603 if (rndm1 >= anE_P_E_isoAng->secondary_energy_cdf[i]
604 && rndm1 < anE_P_E_isoAng->secondary_energy_cdf[i + 1])
605 {
606 G4double sE_value_i = anE_P_E_isoAng->secondary_energy_value[i];
607 G4double sE_pdf_i = anE_P_E_isoAng->secondary_energy_pdf[i];
608 G4double sE_value_i1 = anE_P_E_isoAng->secondary_energy_value[i + 1];
609 G4double sE_pdf_i1 = anE_P_E_isoAng->secondary_energy_pdf[i + 1];
610
611 G4double lambda = 0;
612 G4double alpha = (sE_pdf_i1 - sE_pdf_i) / (sE_pdf_i1 + sE_pdf_i);
613 G4double rndm = rndm1;
614
615 if (std::fabs(alpha) < 1E-8) {
616 lambda = rndm2;
617 }
618 else {
619 G4double beta = 2 * sE_pdf_i / (sE_pdf_i1 + sE_pdf_i);
620 rndm = rndm2;
621 G4double gamma = -rndm;
622 G4double delta = beta * beta - 4 * alpha * gamma;
623
624 if (delta < 0 && std::fabs(delta) < 1.E-8) delta = 0;
625
626 lambda = -beta + std::sqrt(delta);
627 lambda = lambda / (2 * alpha);
628
629 if (lambda > 1)
630 lambda = 1;
631 else if (lambda < 0)
632 lambda = 0;
633 }
634
635 sE_value = sE_value_i + lambda * (sE_value_i1 - sE_value_i);
636
637 break;
638 }
639 }
640
641 return std::pair<G4double, G4int>(sE_value, i);
642 }
643
644 //--------------------------------------------------
645 // New method added by L. Thulliez 2021 (CEA-Saclay)
646 //--------------------------------------------------
647 std::pair<G4double, G4double>
648 G4ParticleHPThermalScattering::sample_inelastic_E_mu(G4double pE,
649 std::vector<E_P_E_isoAng*>* vNEP_EPM)
650 {
651 // Sample primary energy bin
652 std::map<G4double, G4int> map_energy;
653 map_energy.clear();
654 std::vector<G4double> v_energy;
655 v_energy.clear();
656 G4int i = 0;
657 for (auto itv = vNEP_EPM->cbegin(); itv != vNEP_EPM->cend(); ++itv) {
658 v_energy.push_back((*itv)->energy);
659 map_energy.insert(std::pair<G4double, G4int>((*itv)->energy, i));
660 i++;
661 }
662
663 std::pair<G4double, G4double> energyLH = find_LH(pE, &v_energy);
664
665 std::vector<E_P_E_isoAng*> pE_P_E_isoAng_limit(2, nullptr);
666
667 if (energyLH.first != 0.0 && energyLH.second != 0.0) {
668 pE_P_E_isoAng_limit[0] = (*vNEP_EPM)[map_energy.find(energyLH.first)->second];
669 pE_P_E_isoAng_limit[1] = (*vNEP_EPM)[map_energy.find(energyLH.second)->second];
670 }
671 else if (energyLH.first == 0.0) {
672 pE_P_E_isoAng_limit[0] = (*vNEP_EPM)[0];
673 pE_P_E_isoAng_limit[1] = (*vNEP_EPM)[1];
674 }
675 if (energyLH.second == 0.0) {
676 pE_P_E_isoAng_limit[1] = (*vNEP_EPM).back();
677 auto itv = vNEP_EPM->cend();
678 --itv;
679 --itv;
680 pE_P_E_isoAng_limit[0] = *itv;
681 }
682
683 // Compute interpolation factor of the incident neutron energy
684 G4double factor = (energyLH.second - pE) / (energyLH.second - energyLH.first);
685
686 if ((energyLH.second - pE) <= 0. && std::fabs(pE / energyLH.second - 1) < 1E-11) factor = 0.;
687 if ((energyLH.first - pE) >= 0. && std::fabs(energyLH.first / pE - 1) < 1E-11) factor = 1.;
688
689 G4double rndm1 = G4UniformRand();
690 G4double rndm2 = G4UniformRand();
691
692 // Sample secondary neutron energy
693 std::pair<G4double, G4int> sE_lower = sample_inelastic_E(rndm1, rndm2, pE_P_E_isoAng_limit[0]);
694 std::pair<G4double, G4int> sE_upper = sample_inelastic_E(rndm1, rndm2, pE_P_E_isoAng_limit[1]);
695 G4double sE = factor * sE_lower.first + (1 - factor) * sE_upper.first;
696 sE = sE * eV;
697
698 // Sample cosine knowing the secondary neutron energy
699 rndm1 = G4UniformRand();
700 rndm2 = G4UniformRand();
701 G4double mu_lower = getMu(rndm1, rndm2, pE_P_E_isoAng_limit[0]->vE_isoAngle[sE_lower.second]);
702 G4double mu_upper = getMu(rndm1, rndm2, pE_P_E_isoAng_limit[1]->vE_isoAngle[sE_upper.second]);
703 G4double mu = factor * mu_lower + (1 - factor) * mu_upper;
704
705 return std::pair<G4double, G4double>(sE, mu);
706 }
707
708 //--------------------------------------------------
709 // New method added by L. Thulliez 2021 (CEA-Saclay)
710 //--------------------------------------------------
711 G4double G4ParticleHPThermalScattering::getMu(G4double rndm1, G4double rndm2, E_isoAng* anEPM)
712 {
713 G4double result = 0.0;
714
715 auto in = G4int(rndm1 * ((*anEPM).n));
716
717 if (in != 0) {
718 G4double mu_l = (*anEPM).isoAngle[in - 1];
719 G4double mu_h = (*anEPM).isoAngle[in];
720 result = (mu_h - mu_l) * (rndm1 * ((*anEPM).n) - in) + mu_l;
721 }
722 else {
723 G4double x = rndm1 * (*anEPM).n;
724 G4double ratio = 0.5;
725 if (x <= ratio) {
726 G4double mu_l = -1;
727 G4double mu_h = (*anEPM).isoAngle[0];
728 result = (mu_h - mu_l) * rndm2 + mu_l;
729 }
730 else {
731 G4double mu_l = (*anEPM).isoAngle[(*anEPM).n - 1];
732 G4double mu_h = 1;
733 result = (mu_h - mu_l) * rndm2 + mu_l;
734 }
735 }
736
737 return result;
738 }
739
740 //**********************************************************
741 // Geant4 previous algorithm
742 //**********************************************************
743
744 G4double G4ParticleHPThermalScattering::getMu(E_isoAng* anEPM)
745 {
746 G4double random = G4UniformRand();
747 G4double result = 0.0;
748
749 auto in = G4int(random * ((*anEPM).n));
750
751 if (in != 0) {
752 G4double mu_l = (*anEPM).isoAngle[in - 1];
753 G4double mu_h = (*anEPM).isoAngle[in];
754 result = (mu_h - mu_l) * (random * ((*anEPM).n) - in) + mu_l;
755 }
756 else {
757 G4double x = random * (*anEPM).n;
758 // Bugzilla 1971
759 G4double ratio = 0.5;
760 G4double xx = G4UniformRand();
761 if (x <= ratio) {
762 G4double mu_l = -1;
763 G4double mu_h = (*anEPM).isoAngle[0];
764 result = (mu_h - mu_l) * xx + mu_l;
765 }
766 else {
767 G4double mu_l = (*anEPM).isoAngle[(*anEPM).n - 1];
768 G4double mu_h = 1;
769 result = (mu_h - mu_l) * xx + mu_l;
770 }
771 }
772 return result;
773 }
774
775 std::pair<G4double, G4double> G4ParticleHPThermalScattering::find_LH(G4double x,
776 std::vector<G4double>* aVector)
777 {
778 G4double LL = 0.0;
779 G4double H = 0.0;
780
781 // v->size() == 1 --> LL=H=v(0)
782 if (aVector->size() == 1) {
783 LL = aVector->front();
784 H = aVector->front();
785 }
786 else {
787 // 1) temp < v(0) -> LL=0.0 H=v(0)
788 // 2) v(i-1) < temp <= v(i) -> LL=v(i-1) H=v(i)
789 // 3) v(imax) < temp -> LL=v(imax) H=0.0
790 for (auto it = aVector->cbegin(); it != aVector->cend(); ++it) {
791 if (x <= *it) {
792 H = *it;
793 if (it != aVector->cbegin()) {
794 // 2)
795 it--;
796 LL = *it;
797 }
798 else {
799 // 1)
800 LL = 0.0;
801 }
802 break;
803 }
804 }
805 // 3)
806 if (H == 0.0) LL = aVector->back();
807 }
808
809 return std::pair<G4double, G4double>(LL, H);
810 }
811
812 G4double G4ParticleHPThermalScattering::get_linear_interpolated(G4double x,
813 std::pair<G4double, G4double> Low,
814 std::pair<G4double, G4double> High)
815 {
816 G4double y = 0.0;
817 if (High.first - Low.first != 0) {
818 y = (High.second - Low.second) / (High.first - Low.first) * (x - Low.first) + Low.second;
819 }
820 else {
821 if (High.second == Low.second) {
822 y = High.second;
823 }
824 else {
825 G4cout << "G4ParticleHPThermalScattering liner interpolation err!!" << G4endl;
826 }
827 }
828
829 return y;
830 }
831
832 E_isoAng G4ParticleHPThermalScattering::create_E_isoAng_from_energy(G4double energy,
833 std::vector<E_isoAng*>* vEPM)
834 {
835 E_isoAng anEPM_T_E;
836
837 std::vector<G4double> v_e;
838 v_e.clear();
839 for (auto iv = vEPM->cbegin(); iv != vEPM->cend(); ++iv)
840 v_e.push_back((*iv)->energy);
841
842 std::pair<G4double, G4double> energyLH = find_LH(energy, &v_e);
843 // G4cout << " " << energy/eV << " " << energyLH.first/eV << " " << energyLH.second/eV << G4endl;
844
845 E_isoAng* panEPM_T_EL = nullptr;
846 E_isoAng* panEPM_T_EH = nullptr;
847
848 if (energyLH.first != 0.0 && energyLH.second != 0.0) {
849 for (auto iv = vEPM->cbegin(); iv != vEPM->cend(); ++iv) {
850 if (energyLH.first == (*iv)->energy) {
851 panEPM_T_EL = *iv;
852 ++iv;
853 panEPM_T_EH = *iv;
854 break;
855 }
856 }
857 }
858 else if (energyLH.first == 0.0) {
859 panEPM_T_EL = (*vEPM)[0];
860 panEPM_T_EH = (*vEPM)[1];
861 }
862 else if (energyLH.second == 0.0) {
863 panEPM_T_EH = (*vEPM).back();
864 auto iv = vEPM->cend();
865 --iv;
866 --iv;
867 panEPM_T_EL = *iv;
868 }
869
870 if (panEPM_T_EL != nullptr && panEPM_T_EH != nullptr) {
871 // checking isoAng has proper values or not
872 // Inelastic/FS, the first and last entries of *vEPM has all zero values.
873 if (!(check_E_isoAng(panEPM_T_EL))) panEPM_T_EL = panEPM_T_EH;
874 if (!(check_E_isoAng(panEPM_T_EH))) panEPM_T_EH = panEPM_T_EL;
875
876 if (panEPM_T_EL->n == panEPM_T_EH->n) {
877 anEPM_T_E.energy = energy;
878 anEPM_T_E.n = panEPM_T_EL->n;
879
880 for (G4int i = 0; i < panEPM_T_EL->n; ++i) {
881 G4double angle;
882 angle = get_linear_interpolated(
883 energy, std::pair<G4double, G4double>(energyLH.first, panEPM_T_EL->isoAngle[i]),
884 std::pair<G4double, G4double>(energyLH.second, panEPM_T_EH->isoAngle[i]));
885 anEPM_T_E.isoAngle.push_back(angle);
886 }
887 }
888 else {
889 G4Exception("G4ParticleHPThermalScattering::create_E_isoAng_from_energy", "NotSupported",
890 JustWarning,
891 "G4ParticleHPThermalScattering does not support yet EL->n != EH->n.");
892 }
893 }
894 else {
895 G4Exception("G4ParticleHPThermalScattering::create_E_isoAng_from_energy", "HAD_THERM_000",
896 FatalException, "Pointer panEPM_T_EL or panEPM_T_EH is zero");
897 }
898
899 return anEPM_T_E;
900 }
901
902 G4double
903 G4ParticleHPThermalScattering::get_secondary_energy_from_E_P_E_isoAng(G4double random,
904 E_P_E_isoAng* anE_P_E_isoAng)
905 {
906 G4double secondary_energy = 0.0;
907
908 G4int n = anE_P_E_isoAng->n;
909 G4double sum_p = 0.0; // sum_p_H
910 G4double sum_p_L = 0.0;
911
912 G4double total = 0.0;
913
914 /*
915 delete for speed up
916 for ( G4int i = 0 ; i < n-1 ; ++i )
917 {
918 G4double E_L = anE_P_E_isoAng->vE_isoAngle[i]->energy/eV;
919 G4double E_H = anE_P_E_isoAng->vE_isoAngle[i+1]->energy/eV;
920 G4double dE = E_H - E_L;
921 total += ( ( anE_P_E_isoAng->prob[i] ) * dE );
922 }
923
924 if ( std::abs( total - anE_P_E_isoAng->sum_of_probXdEs ) > 1.0e-14 ) G4cout << total -
925 anE_P_E_isoAng->sum_of_probXdEs << G4endl;
926 */
927 total = anE_P_E_isoAng->sum_of_probXdEs;
928
929 for (G4int i = 0; i < n - 1; ++i) {
930 G4double E_L = anE_P_E_isoAng->vE_isoAngle[i]->energy / eV;
931 G4double E_H = anE_P_E_isoAng->vE_isoAngle[i + 1]->energy / eV;
932 G4double dE = E_H - E_L;
933 sum_p += ((anE_P_E_isoAng->prob[i]) * dE);
934
935 if (random <= sum_p / total) {
936 secondary_energy =
937 get_linear_interpolated(random, std::pair<G4double, G4double>(sum_p_L / total, E_L),
938 std::pair<G4double, G4double>(sum_p / total, E_H));
939 secondary_energy = secondary_energy * eV; // need eV
940 break;
941 }
942 sum_p_L = sum_p;
943 }
944
945 return secondary_energy;
946 }
947
948 std::pair<G4double, E_isoAng>
949 G4ParticleHPThermalScattering::create_sE_and_EPM_from_pE_and_vE_P_E_isoAng(
950 G4double rand_for_sE, G4double pE, std::vector<E_P_E_isoAng*>* vNEP_EPM)
951 {
952 std::map<G4double, G4int> map_energy;
953 map_energy.clear();
954 std::vector<G4double> v_energy;
955 v_energy.clear();
956 G4int i = 0;
957 for (auto itv = vNEP_EPM->cbegin(); itv != vNEP_EPM->cend(); ++itv) {
958 v_energy.push_back((*itv)->energy);
959 map_energy.insert(std::pair<G4double, G4int>((*itv)->energy, i));
960 i++;
961 }
962
963 std::pair<G4double, G4double> energyLH = find_LH(pE, &v_energy);
964
965 E_P_E_isoAng* pE_P_E_isoAng_EL = nullptr;
966 E_P_E_isoAng* pE_P_E_isoAng_EH = nullptr;
967
968 if (energyLH.first != 0.0 && energyLH.second != 0.0) {
969 pE_P_E_isoAng_EL = (*vNEP_EPM)[map_energy.find(energyLH.first)->second];
970 pE_P_E_isoAng_EH = (*vNEP_EPM)[map_energy.find(energyLH.second)->second];
971 }
972 else if (energyLH.first == 0.0) {
973 pE_P_E_isoAng_EL = (*vNEP_EPM)[0];
974 pE_P_E_isoAng_EH = (*vNEP_EPM)[1];
975 }
976 if (energyLH.second == 0.0) {
977 pE_P_E_isoAng_EH = (*vNEP_EPM).back();
978 auto itv = vNEP_EPM->cend();
979 --itv;
980 --itv;
981 pE_P_E_isoAng_EL = *itv;
982 }
983
984 G4double sE;
985 G4double sE_L;
986 G4double sE_H;
987
988 sE_L = get_secondary_energy_from_E_P_E_isoAng(rand_for_sE, pE_P_E_isoAng_EL);
989 sE_H = get_secondary_energy_from_E_P_E_isoAng(rand_for_sE, pE_P_E_isoAng_EH);
990
991 sE = get_linear_interpolated(pE, std::pair<G4double, G4double>(energyLH.first, sE_L),
992 std::pair<G4double, G4double>(energyLH.second, sE_H));
993
994 E_isoAng E_isoAng_L = create_E_isoAng_from_energy(sE, &(pE_P_E_isoAng_EL->vE_isoAngle));
995 E_isoAng E_isoAng_H = create_E_isoAng_from_energy(sE, &(pE_P_E_isoAng_EH->vE_isoAngle));
996
997 E_isoAng anE_isoAng;
998 // For defeating warning message from compiler
999 anE_isoAng.n = 1;
1000 anE_isoAng.energy = sE; // never used
1001 if (E_isoAng_L.n == E_isoAng_H.n) {
1002 anE_isoAng.n = E_isoAng_L.n;
1003 for (G4int j = 0; j < anE_isoAng.n; ++j) {
1004 G4double angle;
1005 angle =
1006 get_linear_interpolated(sE, std::pair<G4double, G4double>(sE_L, E_isoAng_L.isoAngle[j]),
1007 std::pair<G4double, G4double>(sE_H, E_isoAng_H.isoAngle[j]));
1008 anE_isoAng.isoAngle.push_back(angle);
1009 }
1010 }
1011 else {
1012 throw G4HadronicException(__FILE__, __LINE__, "Unexpected values!");
1013 }
1014
1015 return std::pair<G4double, E_isoAng>(sE, anE_isoAng);
1016 }
1017
1018 void G4ParticleHPThermalScattering::buildPhysicsTable()
1019 {
1020 // Is rebuild of physics table a necessity
1021 if (nMaterial == G4Material::GetMaterialTable()->size()
1022 && nElement == G4Element::GetElementTable()->size())
1023 {
1024 return;
1025 }
1026 nMaterial = G4Material::GetMaterialTable()->size();
1027 nElement = G4Element::GetElementTable()->size();
1028
1029 dic.clear();
1030 std::map<G4String, G4int> co_dic;
1031
1032 // Searching Nist Materials
1033 static G4ThreadLocal G4MaterialTable* theMaterialTable = nullptr;
1034 if (theMaterialTable == nullptr) theMaterialTable = G4Material::GetMaterialTable();
1035 std::size_t numberOfMaterials = G4Material::GetNumberOfMaterials();
1036 for (std::size_t i = 0; i < numberOfMaterials; ++i) {
1037 G4Material* material = (*theMaterialTable)[i];
1038 auto numberOfElements = (G4int)material->GetNumberOfElements();
1039 for (G4int j = 0; j < numberOfElements; ++j) {
1040 const G4Element* element = material->GetElement(j);
1041 if (names.IsThisThermalElement(material->GetName(), element->GetName())) {
1042 G4int ts_ID_of_this_geometry;
1043 G4String ts_ndl_name = names.GetTS_NDL_Name(material->GetName(), element->GetName());
1044 if (co_dic.find(ts_ndl_name) != co_dic.cend()) {
1045 ts_ID_of_this_geometry = co_dic.find(ts_ndl_name)->second;
1046 }
1047 else {
1048 ts_ID_of_this_geometry = (G4int)co_dic.size();
1049 co_dic.insert(std::pair<G4String, G4int>(ts_ndl_name, ts_ID_of_this_geometry));
1050 }
1051
1052 // G4cout << "Neutron HP Thermal Scattering: Registering a material-element pair of "
1053 // << material->GetName() << " " << element->GetName()
1054 // << " as internal thermal scattering id of " << ts_ID_of_this_geometry << "." <<
1055 // G4endl;
1056
1057 dic.insert(std::pair<std::pair<G4Material*, const G4Element*>, G4int>(
1058 std::pair<G4Material*, const G4Element*>(material, element), ts_ID_of_this_geometry));
1059 }
1060 }
1061 }
1062
1063 // Searching TS Elements
1064 auto theElementTable = G4Element::GetElementTable();
1065 std::size_t numberOfElements = G4Element::GetNumberOfElements();
1066 for (std::size_t i = 0; i < numberOfElements; ++i) {
1067 const G4Element* element = (*theElementTable)[i];
1068 if (names.IsThisThermalElement(element->GetName())) {
1069 if (names.IsThisThermalElement(element->GetName())) {
1070 G4int ts_ID_of_this_geometry;
1071 G4String ts_ndl_name = names.GetTS_NDL_Name(element->GetName());
1072 if (co_dic.find(ts_ndl_name) != co_dic.cend()) {
1073 ts_ID_of_this_geometry = co_dic.find(ts_ndl_name)->second;
1074 }
1075 else {
1076 ts_ID_of_this_geometry = (G4int)co_dic.size();
1077 co_dic.insert(std::pair<G4String, G4int>(ts_ndl_name, ts_ID_of_this_geometry));
1078 }
1079 dic.insert(std::pair<std::pair<const G4Material*, const G4Element*>, G4int>(
1080 std::pair<const G4Material*, const G4Element*>((G4Material*)nullptr, element),
1081 ts_ID_of_this_geometry));
1082 }
1083 }
1084 }
1085
1086 G4cout << G4endl;
1087 G4cout
1088 << "Neutron HP Thermal Scattering: Following material-element pairs or elements are registered."
1089 << G4endl;
1090 for (const auto& it : dic) {
1091 if (it.first.first != nullptr) {
1092 G4cout << "Material " << it.first.first->GetName() << " - Element "
1093 << it.first.second->GetName() << ", internal thermal scattering id " << it.second
1094 << G4endl;
1095 }
1096 else {
1097 G4cout << "Element " << it.first.second->GetName() << ", internal thermal scattering id "
1098 << it.second << G4endl;
1099 }
1100 }
1101 G4cout << G4endl;
1102
1103 // Read Cross Section Data files
1104
1105 G4ParticleHPManager* hpmanager = G4ParticleHPManager::GetInstance();
1106 coherentFSs = hpmanager->GetThermalScatteringCoherentFinalStates();
1107 incoherentFSs = hpmanager->GetThermalScatteringIncoherentFinalStates();
1108 inelasticFSs = hpmanager->GetThermalScatteringInelasticFinalStates();
1109
1110 if (G4Threading::IsMasterThread()) {
1111 clearCurrentFSData();
1112
1113 if (coherentFSs == nullptr)
1114 coherentFSs =
1115 new std::map<G4int, std::map<G4double, std::vector<std::pair<G4double, G4double>*>*>*>;
1116 if (incoherentFSs == nullptr)
1117 incoherentFSs = new std::map<G4int, std::map<G4double, std::vector<E_isoAng*>*>*>;
1118 if (inelasticFSs == nullptr)
1119 inelasticFSs = new std::map<G4int, std::map<G4double, std::vector<E_P_E_isoAng*>*>*>;
1120
1121 G4String dirName;
1122 if (G4FindDataDir("G4NEUTRONHPDATA") == nullptr)
1123 throw G4HadronicException(
1124 __FILE__, __LINE__,
1125 "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
1126 dirName = G4FindDataDir("G4NEUTRONHPDATA");
1127
1128 for (const auto& it : co_dic) {
1129 G4String tsndlName = it.first;
1130 G4int ts_ID = it.second;
1131
1132 // Coherent
1133 G4String fsName = "/ThermalScattering/Coherent/FS/";
1134 G4String fileName = dirName + fsName + tsndlName;
1135 coherentFSs->insert(
1136 std::pair<G4int, std::map<G4double, std::vector<std::pair<G4double, G4double>*>*>*>(
1137 ts_ID, readACoherentFSDATA(fileName)));
1138
1139 // incoherent elastic
1140 fsName = "/ThermalScattering/Incoherent/FS/";
1141 fileName = dirName + fsName + tsndlName;
1142 incoherentFSs->insert(std::pair<G4int, std::map<G4double, std::vector<E_isoAng*>*>*>(
1143 ts_ID, readAnIncoherentFSDATA(fileName)));
1144
1145 // inelastic
1146 fsName = "/ThermalScattering/Inelastic/FS/";
1147 fileName = dirName + fsName + tsndlName;
1148 inelasticFSs->insert(std::pair<G4int, std::map<G4double, std::vector<E_P_E_isoAng*>*>*>(
1149 ts_ID, readAnInelasticFSDATA(fileName)));
1150 }
1151
1152 hpmanager->RegisterThermalScatteringCoherentFinalStates(coherentFSs);
1153 hpmanager->RegisterThermalScatteringIncoherentFinalStates(incoherentFSs);
1154 hpmanager->RegisterThermalScatteringInelasticFinalStates(inelasticFSs);
1155 }
1156
1157 theXSection->BuildPhysicsTable(*(G4Neutron::Neutron()));
1158 }
1159
1160 G4int G4ParticleHPThermalScattering::getTS_ID(const G4Material* material, const G4Element* element)
1161 {
1162 G4int result = -1;
1163 if (dic.find(std::pair<const G4Material*, const G4Element*>(material, element)) != dic.end())
1164 result = dic.find(std::pair<const G4Material*, const G4Element*>(material, element))->second;
1165 return result;
1166 }
1167
1168 const std::pair<G4double, G4double> G4ParticleHPThermalScattering::GetFatalEnergyCheckLevels() const
1169 {
1170 // max energy non-conservation is mass of heavy nucleus
1171 return std::pair<G4double, G4double>(10.0 * perCent, 350.0 * CLHEP::GeV);
1172 }
1173
1174 void G4ParticleHPThermalScattering::AddUserThermalScatteringFile(const G4String& nameG4Element,
1175 const G4String& filename)
1176 {
1177 names.AddThermalElement(nameG4Element, filename);
1178 theXSection->AddUserThermalScatteringFile(nameG4Element, filename);
1179 buildPhysicsTable();
1180 }
1181
1182 G4bool G4ParticleHPThermalScattering::check_E_isoAng(E_isoAng* anE_IsoAng)
1183 {
1184 G4bool result = false;
1185
1186 G4int n = anE_IsoAng->n;
1187 G4double sum = 0.0;
1188 for (G4int i = 0; i < n; ++i) {
1189 sum += anE_IsoAng->isoAngle[i];
1190 }
1191 if (sum != 0.0) result = true;
1192
1193 return result;
1194 }
1195
1196 void G4ParticleHPThermalScattering::ModelDescription(std::ostream& outFile) const
1197 {
1198 outFile << "High Precision model based on thermal scattering data in\n"
1199 << "evaluated nuclear data libraries for neutrons below 5eV\n"
1200 << "on specific materials\n";
1201 }
1202