Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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 //
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 
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 
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 
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);
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);
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;
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);
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;
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 
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
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 
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  //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 
753 std::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 
790 G4double 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 
808 E_isoAng G4ParticleHPThermalScattering::create_E_isoAng_from_energy ( G4double energy , std::vector< E_isoAng* >* vEPM )
809 {
810  E_isoAng anEPM_T_E;
811 
812  std::vector< E_isoAng* >::iterator iv;
813 
814  std::vector< G4double > v_e;
815  v_e.clear();
816  for ( iv = vEPM->begin() ; iv != vEPM->end() ; iv++ )
817  v_e.push_back ( (*iv)->energy );
818 
819  std::pair < G4double , G4double > energyLH = find_LH ( energy , &v_e );
820  //G4cout << " " << energy/eV << " " << energyLH.first/eV << " " << energyLH.second/eV << G4endl;
821 
822  E_isoAng* panEPM_T_EL=0;
823  E_isoAng* panEPM_T_EH=0;
824 
825  if ( energyLH.first != 0.0 && energyLH.second != 0.0 )
826  {
827  for ( iv = vEPM->begin() ; iv != vEPM->end() ; iv++ )
828  {
829  if ( energyLH.first == (*iv)->energy ) {
830  panEPM_T_EL = *iv;
831  iv++;
832  panEPM_T_EH = *iv;
833  break;
834  }
835  }
836  }
837  else if ( energyLH.first == 0.0 )
838  {
839  panEPM_T_EL = (*vEPM)[0];
840  panEPM_T_EH = (*vEPM)[1];
841  }
842  else if ( energyLH.second == 0.0 )
843  {
844  panEPM_T_EH = (*vEPM).back();
845  iv = vEPM->end();
846  iv--;
847  iv--;
848  panEPM_T_EL = *iv;
849  }
850 
851  //checking isoAng has proper values or not
852  // Inelastic/FS, the first and last entries of *vEPM has all zero values.
853  if ( ! ( check_E_isoAng (panEPM_T_EL) ) ) panEPM_T_EL= panEPM_T_EH;
854  if ( ! ( check_E_isoAng (panEPM_T_EH) ) ) panEPM_T_EH= panEPM_T_EL;
855 
856  if ( panEPM_T_EL->n == panEPM_T_EH->n )
857  {
858  anEPM_T_E.energy = energy;
859  anEPM_T_E.n = panEPM_T_EL->n;
860 
861  for ( G4int i=0 ; i < panEPM_T_EL->n ; i++ )
862  {
863  G4double angle;
864  angle = get_linear_interpolated ( energy , std::pair< G4double , G4double > ( energyLH.first , panEPM_T_EL->isoAngle[ i ] ) , std::pair< G4double , G4double > ( energyLH.second , panEPM_T_EH->isoAngle[ i ] ) );
865  anEPM_T_E.isoAngle.push_back( angle );
866  }
867  }
868  else
869  {
870  G4cout << "G4ParticleHPThermalScattering Do not Suuport yet." << G4endl;
871  }
872 
873 
874  return anEPM_T_E;
875 }
876 
877 
878 
879 G4double G4ParticleHPThermalScattering::get_secondary_energy_from_E_P_E_isoAng ( G4double random , E_P_E_isoAng* anE_P_E_isoAng )
880 {
881 
882  G4double secondary_energy = 0.0;
883 
884  G4int n = anE_P_E_isoAng->n;
885  G4double sum_p = 0.0; // sum_p_H
886  G4double sum_p_L = 0.0;
887 
888  G4double total=0.0;
889 
890 /*
891  delete for speed up
892  for ( G4int i = 0 ; i < n-1 ; i++ )
893  {
894  G4double E_L = anE_P_E_isoAng->vE_isoAngle[i]->energy/eV;
895  G4double E_H = anE_P_E_isoAng->vE_isoAngle[i+1]->energy/eV;
896  G4double dE = E_H - E_L;
897  total += ( ( anE_P_E_isoAng->prob[i] ) * dE );
898  }
899 
900  if ( std::abs( total - anE_P_E_isoAng->sum_of_probXdEs ) > 1.0e-14 ) G4cout << total - anE_P_E_isoAng->sum_of_probXdEs << G4endl;
901 */
902  total = anE_P_E_isoAng->sum_of_probXdEs;
903 
904  for ( G4int i = 0 ; i < n-1 ; i++ )
905  {
906  G4double E_L = anE_P_E_isoAng->vE_isoAngle[i]->energy/eV;
907  G4double E_H = anE_P_E_isoAng->vE_isoAngle[i+1]->energy/eV;
908  G4double dE = E_H - E_L;
909  sum_p += ( ( anE_P_E_isoAng->prob[i] ) * dE );
910 
911  if ( random <= sum_p/total )
912  {
913  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 ) );
914  secondary_energy = secondary_energy*eV; //need eV
915  break;
916  }
917  sum_p_L = sum_p;
918  }
919 
920  return secondary_energy;
921 }
922 
923 
924 
925 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 )
926 {
927 
928  std::map< G4double , G4int > map_energy;
929  map_energy.clear();
930  std::vector< G4double > v_energy;
931  v_energy.clear();
932  std::vector< E_P_E_isoAng* >::iterator itv;
933  G4int i = 0;
934  for ( itv = vNEP_EPM->begin(); itv != vNEP_EPM->end(); itv++ )
935  {
936  v_energy.push_back( (*itv)->energy );
937  map_energy.insert( std::pair < G4double , G4int > ( (*itv)->energy , i ) );
938  i++;
939  }
940 
941  std::pair < G4double , G4double > energyLH = find_LH ( pE , &v_energy );
942 
943  E_P_E_isoAng* pE_P_E_isoAng_EL = 0;
944  E_P_E_isoAng* pE_P_E_isoAng_EH = 0;
945 
946  if ( energyLH.first != 0.0 && energyLH.second != 0.0 )
947  {
948  pE_P_E_isoAng_EL = (*vNEP_EPM)[ map_energy.find ( energyLH.first )->second ];
949  pE_P_E_isoAng_EH = (*vNEP_EPM)[ map_energy.find ( energyLH.second )->second ];
950  }
951  else if ( energyLH.first == 0.0 )
952  {
953  pE_P_E_isoAng_EL = (*vNEP_EPM)[ 0 ];
954  pE_P_E_isoAng_EH = (*vNEP_EPM)[ 1 ];
955  }
956  if ( energyLH.second == 0.0 )
957  {
958  pE_P_E_isoAng_EH = (*vNEP_EPM).back();
959  itv = vNEP_EPM->end();
960  itv--;
961  itv--;
962  pE_P_E_isoAng_EL = *itv;
963  }
964 
965 
966  G4double sE;
967  G4double sE_L;
968  G4double sE_H;
969 
970 
971  sE_L = get_secondary_energy_from_E_P_E_isoAng ( rand_for_sE , pE_P_E_isoAng_EL );
972  sE_H = get_secondary_energy_from_E_P_E_isoAng ( rand_for_sE , pE_P_E_isoAng_EH );
973 
974  sE = get_linear_interpolated ( pE , std::pair < G4double , G4double > ( energyLH.first , sE_L ) , std::pair < G4double , G4double > ( energyLH.second , sE_H ) );
975 
976 
977  E_isoAng E_isoAng_L = create_E_isoAng_from_energy ( sE , &(pE_P_E_isoAng_EL->vE_isoAngle) );
978  E_isoAng E_isoAng_H = create_E_isoAng_from_energy ( sE , &(pE_P_E_isoAng_EH->vE_isoAngle) );
979 
980  E_isoAng anE_isoAng;
981  //For defeating warning message from compiler
982  anE_isoAng.n = 1;
983  anE_isoAng.energy = sE; //never used
984  if ( E_isoAng_L.n == E_isoAng_H.n )
985  {
986  anE_isoAng.n = E_isoAng_L.n;
987  for ( G4int j=0 ; j < anE_isoAng.n ; j++ )
988  {
989  G4double angle;
990  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 ] ) );
991  anE_isoAng.isoAngle.push_back( angle );
992  }
993  }
994  else
995  {
996  //G4cout << "Do not Suuport yet." << G4endl;
997  throw G4HadronicException(__FILE__, __LINE__, "Unexpected values!");
998  }
999 
1000 
1001 
1002  return std::pair< G4double , E_isoAng >( sE , anE_isoAng);
1003 }
1004 
1005 void G4ParticleHPThermalScattering::buildPhysicsTable()
1006 {
1007 
1008  //Is rebuild of physics table a necessity
1009  if ( nMaterial == G4Material::GetMaterialTable()->size() && nElement == G4Element::GetElementTable()->size() ) {
1010  return;
1011  } else {
1012  nMaterial = G4Material::GetMaterialTable()->size();
1013  nElement = G4Element::GetElementTable()->size();
1014  }
1015 
1016  dic.clear();
1017  std::map < G4String , G4int > co_dic;
1018 
1019  //Searching Nist Materials
1020  static G4ThreadLocal G4MaterialTable* theMaterialTable = 0 ; if (!theMaterialTable) theMaterialTable= G4Material::GetMaterialTable();
1021  size_t numberOfMaterials = G4Material::GetNumberOfMaterials();
1022  for ( size_t i = 0 ; i < numberOfMaterials ; i++ )
1023  {
1024  G4Material* material = (*theMaterialTable)[i];
1025  size_t numberOfElements = material->GetNumberOfElements();
1026  for ( size_t j = 0 ; j < numberOfElements ; j++ )
1027  {
1028  const G4Element* element = material->GetElement(j);
1029  if ( names.IsThisThermalElement ( material->GetName() , element->GetName() ) )
1030  {
1031  G4int ts_ID_of_this_geometry;
1032  G4String ts_ndl_name = names.GetTS_NDL_Name( material->GetName() , element->GetName() );
1033  if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
1034  {
1035  ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) -> second;
1036  }
1037  else
1038  {
1039  ts_ID_of_this_geometry = co_dic.size();
1040  co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
1041  }
1042 
1043  //G4cout << "Neutron HP Thermal Scattering: Registering a material-element pair of "
1044  // << material->GetName() << " " << element->GetName()
1045  // << " as internal thermal scattering id of " << ts_ID_of_this_geometry << "." << G4endl;
1046 
1047  dic.insert( std::pair < std::pair < G4Material* , const G4Element* > , G4int > ( std::pair < G4Material* , const G4Element* > ( material , element ) , ts_ID_of_this_geometry ) );
1048  }
1049  }
1050  }
1051 
1052  //Searching TS Elements
1053  static G4ThreadLocal G4ElementTable* theElementTable = 0 ; if (!theElementTable) theElementTable= G4Element::GetElementTable();
1054  size_t numberOfElements = G4Element::GetNumberOfElements();
1055  //size_t numberOfThermalElements = 0;
1056  for ( size_t i = 0 ; i < numberOfElements ; i++ )
1057  {
1058  const G4Element* element = (*theElementTable)[i];
1059  if ( names.IsThisThermalElement ( element->GetName() ) )
1060  {
1061  if ( names.IsThisThermalElement ( element->GetName() ) )
1062  {
1063  G4int ts_ID_of_this_geometry;
1064  G4String ts_ndl_name = names.GetTS_NDL_Name( element->GetName() );
1065  if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
1066  {
1067  ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) -> second;
1068  }
1069  else
1070  {
1071  ts_ID_of_this_geometry = co_dic.size();
1072  co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
1073  }
1074 
1075  //G4cout << "Neutron HP Thermal Scattering: Registering an element of "
1076  // << material->GetName() << " " << element->GetName()
1077  // << " as internal thermal scattering id of " << ts_ID_of_this_geometry << "." << G4endl;
1078 
1079  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 ) );
1080  }
1081  }
1082  }
1083 
1084  G4cout << G4endl;
1085  G4cout << "Neutron HP Thermal Scattering: Following material-element pairs or elements are registered." << G4endl;
1086  for ( std::map < std::pair < const G4Material* , const G4Element* > , G4int >::iterator it = dic.begin() ; it != dic.end() ; it++ )
1087  {
1088  if ( it->first.first != NULL )
1089  {
1090  G4cout << "Material " << it->first.first->GetName() << " - Element " << it->first.second->GetName() << ", internal thermal scattering id " << it->second << G4endl;
1091  }
1092  else
1093  {
1094  G4cout << "Element " << it->first.second->GetName() << ", internal thermal scattering id " << it->second << G4endl;
1095  }
1096  }
1097  G4cout << G4endl;
1098 
1099  // Read Cross Section Data files
1100 
1102  coherentFSs = hpmanager->GetThermalScatteringCoherentFinalStates();
1103  incoherentFSs = hpmanager->GetThermalScatteringIncoherentFinalStates();
1104  inelasticFSs = hpmanager->GetThermalScatteringInelasticFinalStates();
1105 
1106  if ( G4Threading::IsMasterThread() ) {
1107 
1108  clearCurrentFSData();
1109 
1110  if ( coherentFSs == NULL ) coherentFSs = new std::map < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* >;
1111  if ( incoherentFSs == NULL ) incoherentFSs = new std::map < G4int , std::map < G4double , std::vector < E_isoAng* >* >* >;
1112  if ( inelasticFSs == NULL ) inelasticFSs = new std::map < G4int , std::map < G4double , std::vector < E_P_E_isoAng* >* >* >;
1113 
1114  G4String dirName;
1115  if ( !getenv( "G4NEUTRONHPDATA" ) )
1116  throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
1117  dirName = getenv( "G4NEUTRONHPDATA" );
1118 
1119  //G4String name;
1120 
1121  for ( std::map < G4String , G4int >::iterator it = co_dic.begin() ; it != co_dic.end() ; it++ )
1122  {
1123  G4String tsndlName = it->first;
1124  G4int ts_ID = it->second;
1125 
1126  // Coherent
1127  G4String fsName = "/ThermalScattering/Coherent/FS/";
1128  G4String fileName = dirName + fsName + tsndlName;
1129  coherentFSs->insert ( std::pair < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* > ( ts_ID , readACoherentFSDATA( fileName ) ) );
1130 
1131  // incoherent elastic
1132  fsName = "/ThermalScattering/Incoherent/FS/";
1133  fileName = dirName + fsName + tsndlName;
1134  incoherentFSs->insert ( std::pair < G4int , std::map < G4double , std::vector < E_isoAng* >* >* > ( ts_ID , readAnIncoherentFSDATA( fileName ) ) );
1135 
1136  // inelastic
1137  fsName = "/ThermalScattering/Inelastic/FS/";
1138  fileName = dirName + fsName + tsndlName;
1139  inelasticFSs->insert ( std::pair < G4int , std::map < G4double , std::vector < E_P_E_isoAng* >* >* > ( ts_ID , readAnInelasticFSDATA( fileName ) ) );
1140  }
1141 
1142  hpmanager->RegisterThermalScatteringCoherentFinalStates( coherentFSs );
1143  hpmanager->RegisterThermalScatteringIncoherentFinalStates( incoherentFSs );
1144  hpmanager->RegisterThermalScatteringInelasticFinalStates( inelasticFSs );
1145  }
1146 
1147  theXSection->BuildPhysicsTable( *(G4Neutron::Neutron()) );
1148 }
1149 
1150 
1151 G4int G4ParticleHPThermalScattering::getTS_ID ( const G4Material* material , const G4Element* element )
1152 {
1153  G4int result = -1;
1154  if ( dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) ) != dic.end() )
1155  result = dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) )->second;
1156  return result;
1157 }
1158 
1159 const std::pair<G4double, G4double> G4ParticleHPThermalScattering::GetFatalEnergyCheckLevels() const
1160 {
1161  //return std::pair<G4double, G4double>(10*perCent,10*GeV);
1162  return std::pair<G4double, G4double>(10*perCent,DBL_MAX);
1163 }
1164 
1166 {
1167  names.AddThermalElement( nameG4Element , filename );
1168  theXSection->AddUserThermalScatteringFile( nameG4Element , filename );
1169  buildPhysicsTable();
1170 }
1171 
1172 
1173 G4bool G4ParticleHPThermalScattering::check_E_isoAng( E_isoAng* anE_IsoAng )
1174 {
1175  G4bool result=false;
1176 
1177  G4int n = anE_IsoAng->n;
1178  G4double sum=0.0;
1179  for ( G4int i = 0 ; i < n ; i++ ) {
1180  sum += anE_IsoAng->isoAngle[ i ];
1181  }
1182  if ( sum != 0.0 ) result = true;
1183 
1184  return result;
1185 }
1186 
1187 void G4ParticleHPThermalScattering::ModelDescription(std::ostream& outFile) const
1188 {
1189  outFile << "High Precision model based on thermal scattering data in evaluated nuclear data libraries for neutrons below 5eV on specific materials\n";
1190 }
G4double G4ParticleHPJENDLHEData::G4double result
static G4ParticleHPManager * GetInstance()
const XML_Char * name
Definition: expat.h:151
std::vector< G4double > isoAngle
G4int first(char) const
void BuildPhysicsTable(const G4ParticleDefinition &)
static constexpr double perCent
Definition: G4SIunits.hh:332
const G4String & GetName() const
Definition: G4Material.hh:178
G4double GetZ() const
Definition: G4Element.hh:131
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
std::vector< G4Material * > G4MaterialTable
void GetDataStream(G4String, std::istringstream &iss)
static G4double angle[DIM]
std::map< G4int, std::map< G4double, std::vector< E_P_E_isoAng * > * > * > * GetThermalScatteringInelasticFinalStates()
tuple x
Definition: test.py:50
static constexpr double second
Definition: G4SIunits.hh:157
#define G4ThreadLocal
Definition: tls.hh:89
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:202
int G4int
Definition: G4Types.hh:78
std::vector< E_isoAng * > vE_isoAngle
std::map< G4int, std::map< G4double, std::vector< E_isoAng * > * > * > * GetThermalScatteringIncoherentFinalStates()
virtual void ModelDescription(std::ostream &outFile) const
G4double GetCoherentCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
std::map< G4int, std::map< G4double, std::vector< std::pair< G4double, G4double > * > * > * > * GetThermalScatteringCoherentFinalStates()
string material
Definition: eplot.py:19
void SetMinEnergy(G4double anEnergy)
static const G4double dE
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
static size_t GetNumberOfElements()
Definition: G4Element.cc:405
const G4ParticleDefinition * GetDefinition() const
bool G4bool
Definition: G4Types.hh:79
G4double GetKineticEnergy() const
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
void BuildPhysicsTable(const G4ParticleDefinition &)
static constexpr double eV
Definition: G4SIunits.hh:215
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
const G4int n
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:594
static constexpr double kelvin
Definition: G4SIunits.hh:281
const G4LorentzVector & Get4Momentum() const
void SetEnergyChange(G4double anEnergy)
G4double total(Particle const *const p1, Particle const *const p2)
void AddUserThermalScatteringFile(G4String, G4String)
G4double energy(const ThreeVector &p, const G4double m)
void RegisterThermalScatteringInelasticFinalStates(std::map< G4int, std::map< G4double, std::vector< E_P_E_isoAng * > * > * > *val)
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
void RegisterThermalScatteringIncoherentFinalStates(std::map< G4int, std::map< G4double, std::vector< E_isoAng * > * > * > *val)
static const G4int LL[nN]
void SetMaxEnergy(const G4double anEnergy)
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus)
G4double GetTemperature() const
Definition: G4Material.hh:182
#define G4endl
Definition: G4ios.hh:61
const G4Material * GetMaterial() const
G4bool IsMasterThread()
Definition: G4Threading.cc:146
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
G4double GetCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
std::vector< G4Element * > G4ElementTable
double G4double
Definition: G4Types.hh:76
G4double GetInelasticCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
const G4String & GetName() const
Definition: G4Element.hh:127
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:398
#define DBL_MAX
Definition: templates.hh:83
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels() const
void BuildPhysicsTable(const G4ParticleDefinition &)
void SetMomentumChange(const G4ThreeVector &aV)
std::vector< G4double > prob
void RegisterThermalScatteringCoherentFinalStates(std::map< G4int, std::map< G4double, std::vector< std::pair< G4double, G4double > * > * > * > *val)