Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4ParticleHPThermalScattering.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// Thermal Neutron Scattering
27// Koi, Tatsumi (SLAC/SCCS)
28//
29// Class Description:
30//
31// Final State Generators for a high precision (based on evaluated data
32// libraries) description of themal neutron scattering below 4 eV;
33// Based on Thermal neutron scattering files
34// from the evaluated nuclear data files ENDF/B-VI, Release2
35// To be used in your physics list in case you need this physics.
36// In this case you want to register an object of this class with
37// the corresponding process.
38
39
40// 070625 Fix memory leaking at destructor by T. Koi
41// 081201 Fix memory leaking at destructor by T. Koi
42// 100729 Add model name in constructor Problem #1116
43// P. Arce, June-2014 Conversion neutron_hp to particle_hp
44//
45#include "G4ParticleHPThermalScattering.hh"
46#include "G4ParticleHPThermalScatteringData.hh"
47#include "G4ParticleHPThermalScatteringNames.hh"
48#include "G4ParticleHPElastic.hh"
49#include "G4ParticleHPManager.hh"
50#include "G4SystemOfUnits.hh"
51#include "G4Neutron.hh"
52#include "G4ElementTable.hh"
53#include "G4MaterialTable.hh"
54#include "G4Threading.hh"
55
56G4ParticleHPThermalScattering::G4ParticleHPThermalScattering()
57 :G4HadronicInteraction("NeutronHPThermalScattering")
58,coherentFSs(NULL)
59,incoherentFSs(NULL)
60,inelasticFSs(NULL)
61{
62 theHPElastic = new G4ParticleHPElastic();
63
64 SetMinEnergy( 0.*eV );
65 SetMaxEnergy( 4*eV );
66 theXSection = new G4ParticleHPThermalScatteringData();
67
68 //sizeOfMaterialTable = G4Material::GetMaterialTable()->size();
69 //buildPhysicsTable();
70 nMaterial = 0;
71 nElement = 0;
72}
73
74
75
76G4ParticleHPThermalScattering::~G4ParticleHPThermalScattering()
77{
78
79/*
80 for ( std::map < G4int , std::map < G4double , std::vector < E_isoAng* >* >* >::iterator it = incoherentFSs->begin() ; it != incoherentFSs->end() ; it++ )
81 {
82 std::map < G4double , std::vector < E_isoAng* >* >::iterator itt;
83 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
84 {
85 std::vector< E_isoAng* >::iterator ittt;
86 for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
87 {
88 delete *ittt;
89 }
90 delete itt->second;
91 }
92 delete it->second;
93 }
94
95 for ( std::map < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* >::iterator it = coherentFSs->begin() ; it != coherentFSs->end() ; it++ )
96 {
97 std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >::iterator itt;
98 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
99 {
100 std::vector < std::pair< G4double , G4double >* >::iterator ittt;
101 for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
102 {
103 delete *ittt;
104 }
105 delete itt->second;
106 }
107 delete it->second;
108 }
109
110 for ( std::map < G4int , std::map < G4double , std::vector < E_P_E_isoAng* >* >* >::iterator it = inelasticFSs->begin() ; it != inelasticFSs->end() ; it++ )
111 {
112 std::map < G4double , std::vector < E_P_E_isoAng* >* >::iterator itt;
113 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
114 {
115 std::vector < E_P_E_isoAng* >::iterator ittt;
116 for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
117 {
118 std::vector < E_isoAng* >::iterator it4;
119 for ( it4 = (*ittt)->vE_isoAngle.begin() ; it4 != (*ittt)->vE_isoAngle.end() ; it4++ )
120 {
121 delete *it4;
122 }
123 delete *ittt;
124 }
125 delete itt->second;
126 }
127 delete it->second;
128 }
129*/
130
131 delete theHPElastic;
132 //TKDB 160506
133 //delete theXSection;
134}
135
136void G4ParticleHPThermalScattering::clearCurrentFSData() {
137
138if ( incoherentFSs != NULL ) {
139 for ( std::map < G4int , std::map < G4double , std::vector < E_isoAng* >* >* >::iterator it = incoherentFSs->begin() ; it != incoherentFSs->end() ; it++ )
140 {
141 std::map < G4double , std::vector < E_isoAng* >* >::iterator itt;
142 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
143 {
144 std::vector< E_isoAng* >::iterator ittt;
145 for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
146 {
147 delete *ittt;
148 }
149 delete itt->second;
150 }
151 delete it->second;
152 }
153}
154
155if ( coherentFSs != NULL ) {
156 for ( std::map < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* >::iterator it = coherentFSs->begin() ; it != coherentFSs->end() ; it++ )
157 {
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}
171
172if ( inelasticFSs != NULL ) {
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}
193
194 incoherentFSs = NULL;
195 coherentFSs = NULL;
196 inelasticFSs = NULL;
197
198}
199
200
201
202void G4ParticleHPThermalScattering::BuildPhysicsTable(const G4ParticleDefinition& particle) {
203 buildPhysicsTable();
204 theHPElastic->BuildPhysicsTable( particle );
205}
206
207
208
209std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* G4ParticleHPThermalScattering::readACoherentFSDATA( G4String name )
210{
211
212 std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* aCoherentFSDATA = new std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >;
213
214 //std::ifstream theChannel( name.c_str() );
215 std::istringstream theChannel(std::ios::in);
216 G4ParticleHPManager::GetInstance()->GetDataStream(name,theChannel);
217
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;
251}
252
253
254
255std::map < G4double , std::vector < E_P_E_isoAng* >* >* G4ParticleHPThermalScattering::readAnInelasticFSDATA ( G4String name )
256{
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* >* >;
258
259 //std::ifstream theChannel( name.c_str() );
260 std::istringstream theChannel(std::ios::in);
261 G4ParticleHPManager::GetInstance()->GetDataStream(name,theChannel);
262
263 G4int dummy;
264 while ( theChannel >> dummy ) // MF // Loop checking, 11.05.2015, T. Koi
265 {
266 theChannel >> dummy; // MT
267 G4double temp;
268 theChannel >> temp;
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++ )
273 {
274 vE_P_E_isoAng->push_back ( readAnE_P_E_isoAng ( &theChannel ) );
275 }
276 anT_E_P_E_isoAng->insert ( std::pair < G4double , std::vector < E_P_E_isoAng* >* > ( temp , vE_P_E_isoAng ) );
277 }
278 //theChannel.close();
279
280 return anT_E_P_E_isoAng;
281}
282
283
284
285E_P_E_isoAng* G4ParticleHPThermalScattering::readAnE_P_E_isoAng( std::istream* file )
286{
287 E_P_E_isoAng* aData = new E_P_E_isoAng;
288
289 G4double dummy;
290 G4double energy;
291 G4int nep , nl;
292 *file >> dummy;
293 *file >> energy;
294 aData->energy = energy*eV;
295 *file >> dummy;
296 *file >> dummy;
297 *file >> nep;
298 *file >> nl;
299 aData->n = nep/nl;
300 for ( G4int i = 0 ; i < aData->n ; i++ )
301 {
302 G4double prob;
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++ )
313 {
314 G4double x;
315 *file >> x;
316 anE_isoAng->isoAngle[j] = x ;
317 //G4cout << "G4ParticleHPThermalScattering inelastic " << x << anE_isoAng->isoAngle[j] << G4endl;
318 }
319 }
320
321 // Calcuate sum_of_provXdEs
322 G4double total = 0;
323 for ( G4int i = 0 ; i < aData->n - 1 ; i++ )
324 {
325 G4double E_L = aData->vE_isoAngle[i]->energy/eV;
326 G4double E_H = aData->vE_isoAngle[i+1]->energy/eV;
327 G4double dE = E_H - E_L;
328 total += ( ( aData->prob[i] ) * dE );
329 }
330 aData->sum_of_probXdEs = total;
331
332 return aData;
333}
334
335
336
337std::map < G4double , std::vector < E_isoAng* >* >* G4ParticleHPThermalScattering::readAnIncoherentFSDATA ( G4String name )
338{
339 std::map < G4double , std::vector < E_isoAng* >* >* T_E = new std::map < G4double , std::vector < E_isoAng* >* >;
340
341 //std::ifstream theChannel( name.c_str() );
342 std::istringstream theChannel(std::ios::in);
343 G4ParticleHPManager::GetInstance()->GetDataStream(name,theChannel);
344
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();
359
360 return T_E;
361}
362
363
364
365E_isoAng* G4ParticleHPThermalScattering::readAnE_isoAng( std::istream* file )
366{
367 E_isoAng* aData = new E_isoAng;
368
369 G4double dummy;
370 G4double energy;
371 G4int n;
372 *file >> dummy;
373 *file >> energy;
374 *file >> dummy;
375 *file >> dummy;
376 *file >> n;
377 *file >> dummy;
378 aData->energy = energy*eV;
379 aData->n = n-2;
380 aData->isoAngle.resize( n );
381
382 *file >> dummy;
383 *file >> dummy;
384 for ( G4int i = 0 ; i < aData->n ; i++ )
385 *file >> aData->isoAngle[i];
386
387 return aData;
388}
389
390
391
392G4HadFinalState* G4ParticleHPThermalScattering::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aNucleus )
393{
394
395/*
396 //Trick for dynamically generated materials
397 if ( sizeOfMaterialTable != G4Material::GetMaterialTable()->size() ) {
398 sizeOfMaterialTable = G4Material::GetMaterialTable()->size();
399 buildPhysicsTable();
400 theXSection->BuildPhysicsTable( *aTrack.GetDefinition() );
401 }
402*/
403// Select Element > Reaction >
404
405 const G4Material * theMaterial = aTrack.GetMaterial();
406 G4double aTemp = theMaterial->GetTemperature();
407 G4int n = theMaterial->GetNumberOfElements();
408 //static const G4ElementTable* theElementTable = G4Element::GetElementTable();
409
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
445
446 G4double random = G4UniformRand();
447 if ( random <= inelastic/total )
448 {
449 // Inelastic
450
451 // T_L and T_H
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 }
459
460// T_L T_H
461 std::pair < G4double , G4double > tempLH = find_LH ( aTemp , &v_temp );
462//
463// For T_L aNEP_EPM_TL and T_H aNEP_EPM_TH
464//
465 std::vector< E_P_E_isoAng* >* vNEP_EPM_TL = 0;
466 std::vector< E_P_E_isoAng* >* vNEP_EPM_TH = 0;
467
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
531
532 G4double E = aTrack.GetKineticEnergy();
533
534 // T_L and T_H
535 std::map < G4double , std::vector< std::pair< G4double , G4double >* >* >::iterator it;
536 std::vector<G4double> v_temp;
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 }
542
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;
551
552 if ( tempLH.first != 0.0 && tempLH.second != 0.0 )
553 {
554 pvE_p_TL = coherentFSs->find( ielement )->second->find ( tempLH.first/kelvin )->second;
555 pvE_p_TH = coherentFSs->find( ielement )->second->find ( tempLH.first/kelvin )->second;
556 }
557 else if ( tempLH.first == 0.0 )
558 {
559 pvE_p_TL = coherentFSs->find( ielement )->second->find ( v_temp[ 0 ] )->second;
560 pvE_p_TH = coherentFSs->find( ielement )->second->find ( v_temp[ 1 ] )->second;
561 tempLH.first = tempLH.second;
562 tempLH.second = v_temp[ 1 ];
563 }
564 else if ( tempLH.second == 0.0 )
565 {
566 pvE_p_TH = coherentFSs->find( ielement )->second->find ( v_temp.back() )->second;
567 std::vector< G4double >::iterator itv;
568 itv = v_temp.end();
569 itv--;
570 itv--;
571 pvE_p_TL = coherentFSs->find( ielement )->second->find ( *itv )->second;
572 tempLH.second = tempLH.first;
573 tempLH.first = *itv;
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;
624
625 theParticleChange.SetEnergyChange( E );
626 theParticleChange.SetMomentumChange( 0.0 , std::sqrt ( 1 - mu*mu ) , mu );
627
628 }
629 else
630 {
631 // InCoherent Elastic
632
633 // T_L and T_H
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 );
644
645//
646// For T_L anEPM_TL and T_H anEPM_TH
647//
648
649 E_isoAng anEPM_TL_E;
650 E_isoAng anEPM_TH_E;
651
652 if ( tempLH.first != 0.0 && tempLH.second != 0.0 ) {
653 //Interpolate TL and TH
654 anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs->find( ielement )->second->find ( tempLH.first/kelvin )->second );
655 anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs->find( ielement )->second->find ( tempLH.second/kelvin )->second );
656 } else if ( tempLH.first == 0.0 ) {
657 //Extrapolate T0 and T1
658 anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs->find( ielement )->second->find ( v_temp[ 0 ] )->second );
659 anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs->find( ielement )->second->find ( v_temp[ 1 ] )->second );
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 }
710
711}
712
713
714
715G4double G4ParticleHPThermalScattering::getMu( E_isoAng* anEPM )
716{
717
718 G4double random = G4UniformRand();
719 G4double result = 0.0;
720
721 G4int in = int ( random * ( (*anEPM).n ) );
722
723 if ( in != 0 )
724 {
725 G4double mu_l = (*anEPM).isoAngle[ in-1 ];
726 G4double mu_h = (*anEPM).isoAngle[ in ];
727 result = ( mu_h - mu_l ) * ( random * ( (*anEPM).n ) - in ) + mu_l;
728 }
729 else
730 {
731 G4double x = random * (*anEPM).n;
732 //Bugzilla 1971
733 G4double ratio = 0.5;
734 G4double xx = G4UniformRand();
735 if ( x <= ratio )
736 {
737 G4double mu_l = -1;
738 G4double mu_h = (*anEPM).isoAngle[ 0 ];
739 result = ( mu_h - mu_l ) * xx + mu_l;
740 }
741 else
742 {
743 G4double mu_l = (*anEPM).isoAngle[ (*anEPM).n - 1 ];
744 G4double mu_h = 1;
745 result = ( mu_h - mu_l ) * xx + mu_l;
746 }
747 }
748 return result;
749}
750
751
752
753std::pair < G4double , G4double > G4ParticleHPThermalScattering::find_LH ( G4double x , std::vector< G4double >* aVector )
754{
755 G4double LL = 0.0;
756 G4double H = 0.0;
757
758 // v->size() == 1 --> LL=H=v(0)
759 if ( aVector->size() == 1 ) {
760 LL = aVector->front();
761 H = aVector->front();
762 } else {
763 // 1) temp < v(0) -> LL=0.0 H=v(0)
764 // 2) v(i-1) < temp <= v(i) -> LL=v(i-1) H=v(i)
765 // 3) v(imax) < temp -> LL=v(imax) H=0.0
766 for ( std::vector< G4double >::iterator
767 it = aVector->begin() ; it != aVector->end() ; it++ ) {
768 if ( x <= *it ) {
769 H = *it;
770 if ( it != aVector->begin() ) {
771 // 2)
772 it--;
773 LL = *it;
774 } else {
775 // 1)
776 LL = 0.0;
777 }
778 break;
779 }
780 }
781 // 3)
782 if ( H == 0.0 ) LL = aVector->back();
783 }
784
785 return std::pair < G4double , G4double > ( LL , H );
786}
787
788
789
790G4double G4ParticleHPThermalScattering::get_linear_interpolated ( G4double x , std::pair< G4double , G4double > Low , std::pair< G4double , G4double > High )
791{
792 G4double y=0.0;
793 if ( High.first - Low.first != 0 ) {
794 y = ( High.second - Low.second ) / ( High.first - Low.first ) * ( x - Low.first ) + Low.second;
795 } else {
796 if ( High.second == Low.second ) {
797 y = High.second;
798 } else {
799 G4cout << "G4ParticleHPThermalScattering liner interpolation err!!" << G4endl;
800 }
801 }
802
803 return y;
804}
805
806
807
808E_isoAng
809G4ParticleHPThermalScattering::create_E_isoAng_from_energy(G4double energy,
810 std::vector<E_isoAng*>* vEPM)
811{
812 E_isoAng anEPM_T_E;
813 std::vector<E_isoAng*>::iterator iv;
814
815 std::vector<G4double> v_e;
816 v_e.clear();
817 for (iv = vEPM->begin(); iv != vEPM->end(); iv++)
818 v_e.push_back( (*iv)->energy );
819
820 std::pair<G4double, G4double> energyLH = find_LH(energy, &v_e);
821 //G4cout << " " << energy/eV << " " << energyLH.first/eV << " " << energyLH.second/eV << G4endl;
822
823 E_isoAng* panEPM_T_EL = 0;
824 E_isoAng* panEPM_T_EH = 0;
825
826 if (energyLH.first != 0.0 && energyLH.second != 0.0) {
827 for (iv = vEPM->begin(); iv != vEPM->end(); iv++) {
828 if (energyLH.first == (*iv)->energy) {
829 panEPM_T_EL = *iv;
830 iv++;
831 panEPM_T_EH = *iv;
832 break;
833 }
834 }
835
836 } else if (energyLH.first == 0.0) {
837 panEPM_T_EL = (*vEPM)[0];
838 panEPM_T_EH = (*vEPM)[1];
839
840 } else if (energyLH.second == 0.0) {
841 panEPM_T_EH = (*vEPM).back();
842 iv = vEPM->end();
843 iv--;
844 iv--;
845 panEPM_T_EL = *iv;
846 }
847
848 if (panEPM_T_EL != 0 && panEPM_T_EH != 0) {
849 //checking isoAng has proper values or not
850 // Inelastic/FS, the first and last entries of *vEPM has all zero values.
851 if ( !(check_E_isoAng(panEPM_T_EL) ) ) panEPM_T_EL = panEPM_T_EH;
852 if ( !(check_E_isoAng(panEPM_T_EH) ) ) panEPM_T_EH = panEPM_T_EL;
853
854 if (panEPM_T_EL->n == panEPM_T_EH->n) {
855 anEPM_T_E.energy = energy;
856 anEPM_T_E.n = panEPM_T_EL->n;
857
858 for (G4int i=0; i < panEPM_T_EL->n; i++) {
859 G4double angle;
860 angle = get_linear_interpolated(energy, std::pair<G4double,G4double>(energyLH.first, panEPM_T_EL->isoAngle[i] ),
861 std::pair<G4double,G4double>(energyLH.second, panEPM_T_EH->isoAngle[i] ) );
862 anEPM_T_E.isoAngle.push_back(angle);
863 }
864
865 } else {
866 G4Exception("G4ParticleHPThermalScattering::create_E_isoAng_from_energy",
867 "NotSupported", JustWarning,
868 "G4ParticleHPThermalScattering does not support yet EL->n != EH->n.");
869 }
870
871 } else {
872 G4Exception("G4ParticleHPThermalScattering::create_E_isoAng_from_energy",
873 "HAD_THERM_000", FatalException,
874 "Pointer panEPM_T_EL or panEPM_T_EH is zero");
875 }
876
877 return anEPM_T_E;
878}
879
880
881
882G4double 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 );
913
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
928std::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);
1006}
1007
1008void G4ParticleHPThermalScattering::buildPhysicsTable()
1009{
1010
1011 //Is rebuild of physics table a necessity
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;
1094 }
1095 else
1096 {
1097 G4cout << "Element " << it->first.second->GetName() << ", internal thermal scattering id " << it->second << G4endl;
1098 }
1099 }
1100 G4cout << G4endl;
1101
1102 // Read Cross Section Data files
1103
1104 G4ParticleHPManager* hpmanager = G4ParticleHPManager::GetInstance();
1105 coherentFSs = hpmanager->GetThermalScatteringCoherentFinalStates();
1106 incoherentFSs = hpmanager->GetThermalScatteringIncoherentFinalStates();
1107 inelasticFSs = hpmanager->GetThermalScatteringInelasticFinalStates();
1108
1109 if ( G4Threading::IsMasterThread() ) {
1110
1111 clearCurrentFSData();
1112
1113 if ( coherentFSs == NULL ) coherentFSs = new std::map < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* >;
1114 if ( incoherentFSs == NULL ) incoherentFSs = new std::map < G4int , std::map < G4double , std::vector < E_isoAng* >* >* >;
1115 if ( inelasticFSs == NULL ) inelasticFSs = new std::map < G4int , std::map < G4double , std::vector < E_P_E_isoAng* >* >* >;
1116
1117 G4String dirName;
1118 if ( !std::getenv( "G4NEUTRONHPDATA" ) )
1119 throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
1120 dirName = std::getenv( "G4NEUTRONHPDATA" );
1121
1122 //G4String name;
1123
1124 for ( std::map < G4String , G4int >::iterator it = co_dic.begin() ; it != co_dic.end() ; it++ )
1125 {
1126 G4String tsndlName = it->first;
1127 G4int ts_ID = it->second;
1128
1129 // Coherent
1130 G4String fsName = "/ThermalScattering/Coherent/FS/";
1131 G4String fileName = dirName + fsName + tsndlName;
1132 coherentFSs->insert ( std::pair < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* > ( ts_ID , readACoherentFSDATA( fileName ) ) );
1133
1134 // incoherent elastic
1135 fsName = "/ThermalScattering/Incoherent/FS/";
1136 fileName = dirName + fsName + tsndlName;
1137 incoherentFSs->insert ( std::pair < G4int , std::map < G4double , std::vector < E_isoAng* >* >* > ( ts_ID , readAnIncoherentFSDATA( fileName ) ) );
1138
1139 // inelastic
1140 fsName = "/ThermalScattering/Inelastic/FS/";
1141 fileName = dirName + fsName + tsndlName;
1142 inelasticFSs->insert ( std::pair < G4int , std::map < G4double , std::vector < E_P_E_isoAng* >* >* > ( ts_ID , readAnInelasticFSDATA( fileName ) ) );
1143 }
1144
1145 hpmanager->RegisterThermalScatteringCoherentFinalStates( coherentFSs );
1146 hpmanager->RegisterThermalScatteringIncoherentFinalStates( incoherentFSs );
1147 hpmanager->RegisterThermalScatteringInelasticFinalStates( inelasticFSs );
1148 }
1149
1150 theXSection->BuildPhysicsTable( *(G4Neutron::Neutron()) );
1151}
1152
1153
1154G4int G4ParticleHPThermalScattering::getTS_ID ( const G4Material* material , const G4Element* element )
1155{
1156 G4int result = -1;
1157 if ( dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) ) != dic.end() )
1158 result = dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) )->second;
1159 return result;
1160}
1161
1162const std::pair<G4double, G4double> G4ParticleHPThermalScattering::GetFatalEnergyCheckLevels() const
1163{
1164 //return std::pair<G4double, G4double>(10*perCent,10*GeV);
1165 return std::pair<G4double, G4double>(10*perCent,DBL_MAX);
1166}
1167
1168void G4ParticleHPThermalScattering::AddUserThermalScatteringFile( G4String nameG4Element , G4String filename)
1169{
1170 names.AddThermalElement( nameG4Element , filename );
1171 theXSection->AddUserThermalScatteringFile( nameG4Element , filename );
1172 buildPhysicsTable();
1173}
1174
1175
1176G4bool G4ParticleHPThermalScattering::check_E_isoAng( E_isoAng* anE_IsoAng )
1177{
1178 G4bool result=false;
1179
1180 G4int n = anE_IsoAng->n;
1181 G4double sum=0.0;
1182 for ( G4int i = 0 ; i < n ; i++ ) {
1183 sum += anE_IsoAng->isoAngle[ i ];
1184 }
1185 if ( sum != 0.0 ) result = true;
1186
1187 return result;
1188}
1189
1190void G4ParticleHPThermalScattering::ModelDescription(std::ostream& outFile) const
1191{
1192 outFile << "High Precision model based on thermal scattering data in\n"
1193 << "evaluated nuclear data libraries for neutrons below 5eV\n"
1194 << "on specific materials\n";
1195}
G4ElementTable
std::vector< G4Element * > G4ElementTable
Definition: G4ElementTable.hh:44
JustWarning
@ JustWarning
Definition: G4ExceptionSeverity.hh:65
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:61
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
G4MaterialTable
std::vector< G4Material * > G4MaterialTable
Definition: G4MaterialTable.hh:44
G4double
double G4double
Definition: G4Types.hh:83
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
G4Element::GetElementTable
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:397
G4Element::GetZ
G4double GetZ() const
Definition: G4Element.hh:130
G4Element::GetNumberOfElements
static size_t GetNumberOfElements()
Definition: G4Element.cc:404
G4Element::GetName
const G4String & GetName() const
Definition: G4Element.hh:126
G4HadFinalState::SetEnergyChange
void SetEnergyChange(G4double anEnergy)
Definition: G4HadFinalState.cc:39
G4HadFinalState::SetMomentumChange
void SetMomentumChange(const G4ThreeVector &aV)
Definition: G4HadFinalState.hh:56
G4HadProjectile::GetMaterial
const G4Material * GetMaterial() const
Definition: G4HadProjectile.hh:81
G4HadProjectile::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4HadProjectile.hh:86
G4HadProjectile::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4HadProjectile.hh:116
G4HadProjectile::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4HadProjectile.hh:91
G4HadronicInteraction::SetMinEnergy
void SetMinEnergy(G4double anEnergy)
Definition: G4HadronicInteraction.hh:89
G4HadronicInteraction::SetMaxEnergy
void SetMaxEnergy(const G4double anEnergy)
Definition: G4HadronicInteraction.hh:102
G4Material::GetNumberOfMaterials
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:644
G4Material::GetTemperature
G4double GetTemperature() const
Definition: G4Material.hh:180
G4Material::GetElement
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:200
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
G4Material::GetMaterialTable
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:637
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:175
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:103
G4Nucleus::GetZ_asInt
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4ParticleHPElastic::BuildPhysicsTable
void BuildPhysicsTable(const G4ParticleDefinition &)
Definition: G4ParticleHPElastic.cc:183
G4ParticleHPManager::RegisterThermalScatteringCoherentFinalStates
void RegisterThermalScatteringCoherentFinalStates(std::map< G4int, std::map< G4double, std::vector< std::pair< G4double, G4double > * > * > * > *val)
Definition: G4ParticleHPManager.hh:125
G4ParticleHPManager::RegisterThermalScatteringIncoherentFinalStates
void RegisterThermalScatteringIncoherentFinalStates(std::map< G4int, std::map< G4double, std::vector< E_isoAng * > * > * > *val)
Definition: G4ParticleHPManager.hh:127
G4ParticleHPManager::RegisterThermalScatteringInelasticFinalStates
void RegisterThermalScatteringInelasticFinalStates(std::map< G4int, std::map< G4double, std::vector< E_P_E_isoAng * > * > * > *val)
Definition: G4ParticleHPManager.hh:129
G4ParticleHPManager::GetThermalScatteringCoherentFinalStates
std::map< G4int, std::map< G4double, std::vector< std::pair< G4double, G4double > * > * > * > * GetThermalScatteringCoherentFinalStates()
Definition: G4ParticleHPManager.hh:124
G4ParticleHPManager::GetInstance
static G4ParticleHPManager * GetInstance()
Definition: G4ParticleHPManager.cc:140
G4ParticleHPManager::GetThermalScatteringIncoherentFinalStates
std::map< G4int, std::map< G4double, std::vector< E_isoAng * > * > * > * GetThermalScatteringIncoherentFinalStates()
Definition: G4ParticleHPManager.hh:126
G4ParticleHPManager::GetDataStream
void GetDataStream(G4String, std::istringstream &iss)
Definition: G4ParticleHPManager.cc:165
G4ParticleHPManager::GetThermalScatteringInelasticFinalStates
std::map< G4int, std::map< G4double, std::vector< E_P_E_isoAng * > * > * > * GetThermalScatteringInelasticFinalStates()
Definition: G4ParticleHPManager.hh:128
G4ParticleHPThermalScatteringData::GetCrossSection
G4double GetCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
Definition: G4ParticleHPThermalScatteringData.cc:450
G4ParticleHPThermalScatteringData::BuildPhysicsTable
void BuildPhysicsTable(const G4ParticleDefinition &)
Definition: G4ParticleHPThermalScatteringData.cc:215
G4ParticleHPThermalScatteringData::GetCoherentCrossSection
G4double GetCoherentCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
Definition: G4ParticleHPThermalScatteringData.cc:481
G4ParticleHPThermalScatteringData::GetInelasticCrossSection
G4double GetInelasticCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
Definition: G4ParticleHPThermalScatteringData.cc:472
G4ParticleHPThermalScattering::ApplyYourself
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus)
Definition: G4ParticleHPThermalScattering.cc:392
G4ParticleHPThermalScattering::AddUserThermalScatteringFile
void AddUserThermalScatteringFile(G4String, G4String)
Definition: G4ParticleHPThermalScattering.cc:1168
G4ParticleHPThermalScattering::ModelDescription
virtual void ModelDescription(std::ostream &outFile) const
Definition: G4ParticleHPThermalScattering.cc:1190
G4ParticleHPThermalScattering::GetFatalEnergyCheckLevels
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels() const
Definition: G4ParticleHPThermalScattering.cc:1162
G4ParticleHPThermalScattering::BuildPhysicsTable
void BuildPhysicsTable(const G4ParticleDefinition &)
Definition: G4ParticleHPThermalScattering.cc:202
G4String::first
std::size_t first(char) const
G4INCL::CrossSections::total
G4double total(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:60
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:158
G4Threading::IsMasterThread
G4bool IsMasterThread()
Definition: G4Threading.cc:124
E_P_E_isoAng
Definition: G4ParticleHPThermalScattering.hh:64
E_P_E_isoAng::prob
std::vector< G4double > prob
Definition: G4ParticleHPThermalScattering.hh:67
E_P_E_isoAng::vE_isoAngle
std::vector< E_isoAng * > vE_isoAngle
Definition: G4ParticleHPThermalScattering.hh:68
E_P_E_isoAng::n
G4int n
Definition: G4ParticleHPThermalScattering.hh:66
E_P_E_isoAng::sum_of_probXdEs
G4double sum_of_probXdEs
Definition: G4ParticleHPThermalScattering.hh:69
E_P_E_isoAng::energy
G4double energy
Definition: G4ParticleHPThermalScattering.hh:65
E_isoAng
Definition: G4ParticleHPThermalScattering.hh:53
E_isoAng::isoAngle
std::vector< G4double > isoAngle
Definition: G4ParticleHPThermalScattering.hh:56
E_isoAng::energy
G4double energy
Definition: G4ParticleHPThermalScattering.hh:54
E_isoAng::n
G4int n
Definition: G4ParticleHPThermalScattering.hh:55
DBL_MAX
#define DBL_MAX
Definition: templates.hh:62
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:77