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