#include <G4NeutronHPThermalScattering.hh>
Inheritance diagram for G4NeutronHPThermalScattering:
Public Member Functions | |
G4NeutronHPThermalScattering () | |
~G4NeutronHPThermalScattering () | |
G4HadFinalState * | ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus) |
virtual const std::pair< G4double, G4double > | GetFatalEnergyCheckLevels () const |
Definition at line 64 of file G4NeutronHPThermalScattering.hh.
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 }
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 }