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