Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPThermalScattering Class Reference

#include <G4ParticleHPThermalScattering.hh>

Inheritance diagram for G4ParticleHPThermalScattering:
Collaboration diagram for G4ParticleHPThermalScattering:

Public Member Functions

 G4ParticleHPThermalScattering ()
 
 ~G4ParticleHPThermalScattering ()
 
G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus)
 
virtual const std::pair
< G4double, G4double
GetFatalEnergyCheckLevels () const
 
void AddUserThermalScatteringFile (G4String, G4String)
 
void BuildPhysicsTable (const G4ParticleDefinition &)
 
virtual void ModelDescription (std::ostream &outFile) const
 
- Public Member Functions inherited from G4HadronicInteraction
 G4HadronicInteraction (const G4String &modelName="HadronicModel")
 
virtual ~G4HadronicInteraction ()
 
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
virtual G4bool IsApplicable (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
G4double GetMinEnergy () const
 
G4double GetMinEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMinEnergy (G4double anEnergy)
 
void SetMinEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMinEnergy (G4double anEnergy, const G4Material *aMaterial)
 
G4double GetMaxEnergy () const
 
G4double GetMaxEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMaxEnergy (const G4double anEnergy)
 
void SetMaxEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMaxEnergy (G4double anEnergy, const G4Material *aMaterial)
 
G4int GetVerboseLevel () const
 
void SetVerboseLevel (G4int value)
 
const G4StringGetModelName () const
 
void DeActivateFor (const G4Material *aMaterial)
 
void ActivateFor (const G4Material *aMaterial)
 
void DeActivateFor (const G4Element *anElement)
 
void ActivateFor (const G4Element *anElement)
 
G4bool IsBlocked (const G4Material *aMaterial) const
 
G4bool IsBlocked (const G4Element *anElement) const
 
void SetRecoilEnergyThreshold (G4double val)
 
G4double GetRecoilEnergyThreshold () const
 
virtual std::pair< G4double,
G4double
GetEnergyMomentumCheckLevels () const
 
void SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
 
virtual void InitialiseModel ()
 

Additional Inherited Members

- Protected Member Functions inherited from G4HadronicInteraction
void SetModelName (const G4String &nam)
 
G4bool IsBlocked () const
 
void Block ()
 
- Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState theParticleChange
 
G4int verboseLevel
 
G4double theMinEnergy
 
G4double theMaxEnergy
 
G4bool isBlocked
 

Detailed Description

Definition at line 77 of file G4ParticleHPThermalScattering.hh.

Constructor & Destructor Documentation

G4ParticleHPThermalScattering::G4ParticleHPThermalScattering ( )

Definition at line 56 of file G4ParticleHPThermalScattering.cc.

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 }
void SetMinEnergy(G4double anEnergy)
static constexpr double eV
Definition: G4SIunits.hh:215
G4HadronicInteraction(const G4String &modelName="HadronicModel")
void SetMaxEnergy(const G4double anEnergy)

Here is the call graph for this function:

G4ParticleHPThermalScattering::~G4ParticleHPThermalScattering ( )

Definition at line 76 of file G4ParticleHPThermalScattering.cc.

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 }

Member Function Documentation

void G4ParticleHPThermalScattering::AddUserThermalScatteringFile ( G4String  nameG4Element,
G4String  filename 
)

Definition at line 1165 of file G4ParticleHPThermalScattering.cc.

1166 {
1167  names.AddThermalElement( nameG4Element , filename );
1168  theXSection->AddUserThermalScatteringFile( nameG4Element , filename );
1169  buildPhysicsTable();
1170 }

Here is the call graph for this function:

G4HadFinalState * G4ParticleHPThermalScattering::ApplyYourself ( const G4HadProjectile aTrack,
G4Nucleus aTargetNucleus 
)
virtual

Implements G4HadronicInteraction.

Definition at line 392 of file G4ParticleHPThermalScattering.cc.

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 }
std::vector< G4double > isoAngle
G4double GetZ() const
Definition: G4Element.hh:131
static G4double angle[DIM]
static constexpr double second
Definition: G4SIunits.hh:157
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:202
int G4int
Definition: G4Types.hh:78
G4double GetCoherentCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
#define G4UniformRand()
Definition: Randomize.hh:97
const G4ParticleDefinition * GetDefinition() const
bool G4bool
Definition: G4Types.hh:79
G4double GetKineticEnergy() const
static constexpr double eV
Definition: G4SIunits.hh:215
const G4int n
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)
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus)
G4double GetTemperature() const
Definition: G4Material.hh:182
const G4Material * GetMaterial() const
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
G4double GetCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
double G4double
Definition: G4Types.hh:76
G4double GetInelasticCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
void SetMomentumChange(const G4ThreeVector &aV)

Here is the call graph for this function:

void G4ParticleHPThermalScattering::BuildPhysicsTable ( const G4ParticleDefinition particle)
virtual

Reimplemented from G4HadronicInteraction.

Definition at line 202 of file G4ParticleHPThermalScattering.cc.

202  {
203  buildPhysicsTable();
204  theHPElastic->BuildPhysicsTable( particle );
205 }
void BuildPhysicsTable(const G4ParticleDefinition &)

Here is the call graph for this function:

const std::pair< G4double, G4double > G4ParticleHPThermalScattering::GetFatalEnergyCheckLevels ( ) const
virtual

Reimplemented from G4HadronicInteraction.

Definition at line 1159 of file G4ParticleHPThermalScattering.cc.

1160 {
1161  //return std::pair<G4double, G4double>(10*perCent,10*GeV);
1162  return std::pair<G4double, G4double>(10*perCent,DBL_MAX);
1163 }
static constexpr double perCent
Definition: G4SIunits.hh:332
#define DBL_MAX
Definition: templates.hh:83
void G4ParticleHPThermalScattering::ModelDescription ( std::ostream &  outFile) const
virtual

Reimplemented from G4HadronicInteraction.

Definition at line 1187 of file G4ParticleHPThermalScattering.cc.

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 }

The documentation for this class was generated from the following files: