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