| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Thermal Neutron Scattering
27 // Koi, Tatsumi (SLAC/SCCS)
28 //
29 // Class Description:
30 //
31 // Final State Generators for a high precision
32 // libraries) description of themal neutron sc
33 // Based on Thermal neutron scattering files
34 // from the evaluated nuclear data files ENDF/
35 // To be used in your physics list in case you
36 // In this case you want to register an object
37 // the corresponding process.
38
39 // 070625 Fix memory leaking at destructor by
40 // 081201 Fix memory leaking at destructor by
41 // 100729 Add model name in constructor Proble
42 // P. Arce, June-2014 Conversion neutron_hp to
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.h
53 #include "G4SystemOfUnits.hh"
54 #include "G4Threading.hh"
55
56 G4ParticleHPThermalScattering::G4ParticleHPThe
57 : G4HadronicInteraction("NeutronHPThermalSca
58 {
59 theHPElastic = new G4ParticleHPElastic();
60
61 SetMinEnergy(0. * eV);
62 SetMaxEnergy(4 * eV);
63 theXSection = new G4ParticleHPThermalScatter
64
65 nMaterial = 0;
66 nElement = 0;
67 }
68
69 G4ParticleHPThermalScattering::~G4ParticleHPTh
70 {
71 delete theHPElastic;
72 }
73
74 void G4ParticleHPThermalScattering::clearCurre
75 {
76 if (incoherentFSs != nullptr) {
77 for (auto it = incoherentFSs->cbegin(); it
78 for (auto itt = it->second->cbegin(); it
79 for (auto ittt = itt->second->cbegin()
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 !
90 for (auto itt = it->second->cbegin(); it
91 for (auto ittt = itt->second->cbegin()
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
102 for (auto itt = it->second->cbegin(); it
103 for (auto ittt = itt->second->cbegin()
104 for (auto it4 = (*ittt)->vE_isoAngle
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::BuildPhysi
122 {
123 buildPhysicsTable();
124 theHPElastic->BuildPhysicsTable(particle);
125 }
126
127 std::map<G4double, std::vector<std::pair<G4dou
128 G4ParticleHPThermalScattering::readACoherentFS
129 {
130 auto aCoherentFSDATA = new std::map<G4double
131
132 std::istringstream theChannel(std::ios::in);
133 G4ParticleHPManager::GetInstance()->GetDataS
134
135 std::vector<G4double> vBraggE;
136
137 G4int dummy;
138 while (theChannel >> dummy) // MF // Loop c
139 {
140 theChannel >> dummy; // MT
141 G4double temp;
142 theChannel >> temp;
143 auto anBragE_P = new std::vector<std::pair
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<G4dou
159 }
160 aCoherentFSDATA->insert(
161 std::pair<G4double, std::vector<std::pai
162 }
163 return aCoherentFSDATA;
164 }
165
166 std::map<G4double, std::vector<E_P_E_isoAng*>*
167 G4ParticleHPThermalScattering::readAnInelastic
168 {
169 auto anT_E_P_E_isoAng = new std::map<G4doubl
170
171 std::istringstream theChannel(std::ios::in);
172 G4ParticleHPManager::GetInstance()->GetDataS
173
174 G4int dummy;
175 while (theChannel >> dummy) // MF // Loop c
176 {
177 theChannel >> dummy; // MT
178 G4double temp;
179 theChannel >> temp;
180 auto vE_P_E_isoAng = new std::vector<E_P_E
181 G4int n;
182 theChannel >> n;
183 for (G4int i = 0; i < n; ++i) {
184 vE_P_E_isoAng->push_back(readAnE_P_E_iso
185 }
186 anT_E_P_E_isoAng->insert(std::pair<G4doubl
187 }
188
189 return anT_E_P_E_isoAng;
190 }
191
192 E_P_E_isoAng*
193 G4ParticleHPThermalScattering::readAnE_P_E_iso
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 << "G4ParticleHPThermalScatterin
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]->ener
232 G4double E_H = aData->vE_isoAngle[i + 1]->
233 G4double dE = E_H - E_L;
234 G4double pdf = (aData->prob[i] + aData->pr
235 total += (pdf);
236 aData->secondary_energy_cdf.push_back(tota
237 aData->secondary_energy_pdf.push_back(pdf)
238 aData->secondary_energy_value.push_back(E_
239 }
240
241 aData->sum_of_probXdEs = total;
242
243 // Normalize CDF
244 aData->secondary_energy_cdf_size = (G4int)aD
245 for (G4int i = 0; i < aData->secondary_energ
246 aData->secondary_energy_cdf[i] /= total;
247 }
248
249 return aData;
250 }
251
252 std::map<G4double, std::vector<E_isoAng*>*>*
253 G4ParticleHPThermalScattering::readAnIncoheren
254 {
255 auto T_E = new std::map<G4double, std::vecto
256
257 // std::ifstream theChannel( name.c_str() );
258 std::istringstream theChannel(std::ios::in);
259 G4ParticleHPManager::GetInstance()->GetDataS
260
261 G4int dummy;
262 while (theChannel >> dummy) // MF // Loop c
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(&the
272 T_E->insert(std::pair<G4double, std::vecto
273 }
274 // theChannel.close();
275
276 return T_E;
277 }
278
279 E_isoAng* G4ParticleHPThermalScattering::readA
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
305
306 {
307 // Select Element > Reaction >
308
309 const G4Material* theMaterial = aTrack.GetMa
310 G4double aTemp = theMaterial->GetTemperature
311 auto n = (G4int)theMaterial->GetNumberOfElem
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)(theEl
320 // Check Applicability of Thermal Scatte
321 if (getTS_ID(nullptr, theElement) != -1)
322 ielement = getTS_ID(nullptr, theElemen
323 findThermalElement = true;
324 break;
325 }
326 if (getTS_ID(theMaterial, theElement) !=
327 ielement = getTS_ID(theMaterial, theEl
328 findThermalElement = true;
329 break;
330 }
331 }
332 }
333
334 if (findThermalElement) {
335 // Select Reaction (Inelastic, coherent,
336 const G4ParticleDefinition* pd = aTrack.Ge
337 auto dp = new G4DynamicParticle(pd, aTrack
338 G4double total = theXSection->GetCrossSect
339 G4double inelastic = theXSection->GetInela
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(ieleme
348 it != inelasticFSs->find(ielement)-
349 {
350 v_temp.push_back(it->first);
351 }
352
353 std::pair<G4double, G4double> tempLH = f
354 //
355 // For T_L aNEP_EPM_TL and T_H aNEP_EPM
356 //
357 std::vector<E_P_E_isoAng*>* vNEP_EPM_TL
358 std::vector<E_P_E_isoAng*>* vNEP_EPM_TH
359
360 if (tempLH.first != 0.0 && tempLH.second
361 vNEP_EPM_TL = inelasticFSs->find(ielem
362 vNEP_EPM_TH = inelasticFSs->find(ielem
363 }
364 else if (tempLH.first == 0.0) {
365 auto itm = inelasticFSs->find(ielement
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
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 temper
385 // (continuous temperature interpolation
386 std::pair<G4double, G4double> secondaryP
387 G4double rand_temp = G4UniformRand();
388 if (rand_temp < (aTemp - tempLH.first) /
389 secondaryParam = sample_inelastic_E_mu
390 else
391 secondaryParam = sample_inelastic_E_mu
392
393 sE = secondaryParam.first;
394 mu = secondaryParam.second;
395
396 // set
397 theParticleChange.SetEnergyChange(sE);
398 G4double phi = CLHEP::twopi * G4UniformR
399 G4double sint = std::sqrt(1 - mu * mu);
400 theParticleChange.SetMomentumChange(sint
401 }
402 else if (random
403 <= (inelastic + theXSection->GetC
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(ielemen
414 it != coherentFSs->find(ielement)->
415 {
416 v_temp.push_back(it->first);
417 }
418
419 // T_L T_H
420 std::pair<G4double, G4double> tempLH = f
421 //
422 //
423 // For T_L anEPM_TL and T_H anEPM_TH
424 //
425 std::vector<std::pair<G4double, G4double
426 std::vector<std::pair<G4double, G4double
427
428 if (tempLH.first != 0.0 && tempLH.second
429 pvE_p_TL = coherentFSs->find(ielement)
430 pvE_p_TH = coherentFSs->find(ielement)
431 }
432 else if (tempLH.first == 0.0) {
433 pvE_p_TL = coherentFSs->find(ielement)
434 pvE_p_TH = coherentFSs->find(ielement)
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)
440 auto itv = v_temp.cend();
441 --itv;
442 --itv;
443 pvE_p_TL = coherentFSs->find(ielement)
444 tempLH.second = tempLH.first;
445 tempLH.first = *itv;
446 }
447 else {
448 // tempLH.first == 0.0 && tempLH.secon
449 throw G4HadronicException(
450 __FILE__, __LINE__,
451 "A problem is found in Thermal Scatt
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 interp
460 std::vector<std::pair<G4double, G4double
461 G4double rand_temp = G4UniformRand();
462 if (rand_temp < (aTemp - 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
468 for (G4int i = 0; i < n1; ++i) {
469 vE_T.push_back((*pvE_p_T_sampled)[i]->
470 vp_T.push_back((*pvE_p_T_sampled)[i]->
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 * G4UniformR
500 G4double sint = std::sqrt(1 - mu * mu);
501 theParticleChange.SetMomentumChange(sint
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(ielem
510 it != incoherentFSs->find(ielement)
511 {
512 v_temp.push_back(it->first);
513 }
514
515 // T_L T_H
516 std::pair<G4double, G4double> tempLH = f
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
526 // Interpolate TL and TH
527 anEPM_TL_E = create_E_isoAng_from_ener
528 aTrack.GetKineticEnergy(),
529 incoherentFSs->find(ielement)->secon
530 anEPM_TH_E = create_E_isoAng_from_ener
531 aTrack.GetKineticEnergy(),
532 incoherentFSs->find(ielement)->secon
533 }
534 else if (tempLH.first == 0.0) {
535 // Extrapolate T0 and T1
536 anEPM_TL_E = create_E_isoAng_from_ener
537 aTrack.GetKineticEnergy(),
538 incoherentFSs->find(ielement)->secon
539 anEPM_TH_E = create_E_isoAng_from_ener
540 aTrack.GetKineticEnergy(),
541 incoherentFSs->find(ielement)->secon
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_ener
548 aTrack.GetKineticEnergy(),
549 incoherentFSs->find(ielement)->secon
550 auto itv = v_temp.cend();
551 --itv;
552 --itv;
553 anEPM_TL_E = create_E_isoAng_from_ener
554 aTrack.GetKineticEnergy(), incoheren
555 tempLH.second = tempLH.first;
556 tempLH.first = *itv;
557 }
558
559 // E_isoAng for aTemp and aTrack.GetKine
560 G4double mu = 1.0;
561
562 // New Geant4 method - Stochastic interp
563 E_isoAng anEPM_T_E_sampled;
564 G4double rand_temp = G4UniformRand();
565 if (rand_temp < (aTemp - tempLH.first) /
566 anEPM_T_E_sampled = std::move(anEPM_TH
567 else
568 anEPM_T_E_sampled = std::move(anEPM_TL
569
570 mu = getMu(&anEPM_T_E_sampled);
571
572 // Set Final State
573 theParticleChange.SetEnergyChange(aTrack
574 G4double phi = CLHEP::twopi * G4UniformR
575 G4double sint = std::sqrt(1 - mu * mu);
576 theParticleChange.SetMomentumChange(sint
577 }
578 delete dp;
579
580 return &theParticleChange;
581 }
582 // Not thermal element
583 // Neutron HP will handle
584 return theHPElastic->ApplyYourself(aTrack, a
585 true); /
586 }
587
588 //********************************************
589 // Geant4 new algorithm
590 //********************************************
591
592 //--------------------------------------------
593 // New method added by L. Thulliez 2021 (CEA-S
594 //--------------------------------------------
595 std::pair<G4double, G4int>
596 G4ParticleHPThermalScattering::sample_inelasti
597
598 {
599 G4int i = 0;
600 G4double sE_value = 0;
601
602 for (; i < anE_P_E_isoAng->secondary_energy_
603 if (rndm1 >= anE_P_E_isoAng->secondary_ene
604 && rndm1 < anE_P_E_isoAng->secondary_e
605 {
606 G4double sE_value_i = anE_P_E_isoAng->se
607 G4double sE_pdf_i = anE_P_E_isoAng->seco
608 G4double sE_value_i1 = anE_P_E_isoAng->s
609 G4double sE_pdf_i1 = anE_P_E_isoAng->sec
610
611 G4double lambda = 0;
612 G4double alpha = (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
620 rndm = rndm2;
621 G4double gamma = -rndm;
622 G4double delta = beta * beta - 4 * alp
623
624 if (delta < 0 && std::fabs(delta) < 1.
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_val
636
637 break;
638 }
639 }
640
641 return std::pair<G4double, G4int>(sE_value,
642 }
643
644 //--------------------------------------------
645 // New method added by L. Thulliez 2021 (CEA-S
646 //--------------------------------------------
647 std::pair<G4double, G4double>
648 G4ParticleHPThermalScattering::sample_inelasti
649
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 != v
658 v_energy.push_back((*itv)->energy);
659 map_energy.insert(std::pair<G4double, G4in
660 i++;
661 }
662
663 std::pair<G4double, G4double> energyLH = fin
664
665 std::vector<E_P_E_isoAng*> pE_P_E_isoAng_lim
666
667 if (energyLH.first != 0.0 && energyLH.second
668 pE_P_E_isoAng_limit[0] = (*vNEP_EPM)[map_e
669 pE_P_E_isoAng_limit[1] = (*vNEP_EPM)[map_e
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 incid
684 G4double factor = (energyLH.second - pE) / (
685
686 if ((energyLH.second - pE) <= 0. && std::fab
687 if ((energyLH.first - pE) >= 0. && std::fabs
688
689 G4double rndm1 = G4UniformRand();
690 G4double rndm2 = G4UniformRand();
691
692 // Sample secondary neutron energy
693 std::pair<G4double, G4int> sE_lower = sample
694 std::pair<G4double, G4int> sE_upper = sample
695 G4double sE = factor * sE_lower.first + (1 -
696 sE = sE * eV;
697
698 // Sample cosine knowing the secondary neutr
699 rndm1 = G4UniformRand();
700 rndm2 = G4UniformRand();
701 G4double mu_lower = getMu(rndm1, rndm2, pE_P
702 G4double mu_upper = getMu(rndm1, rndm2, pE_P
703 G4double mu = factor * mu_lower + (1 - facto
704
705 return std::pair<G4double, G4double>(sE, mu)
706 }
707
708 //--------------------------------------------
709 // New method added by L. Thulliez 2021 (CEA-S
710 //--------------------------------------------
711 G4double G4ParticleHPThermalScattering::getMu(
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
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[(*anEP
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(
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 * ((*anEP
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[(*anEP
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> G4ParticleHPTher
776
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
789 // 3) v(imax) < temp -> LL=v(imax) H=0.0
790 for (auto it = aVector->cbegin(); it != aV
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_li
813
814
815 {
816 G4double y = 0.0;
817 if (High.first - Low.first != 0) {
818 y = (High.second - Low.second) / (High.fir
819 }
820 else {
821 if (High.second == Low.second) {
822 y = High.second;
823 }
824 else {
825 G4cout << "G4ParticleHPThermalScattering
826 }
827 }
828
829 return y;
830 }
831
832 E_isoAng G4ParticleHPThermalScattering::create
833
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->c
840 v_e.push_back((*iv)->energy);
841
842 std::pair<G4double, G4double> energyLH = fin
843 // G4cout << " " << energy/eV << " " << ener
844
845 E_isoAng* panEPM_T_EL = nullptr;
846 E_isoAng* panEPM_T_EH = nullptr;
847
848 if (energyLH.first != 0.0 && energyLH.second
849 for (auto iv = vEPM->cbegin(); iv != vEPM-
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 !=
871 // checking isoAng has proper values or no
872 // Inelastic/FS, the first and last entri
873 if (!(check_E_isoAng(panEPM_T_EL))) panEPM
874 if (!(check_E_isoAng(panEPM_T_EH))) panEPM
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; ++
881 G4double angle;
882 angle = get_linear_interpolated(
883 energy, std::pair<G4double, G4double
884 std::pair<G4double, G4double>(energy
885 anEPM_T_E.isoAngle.push_back(angle);
886 }
887 }
888 else {
889 G4Exception("G4ParticleHPThermalScatteri
890 JustWarning,
891 "G4ParticleHPThermalScatteri
892 }
893 }
894 else {
895 G4Exception("G4ParticleHPThermalScattering
896 FatalException, "Pointer panEP
897 }
898
899 return anEPM_T_E;
900 }
901
902 G4double
903 G4ParticleHPThermalScattering::get_secondary_e
904
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_isoA
919 G4double E_H = anE_P_E_isoAng->vE_isoA
920 G4double dE = E_H - E_L;
921 total += ( ( anE_P_E_isoAng->prob[i] )
922 }
923
924 if ( std::abs( total - anE_P_E_isoAng->su
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
931 G4double E_H = anE_P_E_isoAng->vE_isoAngle
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::p
938 std::pair<G4do
939 secondary_energy = secondary_energy * 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_E
950 G4double rand_for_sE, G4double pE, std::vect
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 != v
958 v_energy.push_back((*itv)->energy);
959 map_energy.insert(std::pair<G4double, G4in
960 i++;
961 }
962
963 std::pair<G4double, G4double> energyLH = fin
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
969 pE_P_E_isoAng_EL = (*vNEP_EPM)[map_energy.
970 pE_P_E_isoAng_EH = (*vNEP_EPM)[map_energy.
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_isoAn
989 sE_H = get_secondary_energy_from_E_P_E_isoAn
990
991 sE = get_linear_interpolated(pE, std::pair<G
992 std::pair<G4dou
993
994 E_isoAng E_isoAng_L = create_E_isoAng_from_e
995 E_isoAng E_isoAng_H = create_E_isoAng_from_e
996
997 E_isoAng anE_isoAng;
998 // For defeating warning message from compil
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
1007 std::pair<G4d
1008 anE_isoAng.isoAngle.push_back(angle);
1009 }
1010 }
1011 else {
1012 throw G4HadronicException(__FILE__, __LIN
1013 }
1014
1015 return std::pair<G4double, E_isoAng>(sE, an
1016 }
1017
1018 void G4ParticleHPThermalScattering::buildPhys
1019 {
1020 // Is rebuild of physics table a necessity
1021 if (nMaterial == G4Material::GetMaterialTab
1022 && nElement == G4Element::GetElementTab
1023 {
1024 return;
1025 }
1026 nMaterial = G4Material::GetMaterialTable()-
1027 nElement = G4Element::GetElementTable()->si
1028
1029 dic.clear();
1030 std::map<G4String, G4int> co_dic;
1031
1032 // Searching Nist Materials
1033 static G4ThreadLocal G4MaterialTable* theMa
1034 if (theMaterialTable == nullptr) theMateria
1035 std::size_t numberOfMaterials = G4Material:
1036 for (std::size_t i = 0; i < numberOfMateria
1037 G4Material* material = (*theMaterialTable
1038 auto numberOfElements = (G4int)material->
1039 for (G4int j = 0; j < numberOfElements; +
1040 const G4Element* element = material->Ge
1041 if (names.IsThisThermalElement(material
1042 G4int ts_ID_of_this_geometry;
1043 G4String ts_ndl_name = names.GetTS_ND
1044 if (co_dic.find(ts_ndl_name) != co_di
1045 ts_ID_of_this_geometry = co_dic.fin
1046 }
1047 else {
1048 ts_ID_of_this_geometry = (G4int)co_
1049 co_dic.insert(std::pair<G4String, G
1050 }
1051
1052 // G4cout << "Neutron HP Thermal Scat
1053 // << material->GetName() << "
1054 // << " as internal thermal sc
1055 // G4endl;
1056
1057 dic.insert(std::pair<std::pair<G4Mate
1058 std::pair<G4Material*, const G4Elem
1059 }
1060 }
1061 }
1062
1063 // Searching TS Elements
1064 auto theElementTable = G4Element::GetElemen
1065 std::size_t numberOfElements = G4Element::G
1066 for (std::size_t i = 0; i < numberOfElement
1067 const G4Element* element = (*theElementTa
1068 if (names.IsThisThermalElement(element->G
1069 if (names.IsThisThermalElement(element-
1070 G4int ts_ID_of_this_geometry;
1071 G4String ts_ndl_name = names.GetTS_ND
1072 if (co_dic.find(ts_ndl_name) != co_di
1073 ts_ID_of_this_geometry = co_dic.fin
1074 }
1075 else {
1076 ts_ID_of_this_geometry = (G4int)co_
1077 co_dic.insert(std::pair<G4String, G
1078 }
1079 dic.insert(std::pair<std::pair<const
1080 std::pair<const G4Material*, const
1081 ts_ID_of_this_geometry));
1082 }
1083 }
1084 }
1085
1086 G4cout << G4endl;
1087 G4cout
1088 << "Neutron HP Thermal Scattering: Follow
1089 << G4endl;
1090 for (const auto& it : dic) {
1091 if (it.first.first != nullptr) {
1092 G4cout << "Material " << it.first.first
1093 << it.first.second->GetName() <<
1094 << G4endl;
1095 }
1096 else {
1097 G4cout << "Element " << it.first.second
1098 << it.second << G4endl;
1099 }
1100 }
1101 G4cout << G4endl;
1102
1103 // Read Cross Section Data files
1104
1105 G4ParticleHPManager* hpmanager = G4Particle
1106 coherentFSs = hpmanager->GetThermalScatteri
1107 incoherentFSs = hpmanager->GetThermalScatte
1108 inelasticFSs = hpmanager->GetThermalScatter
1109
1110 if (G4Threading::IsMasterThread()) {
1111 clearCurrentFSData();
1112
1113 if (coherentFSs == nullptr)
1114 coherentFSs =
1115 new std::map<G4int, std::map<G4double
1116 if (incoherentFSs == nullptr)
1117 incoherentFSs = new std::map<G4int, std
1118 if (inelasticFSs == nullptr)
1119 inelasticFSs = new std::map<G4int, std:
1120
1121 G4String dirName;
1122 if (G4FindDataDir("G4NEUTRONHPDATA") == n
1123 throw G4HadronicException(
1124 __FILE__, __LINE__,
1125 "Please setenv G4NEUTRONHPDATA to poi
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/C
1134 G4String fileName = dirName + fsName +
1135 coherentFSs->insert(
1136 std::pair<G4int, std::map<G4double, s
1137 ts_ID, readACoherentFSDATA(fileName
1138
1139 // incoherent elastic
1140 fsName = "/ThermalScattering/Incoherent
1141 fileName = dirName + fsName + tsndlName
1142 incoherentFSs->insert(std::pair<G4int,
1143 ts_ID, readAnIncoherentFSDATA(fileNam
1144
1145 // inelastic
1146 fsName = "/ThermalScattering/Inelastic/
1147 fileName = dirName + fsName + tsndlName
1148 inelasticFSs->insert(std::pair<G4int, s
1149 ts_ID, readAnInelasticFSDATA(fileName
1150 }
1151
1152 hpmanager->RegisterThermalScatteringCoher
1153 hpmanager->RegisterThermalScatteringIncoh
1154 hpmanager->RegisterThermalScatteringInela
1155 }
1156
1157 theXSection->BuildPhysicsTable(*(G4Neutron:
1158 }
1159
1160 G4int G4ParticleHPThermalScattering::getTS_ID
1161 {
1162 G4int result = -1;
1163 if (dic.find(std::pair<const G4Material*, c
1164 result = dic.find(std::pair<const G4Mater
1165 return result;
1166 }
1167
1168 const std::pair<G4double, G4double> G4Particl
1169 {
1170 // max energy non-conservation is mass of h
1171 return std::pair<G4double, G4double>(10.0 *
1172 }
1173
1174 void G4ParticleHPThermalScattering::AddUserTh
1175
1176 {
1177 names.AddThermalElement(nameG4Element, file
1178 theXSection->AddUserThermalScatteringFile(n
1179 buildPhysicsTable();
1180 }
1181
1182 G4bool G4ParticleHPThermalScattering::check_E
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::ModelDesc
1197 {
1198 outFile << "High Precision model based on t
1199 << "evaluated nuclear data librarie
1200 << "on specific materials\n";
1201 }
1202