Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 
56 G4ParticleHPThermalScattering::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 
76 G4ParticleHPThermalScattering::~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 
136 void G4ParticleHPThermalScattering::clearCurrentFSData() {
137 
138 if ( 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 
155 if ( 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 
172 if ( 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 
202 void G4ParticleHPThermalScattering::BuildPhysicsTable(const G4ParticleDefinition& particle) {
203  buildPhysicsTable();
204  theHPElastic->BuildPhysicsTable( particle );
205 }
206 
207 
208 
209 std::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 
255 std::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 
285 E_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 
337 std::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 
365 E_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 
392 G4HadFinalState* 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  for ( G4int i=1 ; i < n1 ; i++ )
588  {
589  if ( (*pvE_p_TL)[i]->first != (*pvE_p_TH)[i]->first ) throw G4HadronicException(__FILE__, __LINE__, "A problem is found in Thermal Scattering Data!");
590  vE_T.push_back ( (*pvE_p_TL)[i]->first );
591  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 ) ) );
592  }
593 
594  G4int j = 0;
595  for ( G4int i = 1 ; i < n ; i++ )
596  {
597  if ( E/eV < vE_T[ i ] )
598  {
599  j = i-1;
600  break;
601  }
602  }
603 
604  G4double rand_for_mu = G4UniformRand();
605 
606  G4int k = 0;
607  for ( G4int i = 1 ; i < j ; i++ )
608  {
609  G4double Pi = vp_T[ i ] / vp_T[ j ];
610  if ( rand_for_mu < Pi )
611  {
612  k = i-1;
613  break;
614  }
615  }
616 
617  //G4double Ei = vE_T[ j ];
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 
624  theParticleChange.SetEnergyChange( E );
625  theParticleChange.SetMomentumChange( 0.0 , std::sqrt ( 1 - mu*mu ) , mu );
626 
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 
715 G4double 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  G4double D = ( (*anEPM).isoAngle[ 0 ] - ( -1 ) ) + ( 1 - (*anEPM).isoAngle[ (*anEPM).n - 1 ] );
733  G4double ratio = ( (*anEPM).isoAngle[ 0 ] - ( -1 ) ) / D;
734  if ( x <= ratio )
735  {
736  G4double mu_l = -1;
737  G4double mu_h = (*anEPM).isoAngle[ 0 ];
738  result = ( mu_h - mu_l ) * x + mu_l;
739  }
740  else
741  {
742  G4double mu_l = (*anEPM).isoAngle[ (*anEPM).n - 1 ];
743  G4double mu_h = 1;
744  result = ( mu_h - mu_l ) * x + mu_l;
745  }
746  }
747  return result;
748 }
749 
750 
751 
752 std::pair < G4double , G4double > G4ParticleHPThermalScattering::find_LH ( G4double x , std::vector< G4double >* aVector )
753 {
754  G4double LL = 0.0;
755  G4double H = 0.0;
756 
757  // v->size() == 1 --> LL=H=v(0)
758  if ( aVector->size() == 1 ) {
759  LL = aVector->front();
760  H = aVector->front();
761  } else {
762  // 1) temp < v(0) -> LL=0.0 H=v(0)
763  // 2) v(i-1) < temp <= v(i) -> LL=v(i-1) H=v(i)
764  // 3) v(imax) < temp -> LL=v(imax) H=0.0
765  for ( std::vector< G4double >::iterator
766  it = aVector->begin() ; it != aVector->end() ; it++ ) {
767  if ( x <= *it ) {
768  H = *it;
769  if ( it != aVector->begin() ) {
770  // 2)
771  it--;
772  LL = *it;
773  } else {
774  // 1)
775  LL = 0.0;
776  }
777  break;
778  }
779  }
780  // 3)
781  if ( H == 0.0 ) LL = aVector->back();
782  }
783 
784  return std::pair < G4double , G4double > ( LL , H );
785 }
786 
787 
788 
789 G4double G4ParticleHPThermalScattering::get_linear_interpolated ( G4double x , std::pair< G4double , G4double > Low , std::pair< G4double , G4double > High )
790 {
791  G4double y=0.0;
792  if ( High.first - Low.first != 0 ) {
793  y = ( High.second - Low.second ) / ( High.first - Low.first ) * ( x - Low.first ) + Low.second;
794  } else {
795  if ( High.second == Low.second ) {
796  y = High.second;
797  } else {
798  G4cout << "G4ParticleHPThermalScattering liner interpolation err!!" << G4endl;
799  }
800  }
801 
802  return y;
803 }
804 
805 
806 
807 E_isoAng G4ParticleHPThermalScattering::create_E_isoAng_from_energy ( G4double energy , std::vector< E_isoAng* >* vEPM )
808 {
809  E_isoAng anEPM_T_E;
810 
811  std::vector< E_isoAng* >::iterator iv;
812 
813  std::vector< G4double > v_e;
814  v_e.clear();
815  for ( iv = vEPM->begin() ; iv != vEPM->end() ; iv++ )
816  v_e.push_back ( (*iv)->energy );
817 
818  std::pair < G4double , G4double > energyLH = find_LH ( energy , &v_e );
819  //G4cout << " " << energy/eV << " " << energyLH.first/eV << " " << energyLH.second/eV << G4endl;
820 
821  E_isoAng* panEPM_T_EL=0;
822  E_isoAng* panEPM_T_EH=0;
823 
824  if ( energyLH.first != 0.0 && energyLH.second != 0.0 )
825  {
826  for ( iv = vEPM->begin() ; iv != vEPM->end() ; iv++ )
827  {
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  {
838  panEPM_T_EL = (*vEPM)[0];
839  panEPM_T_EH = (*vEPM)[1];
840  }
841  else if ( energyLH.second == 0.0 )
842  {
843  panEPM_T_EH = (*vEPM).back();
844  iv = vEPM->end();
845  iv--;
846  iv--;
847  panEPM_T_EL = *iv;
848  }
849 
850  //checking isoAng has proper values or not
851  // Inelastic/FS, the first and last entries of *vEPM has all zero values.
852  if ( ! ( check_E_isoAng (panEPM_T_EL) ) ) panEPM_T_EL= panEPM_T_EH;
853  if ( ! ( check_E_isoAng (panEPM_T_EH) ) ) panEPM_T_EH= panEPM_T_EL;
854 
855  if ( panEPM_T_EL->n == panEPM_T_EH->n )
856  {
857  anEPM_T_E.energy = energy;
858  anEPM_T_E.n = panEPM_T_EL->n;
859 
860  for ( G4int i=0 ; i < panEPM_T_EL->n ; i++ )
861  {
862  G4double angle;
863  angle = get_linear_interpolated ( energy , std::pair< G4double , G4double > ( energyLH.first , panEPM_T_EL->isoAngle[ i ] ) , std::pair< G4double , G4double > ( energyLH.second , panEPM_T_EH->isoAngle[ i ] ) );
864  anEPM_T_E.isoAngle.push_back( angle );
865  }
866  }
867  else
868  {
869  G4cout << "G4ParticleHPThermalScattering Do not Suuport yet." << G4endl;
870  }
871 
872 
873  return anEPM_T_E;
874 }
875 
876 
877 
878 G4double G4ParticleHPThermalScattering::get_secondary_energy_from_E_P_E_isoAng ( G4double random , E_P_E_isoAng* anE_P_E_isoAng )
879 {
880 
881  G4double secondary_energy = 0.0;
882 
883  G4int n = anE_P_E_isoAng->n;
884  G4double sum_p = 0.0; // sum_p_H
885  G4double sum_p_L = 0.0;
886 
887  G4double total=0.0;
888 
889 /*
890  delete for speed up
891  for ( G4int i = 0 ; i < n-1 ; i++ )
892  {
893  G4double E_L = anE_P_E_isoAng->vE_isoAngle[i]->energy/eV;
894  G4double E_H = anE_P_E_isoAng->vE_isoAngle[i+1]->energy/eV;
895  G4double dE = E_H - E_L;
896  total += ( ( anE_P_E_isoAng->prob[i] ) * dE );
897  }
898 
899  if ( std::abs( total - anE_P_E_isoAng->sum_of_probXdEs ) > 1.0e-14 ) G4cout << total - anE_P_E_isoAng->sum_of_probXdEs << G4endl;
900 */
901  total = anE_P_E_isoAng->sum_of_probXdEs;
902 
903  for ( G4int i = 0 ; i < n-1 ; i++ )
904  {
905  G4double E_L = anE_P_E_isoAng->vE_isoAngle[i]->energy/eV;
906  G4double E_H = anE_P_E_isoAng->vE_isoAngle[i+1]->energy/eV;
907  G4double dE = E_H - E_L;
908  sum_p += ( ( anE_P_E_isoAng->prob[i] ) * dE );
909 
910  if ( random <= sum_p/total )
911  {
912  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 ) );
913  secondary_energy = secondary_energy*eV; //need eV
914  break;
915  }
916  sum_p_L = sum_p;
917  }
918 
919  return secondary_energy;
920 }
921 
922 
923 
924 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 )
925 {
926 
927  std::map< G4double , G4int > map_energy;
928  map_energy.clear();
929  std::vector< G4double > v_energy;
930  v_energy.clear();
931  std::vector< E_P_E_isoAng* >::iterator itv;
932  G4int i = 0;
933  for ( itv = vNEP_EPM->begin(); itv != vNEP_EPM->end(); itv++ )
934  {
935  v_energy.push_back( (*itv)->energy );
936  map_energy.insert( std::pair < G4double , G4int > ( (*itv)->energy , i ) );
937  i++;
938  }
939 
940  std::pair < G4double , G4double > energyLH = find_LH ( pE , &v_energy );
941 
942  E_P_E_isoAng* pE_P_E_isoAng_EL = 0;
943  E_P_E_isoAng* pE_P_E_isoAng_EH = 0;
944 
945  if ( energyLH.first != 0.0 && energyLH.second != 0.0 )
946  {
947  pE_P_E_isoAng_EL = (*vNEP_EPM)[ map_energy.find ( energyLH.first )->second ];
948  pE_P_E_isoAng_EH = (*vNEP_EPM)[ map_energy.find ( energyLH.second )->second ];
949  }
950  else if ( energyLH.first == 0.0 )
951  {
952  pE_P_E_isoAng_EL = (*vNEP_EPM)[ 0 ];
953  pE_P_E_isoAng_EH = (*vNEP_EPM)[ 1 ];
954  }
955  if ( energyLH.second == 0.0 )
956  {
957  pE_P_E_isoAng_EH = (*vNEP_EPM).back();
958  itv = vNEP_EPM->end();
959  itv--;
960  itv--;
961  pE_P_E_isoAng_EL = *itv;
962  }
963 
964 
965  G4double sE;
966  G4double sE_L;
967  G4double sE_H;
968 
969 
970  sE_L = get_secondary_energy_from_E_P_E_isoAng ( rand_for_sE , pE_P_E_isoAng_EL );
971  sE_H = get_secondary_energy_from_E_P_E_isoAng ( rand_for_sE , pE_P_E_isoAng_EH );
972 
973  sE = get_linear_interpolated ( pE , std::pair < G4double , G4double > ( energyLH.first , sE_L ) , std::pair < G4double , G4double > ( energyLH.second , sE_H ) );
974 
975 
976  E_isoAng E_isoAng_L = create_E_isoAng_from_energy ( sE , &(pE_P_E_isoAng_EL->vE_isoAngle) );
977  E_isoAng E_isoAng_H = create_E_isoAng_from_energy ( sE , &(pE_P_E_isoAng_EH->vE_isoAngle) );
978 
979  E_isoAng anE_isoAng;
980  //For defeating warning message from compiler
981  anE_isoAng.n = 1;
982  anE_isoAng.energy = sE; //never used
983  if ( E_isoAng_L.n == E_isoAng_H.n )
984  {
985  anE_isoAng.n = E_isoAng_L.n;
986  for ( G4int j=0 ; j < anE_isoAng.n ; j++ )
987  {
988  G4double angle;
989  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 ] ) );
990  anE_isoAng.isoAngle.push_back( angle );
991  }
992  }
993  else
994  {
995  //G4cout << "Do not Suuport yet." << G4endl;
996  throw G4HadronicException(__FILE__, __LINE__, "Unexpected values!");
997  }
998 
999 
1000 
1001  return std::pair< G4double , E_isoAng >( sE , anE_isoAng);
1002 }
1003 
1004 void G4ParticleHPThermalScattering::buildPhysicsTable()
1005 {
1006 
1007  //Is rebuild of physics table a necessity
1008  if ( nMaterial == G4Material::GetMaterialTable()->size() && nElement == G4Element::GetElementTable()->size() ) {
1009  return;
1010  } else {
1011  nMaterial = G4Material::GetMaterialTable()->size();
1012  nElement = G4Element::GetElementTable()->size();
1013  }
1014 
1015  dic.clear();
1016  std::map < G4String , G4int > co_dic;
1017 
1018  //Searching Nist Materials
1019  static G4ThreadLocal G4MaterialTable* theMaterialTable = 0 ; if (!theMaterialTable) theMaterialTable= G4Material::GetMaterialTable();
1020  size_t numberOfMaterials = G4Material::GetNumberOfMaterials();
1021  for ( size_t i = 0 ; i < numberOfMaterials ; i++ )
1022  {
1023  G4Material* material = (*theMaterialTable)[i];
1024  size_t numberOfElements = material->GetNumberOfElements();
1025  for ( size_t j = 0 ; j < numberOfElements ; j++ )
1026  {
1027  const G4Element* element = material->GetElement(j);
1028  if ( names.IsThisThermalElement ( material->GetName() , element->GetName() ) )
1029  {
1030  G4int ts_ID_of_this_geometry;
1031  G4String ts_ndl_name = names.GetTS_NDL_Name( material->GetName() , element->GetName() );
1032  if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
1033  {
1034  ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) -> second;
1035  }
1036  else
1037  {
1038  ts_ID_of_this_geometry = co_dic.size();
1039  co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
1040  }
1041 
1042  //G4cout << "Neutron HP Thermal Scattering: Registering a material-element pair of "
1043  // << material->GetName() << " " << element->GetName()
1044  // << " as internal thermal scattering id of " << ts_ID_of_this_geometry << "." << G4endl;
1045 
1046  dic.insert( std::pair < std::pair < G4Material* , const G4Element* > , G4int > ( std::pair < G4Material* , const G4Element* > ( material , element ) , ts_ID_of_this_geometry ) );
1047  }
1048  }
1049  }
1050 
1051  //Searching TS Elements
1052  static G4ThreadLocal G4ElementTable* theElementTable = 0 ; if (!theElementTable) theElementTable= G4Element::GetElementTable();
1053  size_t numberOfElements = G4Element::GetNumberOfElements();
1054  //size_t numberOfThermalElements = 0;
1055  for ( size_t i = 0 ; i < numberOfElements ; i++ )
1056  {
1057  const G4Element* element = (*theElementTable)[i];
1058  if ( names.IsThisThermalElement ( element->GetName() ) )
1059  {
1060  if ( names.IsThisThermalElement ( element->GetName() ) )
1061  {
1062  G4int ts_ID_of_this_geometry;
1063  G4String ts_ndl_name = names.GetTS_NDL_Name( element->GetName() );
1064  if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
1065  {
1066  ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) -> second;
1067  }
1068  else
1069  {
1070  ts_ID_of_this_geometry = co_dic.size();
1071  co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
1072  }
1073 
1074  //G4cout << "Neutron HP Thermal Scattering: Registering an element of "
1075  // << material->GetName() << " " << element->GetName()
1076  // << " as internal thermal scattering id of " << ts_ID_of_this_geometry << "." << G4endl;
1077 
1078  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 ) );
1079  }
1080  }
1081  }
1082 
1083  G4cout << G4endl;
1084  G4cout << "Neutron HP Thermal Scattering: Following material-element pairs or elements are registered." << G4endl;
1085  for ( std::map < std::pair < const G4Material* , const G4Element* > , G4int >::iterator it = dic.begin() ; it != dic.end() ; it++ )
1086  {
1087  if ( it->first.first != NULL )
1088  {
1089  G4cout << "Material " << it->first.first->GetName() << " - Element " << it->first.second->GetName() << ", internal thermal scattering id " << it->second << G4endl;
1090  }
1091  else
1092  {
1093  G4cout << "Element " << it->first.second->GetName() << ", internal thermal scattering id " << it->second << G4endl;
1094  }
1095  }
1096  G4cout << G4endl;
1097 
1098  // Read Cross Section Data files
1099 
1100  G4ParticleHPManager* hpmanager = G4ParticleHPManager::GetInstance();
1101  coherentFSs = hpmanager->GetThermalScatteringCoherentFinalStates();
1102  incoherentFSs = hpmanager->GetThermalScatteringIncoherentFinalStates();
1103  inelasticFSs = hpmanager->GetThermalScatteringInelasticFinalStates();
1104 
1105  if ( G4Threading::IsMasterThread() ) {
1106 
1107  clearCurrentFSData();
1108 
1109  if ( coherentFSs == NULL ) coherentFSs = new std::map < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* >;
1110  if ( incoherentFSs == NULL ) incoherentFSs = new std::map < G4int , std::map < G4double , std::vector < E_isoAng* >* >* >;
1111  if ( inelasticFSs == NULL ) inelasticFSs = new std::map < G4int , std::map < G4double , std::vector < E_P_E_isoAng* >* >* >;
1112 
1113  G4String dirName;
1114  if ( !getenv( "G4NEUTRONHPDATA" ) )
1115  throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
1116  dirName = getenv( "G4NEUTRONHPDATA" );
1117 
1118  //G4String name;
1119 
1120  for ( std::map < G4String , G4int >::iterator it = co_dic.begin() ; it != co_dic.end() ; it++ )
1121  {
1122  G4String tsndlName = it->first;
1123  G4int ts_ID = it->second;
1124 
1125  // Coherent
1126  G4String fsName = "/ThermalScattering/Coherent/FS/";
1127  G4String fileName = dirName + fsName + tsndlName;
1128  coherentFSs->insert ( std::pair < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* > ( ts_ID , readACoherentFSDATA( fileName ) ) );
1129 
1130  // incoherent elastic
1131  fsName = "/ThermalScattering/Incoherent/FS/";
1132  fileName = dirName + fsName + tsndlName;
1133  incoherentFSs->insert ( std::pair < G4int , std::map < G4double , std::vector < E_isoAng* >* >* > ( ts_ID , readAnIncoherentFSDATA( fileName ) ) );
1134 
1135  // inelastic
1136  fsName = "/ThermalScattering/Inelastic/FS/";
1137  fileName = dirName + fsName + tsndlName;
1138  inelasticFSs->insert ( std::pair < G4int , std::map < G4double , std::vector < E_P_E_isoAng* >* >* > ( ts_ID , readAnInelasticFSDATA( fileName ) ) );
1139  }
1140 
1141  hpmanager->RegisterThermalScatteringCoherentFinalStates( coherentFSs );
1142  hpmanager->RegisterThermalScatteringIncoherentFinalStates( incoherentFSs );
1143  hpmanager->RegisterThermalScatteringInelasticFinalStates( inelasticFSs );
1144  }
1145 
1146  theXSection->BuildPhysicsTable( *(G4Neutron::Neutron()) );
1147 }
1148 
1149 
1150 G4int G4ParticleHPThermalScattering::getTS_ID ( const G4Material* material , const G4Element* element )
1151 {
1152  G4int result = -1;
1153  if ( dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) ) != dic.end() )
1154  result = dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) )->second;
1155  return result;
1156 }
1157 
1158 const std::pair<G4double, G4double> G4ParticleHPThermalScattering::GetFatalEnergyCheckLevels() const
1159 {
1160  //return std::pair<G4double, G4double>(10*perCent,10*GeV);
1161  return std::pair<G4double, G4double>(10*perCent,DBL_MAX);
1162 }
1163 
1164 void G4ParticleHPThermalScattering::AddUserThermalScatteringFile( G4String nameG4Element , G4String filename)
1165 {
1166  names.AddThermalElement( nameG4Element , filename );
1167  theXSection->AddUserThermalScatteringFile( nameG4Element , filename );
1168  buildPhysicsTable();
1169 }
1170 
1171 
1172 G4bool G4ParticleHPThermalScattering::check_E_isoAng( E_isoAng* anE_IsoAng )
1173 {
1174  G4bool result=false;
1175 
1176  G4int n = anE_IsoAng->n;
1177  G4double sum=0.0;
1178  for ( G4int i = 0 ; i < n ; i++ ) {
1179  sum += anE_IsoAng->isoAngle[ i ];
1180  }
1181  if ( sum != 0.0 ) result = true;
1182 
1183  return result;
1184 }
1185 
1186 void G4ParticleHPThermalScattering::ModelDescription(std::ostream& outFile) const
1187 {
1188  outFile << "High Precision model based on thermal scattering data in evaluated nuclear data libraries for neutrons below 5eV on specific materials\n";
1189 }
result
G4double G4ParticleHPJENDLHEData::G4double result
Definition: G4ParticleHPJENDLHEData.cc:257
G4ParticleHPManager::GetInstance
static G4ParticleHPManager * GetInstance()
Definition: G4ParticleHPManager.hh:56
name
const XML_Char * name
Definition: expat.h:151
E_P_E_isoAng::energy
G4double energy
Definition: G4ParticleHPThermalScattering.hh:65
E_isoAng::isoAngle
std::vector< G4double > isoAngle
Definition: G4ParticleHPThermalScattering.hh:56
G4ParticleHPThermalScatteringNames.hh
G4String::first
G4int first(char) const
G4ParticleHPElastic::BuildPhysicsTable
void BuildPhysicsTable(const G4ParticleDefinition &)
Definition: G4ParticleHPElastic.cc:180
perCent
static constexpr double perCent
Definition: G4SIunits.hh:332
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:178
G4Element::GetZ
G4double GetZ() const
Definition: G4Element.hh:131
G4Material::GetMaterialTable
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
G4MaterialTable
std::vector< G4Material * > G4MaterialTable
Definition: G4MaterialTable.hh:41
first
G4int first
Definition: G4ParticleHPJENDLHEData.cc:280
E_isoAng
Definition: G4ParticleHPThermalScattering.hh:52
G4ParticleHPManager::GetDataStream
void GetDataStream(G4String, std::istringstream &iss)
Definition: G4ParticleHPManager.cc:100
angle
static G4double angle[DIM]
Definition: G4ScreeningMottCrossSection.cc:80
E_P_E_isoAng::n
G4int n
Definition: G4ParticleHPThermalScattering.hh:66
G4ParticleHPManager::GetThermalScatteringInelasticFinalStates
std::map< G4int, std::map< G4double, std::vector< E_P_E_isoAng * > * > * > * GetThermalScatteringInelasticFinalStates()
Definition: G4ParticleHPManager.hh:122
E_isoAng::energy
G4double energy
Definition: G4ParticleHPThermalScattering.hh:54
second
static constexpr double second
Definition: G4SIunits.hh:157
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:89
G4Material::GetElement
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:202
G4int
int G4int
Definition: G4Types.hh:78
E_P_E_isoAng::vE_isoAngle
std::vector< E_isoAng * > vE_isoAngle
Definition: G4ParticleHPThermalScattering.hh:68
G4ParticleHPManager::GetThermalScatteringIncoherentFinalStates
std::map< G4int, std::map< G4double, std::vector< E_isoAng * > * > * > * GetThermalScatteringIncoherentFinalStates()
Definition: G4ParticleHPManager.hh:120
G4ParticleHPThermalScattering::ModelDescription
virtual void ModelDescription(std::ostream &outFile) const
Definition: G4ParticleHPThermalScattering.cc:1186
G4ParticleHPThermalScatteringData::GetCoherentCrossSection
G4double GetCoherentCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
Definition: G4ParticleHPThermalScatteringData.cc:481
G4ParticleHPManager::GetThermalScatteringCoherentFinalStates
std::map< G4int, std::map< G4double, std::vector< std::pair< G4double, G4double > * > * > * > * GetThermalScatteringCoherentFinalStates()
Definition: G4ParticleHPManager.hh:118
G4HadronicInteraction::SetMinEnergy
void SetMinEnergy(G4double anEnergy)
Definition: G4HadronicInteraction.hh:90
dE
static const G4double dE
Definition: G4PhotoNuclearCrossSection.cc:47
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
G4ElementTable.hh
E_P_E_isoAng
Definition: G4ParticleHPThermalScattering.hh:63
G4Element::GetNumberOfElements
static size_t GetNumberOfElements()
Definition: G4Element.cc:405
G4HadProjectile::GetDefinition
const G4ParticleDefinition * GetDefinition() const
Definition: G4HadProjectile.hh:81
G4Threading.hh
G4Neutron.hh
G4bool
bool G4bool
Definition: G4Types.hh:79
G4HadProjectile::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4HadProjectile.hh:106
int
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
G4ParticleHPThermalScattering::BuildPhysicsTable
void BuildPhysicsTable(const G4ParticleDefinition &)
Definition: G4ParticleHPThermalScattering.cc:202
eV
static constexpr double eV
Definition: G4SIunits.hh:215
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4Material::GetNumberOfMaterials
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:594
kelvin
static constexpr double kelvin
Definition: G4SIunits.hh:281
G4HadProjectile::Get4Momentum
const G4LorentzVector & Get4Momentum() const
Definition: G4HadProjectile.hh:86
G4ParticleHPThermalScatteringData.hh
G4HadFinalState::SetEnergyChange
void SetEnergyChange(G4double anEnergy)
Definition: G4HadFinalState.cc:39
G4INCL::CrossSections::total
G4double total(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:59
G4ParticleHPThermalScattering::AddUserThermalScatteringFile
void AddUserThermalScatteringFile(G4String, G4String)
Definition: G4ParticleHPThermalScattering.cc:1164
E_isoAng::n
G4int n
Definition: G4ParticleHPThermalScattering.hh:55
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition: G4INCLKinematicsUtils.cc:157
G4ParticleHPManager::RegisterThermalScatteringInelasticFinalStates
void RegisterThermalScatteringInelasticFinalStates(std::map< G4int, std::map< G4double, std::vector< E_P_E_isoAng * > * > * > *val)
Definition: G4ParticleHPManager.hh:123
G4ParticleHPThermalScattering.hh
G4Nucleus::GetZ_asInt
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4ParticleHPManager::RegisterThermalScatteringIncoherentFinalStates
void RegisterThermalScatteringIncoherentFinalStates(std::map< G4int, std::map< G4double, std::vector< E_isoAng * > * > * > *val)
Definition: G4ParticleHPManager.hh:121
LL
static const G4int LL[nN]
Definition: G4ElectroNuclearCrossSection.cc:113
D
double D(double temp)
Definition: G4DNAElectronHoleRecombination.cc:91
G4HadronicInteraction::SetMaxEnergy
void SetMaxEnergy(const G4double anEnergy)
Definition: G4HadronicInteraction.hh:103
G4ParticleHPThermalScattering::ApplyYourself
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus)
Definition: G4ParticleHPThermalScattering.cc:392
G4ParticleHPElastic.hh
G4Material::GetTemperature
G4double GetTemperature() const
Definition: G4Material.hh:182
G4endl
#define G4endl
Definition: G4ios.hh:61
G4HadProjectile::GetMaterial
const G4Material * GetMaterial() const
Definition: G4HadProjectile.hh:76
G4MaterialTable.hh
G4Threading::IsMasterThread
G4bool IsMasterThread()
Definition: G4Threading.cc:146
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
G4ParticleHPThermalScatteringData::GetCrossSection
G4double GetCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
Definition: G4ParticleHPThermalScatteringData.cc:450
E_P_E_isoAng::sum_of_probXdEs
G4double sum_of_probXdEs
Definition: G4ParticleHPThermalScattering.hh:69
G4ElementTable
std::vector< G4Element * > G4ElementTable
Definition: G4ElementTable.hh:41
G4double
double G4double
Definition: G4Types.hh:76
G4ParticleHPThermalScatteringData::GetInelasticCrossSection
G4double GetInelasticCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
Definition: G4ParticleHPThermalScatteringData.cc:472
G4Element::GetName
const G4String & GetName() const
Definition: G4Element.hh:127
G4SystemOfUnits.hh
G4Element::GetElementTable
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:398
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
G4ParticleHPThermalScattering::GetFatalEnergyCheckLevels
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels() const
Definition: G4ParticleHPThermalScattering.cc:1158
G4ParticleHPThermalScatteringData::BuildPhysicsTable
void BuildPhysicsTable(const G4ParticleDefinition &)
Definition: G4ParticleHPThermalScatteringData.cc:215
G4HadFinalState::SetMomentumChange
void SetMomentumChange(const G4ThreeVector &aV)
Definition: G4HadFinalState.hh:56
E_P_E_isoAng::prob
std::vector< G4double > prob
Definition: G4ParticleHPThermalScattering.hh:67
G4ParticleHPManager::RegisterThermalScatteringCoherentFinalStates
void RegisterThermalScatteringCoherentFinalStates(std::map< G4int, std::map< G4double, std::vector< std::pair< G4double, G4double > * > * > * > *val)
Definition: G4ParticleHPManager.hh:119
G4ParticleHPManager.hh