G4NeutronHPThermalScattering Class Reference

#include <G4NeutronHPThermalScattering.hh>

Inheritance diagram for G4NeutronHPThermalScattering:

G4HadronicInteraction

Public Member Functions

 G4NeutronHPThermalScattering ()
 ~G4NeutronHPThermalScattering ()
G4HadFinalState * ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus)
virtual const std::pair< G4double,
G4double > 		GetFatalEnergyCheckLevels () const

Detailed Description

Definition at line 64 of file G4NeutronHPThermalScattering.hh.

Constructor & Destructor Documentation

G4NeutronHPThermalScattering::G4NeutronHPThermalScattering (  ) 

Definition at line 49 of file G4NeutronHPThermalScattering.cc.

References G4NeutronHPThermalScatteringData::BuildPhysicsTable(), G4Neutron::Neutron(), G4HadronicInteraction::SetMaxEnergy(), and G4HadronicInteraction::SetMinEnergy().

00050                              :G4HadronicInteraction("NeutronHPThermalScattering")
00051 {
00052 
00053    theHPElastic = new G4NeutronHPElastic();
00054 
00055    SetMinEnergy( 0.*eV );
00056    SetMaxEnergy( 4*eV );
00057    theXSection = new G4NeutronHPThermalScatteringData();
00058    theXSection->BuildPhysicsTable( *(G4Neutron::Neutron()) );
00059 
00060    buildPhysicsTable();
00061 
00062 }

G4NeutronHPThermalScattering::~G4NeutronHPThermalScattering (  ) 

Definition at line 66 of file G4NeutronHPThermalScattering.cc.

00067 {
00068 
00069    for ( std::map < G4int , std::map < G4double , std::vector < E_isoAng* >* >* >::iterator it = incoherentFSs.begin() ; it != incoherentFSs.end() ; it++ )
00070    {
00071       std::map < G4double , std::vector < E_isoAng* >* >::iterator itt;
00072       for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
00073       {
00074          std::vector< E_isoAng* >::iterator ittt;
00075          for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
00076          {
00077             delete *ittt;
00078          }
00079          delete itt->second;
00080       }
00081       delete it->second;
00082    }
00083 
00084    for ( std::map < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* >::iterator it = coherentFSs.begin() ; it != coherentFSs.end() ; it++ )
00085    {
00086       std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >::iterator itt;
00087       for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
00088       {
00089          std::vector < std::pair< G4double , G4double >* >::iterator ittt;
00090          for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
00091          {
00092             delete *ittt;
00093          }
00094          delete itt->second;
00095       }
00096       delete it->second;
00097    }
00098 
00099    for ( std::map < G4int ,  std::map < G4double , std::vector < E_P_E_isoAng* >* >* >::iterator it = inelasticFSs.begin() ; it != inelasticFSs.end() ; it++ )
00100    {
00101       std::map < G4double , std::vector < E_P_E_isoAng* >* >::iterator itt;
00102       for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
00103       {
00104          std::vector < E_P_E_isoAng* >::iterator ittt;
00105          for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
00106          {
00107             std::vector < E_isoAng* >::iterator it4;
00108             for ( it4 = (*ittt)->vE_isoAngle.begin() ; it4 != (*ittt)->vE_isoAngle.end() ; it4++ )
00109             {
00110                delete *it4;
00111             }
00112             delete *ittt;
00113          }
00114          delete itt->second;
00115       }
00116       delete it->second;
00117    }
00118 
00119    delete theHPElastic;
00120    delete theXSection;
00121 }

Member Function Documentation

G4HadFinalState * G4NeutronHPThermalScattering::ApplyYourself ( const G4HadProjectile &  aTrack,
G4Nucleus &  aTargetNucleus 
) [virtual]

Implements G4HadronicInteraction.

Definition at line 302 of file G4NeutronHPThermalScattering.cc.

References E_isoAng::energy, G4UniformRand, G4HadProjectile::Get4Momentum(), G4NeutronHPThermalScatteringData::GetCoherentCrossSection(), G4NeutronHPThermalScatteringData::GetCrossSection(), G4HadProjectile::GetDefinition(), G4Material::GetElement(), G4NeutronHPThermalScatteringData::GetInelasticCrossSection(), G4HadProjectile::GetKineticEnergy(), G4HadProjectile::GetMaterial(), G4Material::GetNumberOfElements(), G4Material::GetTemperature(), G4Element::GetZ(), G4Nucleus::GetZ_asInt(), E_isoAng::isoAngle, E_isoAng::n, G4HadFinalState::SetEnergyChange(), G4HadFinalState::SetMomentumChange(), and G4HadronicInteraction::theParticleChange.

00303 {
00304 
00305 // Select Element > Reaction >
00306 
00307    const G4Material * theMaterial = aTrack.GetMaterial();
00308    G4double aTemp = theMaterial->GetTemperature();
00309    G4int n = theMaterial->GetNumberOfElements();
00310    //static const G4ElementTable* theElementTable = G4Element::GetElementTable();
00311 
00312    G4bool findThermalElement = false;
00313    G4int ielement;
00314    const G4Element* theElement = NULL;
00315    for ( G4int i = 0; i < n ; i++ )
00316    {
00317       theElement = theMaterial->GetElement(i);
00318       //Select target element 
00319       if ( aNucleus.GetZ_asInt() == (G4int)(theElement->GetZ() + 0.5 ) )
00320       {
00321          //Check Applicability of Thermal Scattering 
00322          if (  getTS_ID( NULL , theElement ) != -1 )
00323          {
00324             ielement = getTS_ID( NULL , theElement );
00325             findThermalElement = true;
00326             break;
00327          }
00328          else if (  getTS_ID( theMaterial , theElement ) != -1 )
00329          {
00330             ielement = getTS_ID( theMaterial , theElement );
00331             findThermalElement = true;
00332             break;
00333          }
00334       }       
00335    } 
00336 
00337    if ( findThermalElement == true )
00338    {
00339 
00340 //    Select Reaction  (Inelastic, coherent, incoherent)  
00341 
00342       G4ParticleDefinition* pd = const_cast< G4ParticleDefinition* >( aTrack.GetDefinition() );
00343       G4DynamicParticle* dp = new G4DynamicParticle ( pd , aTrack.Get4Momentum() );
00344       G4double total = theXSection->GetCrossSection( dp , theElement , theMaterial );
00345       G4double inelastic = theXSection->GetInelasticCrossSection( dp , theElement , theMaterial );
00346    
00347 
00348       G4double random = G4UniformRand();
00349       if ( random <= inelastic/total ) 
00350       {
00351          // Inelastic
00352 
00353          // T_L and T_H 
00354          std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator it; 
00355          std::vector<G4double> v_temp;
00356          v_temp.clear();
00357          for ( it = inelasticFSs.find( ielement )->second->begin() ; it != inelasticFSs.find( ielement )->second->end() ; it++ )
00358          {
00359             v_temp.push_back( it->first );
00360          }
00361 
00362 //                   T_L         T_H 
00363          std::pair < G4double , G4double > tempLH = find_LH ( aTemp , &v_temp );
00364 //
00365 //       For T_L aNEP_EPM_TL  and T_H aNEP_EPM_TH
00366 //
00367          std::vector< E_P_E_isoAng* >* vNEP_EPM_TL = 0;
00368          std::vector< E_P_E_isoAng* >* vNEP_EPM_TH = 0;
00369 
00370          if ( tempLH.first != 0.0 && tempLH.second != 0.0 ) 
00371          {
00372             vNEP_EPM_TL = inelasticFSs.find( ielement )->second->find ( tempLH.first/kelvin )->second;
00373             vNEP_EPM_TH = inelasticFSs.find( ielement )->second->find ( tempLH.second/kelvin )->second;
00374          }
00375          else if ( tempLH.first == 0.0 )
00376          {
00377             std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator itm;  
00378             itm = inelasticFSs.find( ielement )->second->begin();
00379             vNEP_EPM_TL = itm->second;
00380             itm++;
00381             vNEP_EPM_TH = itm->second;
00382          }
00383          else if (  tempLH.second == 0.0 )
00384          {
00385             std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator itm;  
00386             itm = inelasticFSs.find( ielement )->second->end();
00387             itm--;
00388             vNEP_EPM_TH = itm->second;
00389             itm--;
00390             vNEP_EPM_TL = itm->second;
00391          } 
00392 
00393 //
00394 
00395          G4double rand_for_sE = G4UniformRand();
00396 
00397          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 );
00398          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 );
00399 
00400          G4double sE;
00401          sE = get_linear_interpolated ( aTemp , std::pair < G4double , G4double > ( tempLH.first , TL.first ) , std::pair < G4double , G4double > ( tempLH.second , TH.first ) );  
00402          E_isoAng anE_isoAng; 
00403          if ( TL.second.n == TH.second.n ) 
00404          {
00405             anE_isoAng.energy = sE; 
00406             anE_isoAng.n =  TL.second.n; 
00407             for ( G4int i=0 ; i < anE_isoAng.n ; i++ )
00408             { 
00409                G4double angle;
00410                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 ] ) );  
00411                anE_isoAng.isoAngle.push_back( angle ); 
00412             }
00413          }
00414          else
00415          {
00416             std::cout << "Do not Suuport yet." << std::endl; 
00417          }
00418      
00419          //set 
00420          theParticleChange.SetEnergyChange( sE );
00421          G4double mu = getMu( &anE_isoAng );
00422          theParticleChange.SetMomentumChange( 0.0 , std::sqrt ( 1 - mu*mu ) , mu );
00423 
00424       } 
00425       //else if ( random <= ( inelastic + theXSection->GetCoherentCrossSection( dp , (*theElementTable)[ ielement ] , aTemp ) ) / total )
00426       else if ( random <= ( inelastic + theXSection->GetCoherentCrossSection( dp , theElement , theMaterial ) ) / total )
00427       {
00428          // Coherent Elastic 
00429 
00430          G4double E = aTrack.GetKineticEnergy();
00431 
00432          // T_L and T_H 
00433          std::map < G4double , std::vector< std::pair< G4double , G4double >* >* >::iterator it; 
00434          std::vector<G4double> v_temp;
00435          v_temp.clear();
00436          for ( it = coherentFSs.find( ielement )->second->begin() ; it != coherentFSs.find( ielement )->second->end() ; it++ )
00437          {
00438             v_temp.push_back( it->first );
00439          }
00440 
00441 //                   T_L         T_H 
00442          std::pair < G4double , G4double > tempLH = find_LH ( aTemp , &v_temp );
00443 //
00444 //
00445 //       For T_L anEPM_TL  and T_H anEPM_TH
00446 //
00447          std::vector< std::pair< G4double , G4double >* >* pvE_p_TL = 0; 
00448          std::vector< std::pair< G4double , G4double >* >* pvE_p_TH = 0; 
00449 
00450          if ( tempLH.first != 0.0 && tempLH.second != 0.0 ) 
00451          {
00452             pvE_p_TL = coherentFSs.find( ielement )->second->find ( tempLH.first/kelvin )->second;
00453             pvE_p_TH = coherentFSs.find( ielement )->second->find ( tempLH.first/kelvin )->second;
00454          }
00455          else if ( tempLH.first == 0.0 )
00456          {
00457             pvE_p_TL = coherentFSs.find( ielement )->second->find ( v_temp[ 0 ] )->second;
00458             pvE_p_TH = coherentFSs.find( ielement )->second->find ( v_temp[ 1 ] )->second;
00459          }
00460          else if (  tempLH.second == 0.0 )
00461          {
00462             pvE_p_TL = coherentFSs.find( ielement )->second->find ( v_temp.back() )->second;
00463             std::vector< G4double >::iterator itv;
00464             itv = v_temp.end();
00465             itv--;
00466             itv--;
00467             pvE_p_TL = coherentFSs.find( ielement )->second->find ( *itv )->second;
00468          }
00469 
00470 
00471          std::vector< G4double > vE_T;
00472          std::vector< G4double > vp_T;
00473 
00474          G4int n1 = pvE_p_TL->size();  
00475          //G4int n2 = pvE_p_TH->size();  
00476 
00477          for ( G4int i=1 ; i < n1 ; i++ ) 
00478          {
00479             if ( (*pvE_p_TL)[i]->first != (*pvE_p_TH)[i]->first ) G4HadronicException(__FILE__, __LINE__, "A problem is found in Thermal Scattering Data!");
00480             vE_T.push_back ( (*pvE_p_TL)[i]->first );
00481             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 ) ) );  
00482          }
00483 
00484          G4int j = 0;  
00485          for ( G4int i = 1 ; i < n ; i++ ) 
00486          {
00487             if ( E/eV < vE_T[ i ] ) 
00488             {
00489                j = i-1;
00490                break;
00491             }
00492          }
00493 
00494          G4double rand_for_mu = G4UniformRand();
00495 
00496          G4int k = 0;
00497          for ( G4int i = 1 ; i < j ; i++ )
00498          {
00499              G4double Pi = vp_T[ i ] / vp_T[ j ]; 
00500              if ( rand_for_mu < Pi )
00501              {
00502                 k = i-1; 
00503                 break;
00504              }
00505          }
00506 
00507          //G4double Ei = vE_T[ j ];
00508          G4double Ei = vE_T[ k ];
00509 
00510          G4double mu = 1 - 2 * Ei / (E/eV) ;  
00511          //111102
00512          if ( mu < -1.0 ) mu = -1.0;
00513 
00514          theParticleChange.SetEnergyChange( E );
00515          theParticleChange.SetMomentumChange( 0.0 , std::sqrt ( 1 - mu*mu ) , mu );
00516 
00517 
00518       }
00519       else
00520       {
00521          // InCoherent Elastic
00522 
00523          // T_L and T_H 
00524          std::map < G4double , std::vector < E_isoAng* >* >::iterator it; 
00525          std::vector<G4double> v_temp;
00526          v_temp.clear();
00527          for ( it = incoherentFSs.find( ielement )->second->begin() ; it != incoherentFSs.find( ielement )->second->end() ; it++ )
00528          {
00529             v_temp.push_back( it->first );
00530          }
00531               
00532 //                   T_L         T_H 
00533          std::pair < G4double , G4double > tempLH = find_LH ( aTemp , &v_temp );
00534 
00535 //
00536 //       For T_L anEPM_TL  and T_H anEPM_TH
00537 //
00538 
00539          E_isoAng anEPM_TL_E;
00540          E_isoAng anEPM_TH_E;
00541 
00542          if ( tempLH.first != 0.0 && tempLH.second != 0.0 ) 
00543          {
00544             anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( tempLH.first/kelvin )->second );
00545             anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( tempLH.second/kelvin )->second );
00546          }
00547          else if ( tempLH.first == 0.0 )
00548          {
00549             anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( v_temp[ 0 ] )->second );
00550             anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( v_temp[ 1 ] )->second );
00551          }
00552          else if (  tempLH.second == 0.0 )
00553          {
00554             anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( v_temp.back() )->second );
00555             std::vector< G4double >::iterator itv;
00556             itv = v_temp.end();
00557             itv--;
00558             itv--;
00559             anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( *itv )->second );
00560          } 
00561         
00562          // E_isoAng for aTemp and aTrack.GetKineticEnergy() 
00563          E_isoAng anEPM_T_E;  
00564 
00565          if ( anEPM_TL_E.n == anEPM_TH_E.n ) 
00566          {
00567             anEPM_T_E.n = anEPM_TL_E.n; 
00568             for ( G4int i=0 ; i < anEPM_TL_E.n ; i++ )
00569             { 
00570                G4double angle;
00571                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 ] ) );  
00572                anEPM_T_E.isoAngle.push_back( angle ); 
00573             }
00574          }
00575          else
00576          {
00577             std::cout << "Do not Suuport yet." << std::endl; 
00578          }
00579 
00580          // Decide mu 
00581          G4double mu = getMu ( &anEPM_T_E );
00582 
00583          // Set Final State
00584          theParticleChange.SetEnergyChange( aTrack.GetKineticEnergy() );  // No energy change in Elastic
00585          theParticleChange.SetMomentumChange( 0.0 , std::sqrt ( 1 - mu*mu ) , mu ); 
00586 
00587       } 
00588       delete dp;
00589 
00590       return &theParticleChange;
00591       
00592    }
00593    else 
00594    {
00595       // Not thermal element   
00596       // Neutron HP will handle
00597       return theHPElastic -> ApplyYourself( aTrack, aNucleus );
00598    }
00599 
00600 }

const std::pair< G4double, G4double > G4NeutronHPThermalScattering::GetFatalEnergyCheckLevels (  )  const [virtual]

Reimplemented from G4HadronicInteraction.

Definition at line 999 of file G4NeutronHPThermalScattering.cc.

References DBL_MAX.

01000 {
01001    //return std::pair<G4double, G4double>(10*perCent,10*GeV);
01002    return std::pair<G4double, G4double>(10*perCent,DBL_MAX);
01003 }

The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:43 2013 for Geant4 by  doxygen 1.4.7