74 std::map < G4double , std::vector < E_isoAng* >* >::iterator itt;
75 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
77 std::vector< E_isoAng* >::iterator ittt;
78 for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
87 for ( std::map <
G4int , std::map <
G4double , std::vector < std::pair< G4double , G4double >* >* >* >::iterator it =
coherentFSs.begin() ; it !=
coherentFSs.end() ; it++ )
89 std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >::iterator itt;
90 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
92 std::vector < std::pair< G4double , G4double >* >::iterator ittt;
93 for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
104 std::map < G4double , std::vector < E_P_E_isoAng* >* >::iterator itt;
105 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
107 std::vector < E_P_E_isoAng* >::iterator ittt;
108 for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
110 std::vector < E_isoAng* >::iterator it4;
111 for ( it4 = (*ittt)->vE_isoAngle.begin() ; it4 != (*ittt)->vE_isoAngle.end() ; it4++ )
131 std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* aCoherentFSDATA =
new std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >;
134 std::istringstream theChannel(std::ios::in);
137 std::vector< G4double > vBraggE;
140 while ( theChannel >> dummy )
145 std::vector < std::pair< G4double , G4double >* >* anBragE_P =
new std::vector < std::pair< G4double , G4double >* >;
149 for (
G4int i = 0 ; i <
n ; i++ )
153 if ( aCoherentFSDATA->size() == 0 )
156 vBraggE.push_back( Ei );
163 anBragE_P->push_back (
new std::pair < G4double , G4double > ( Ei , Pi ) );
166 aCoherentFSDATA->insert ( std::pair <
G4double , std::vector < std::pair< G4double , G4double >* >* > ( temp , anBragE_P ) );
169 return aCoherentFSDATA;
176 std::map < G4double , std::vector < E_P_E_isoAng* >* >* anT_E_P_E_isoAng =
new std::map < G4double , std::vector < E_P_E_isoAng* >* >;
179 std::istringstream theChannel(std::ios::in);
183 while ( theChannel >> dummy )
188 std::vector < E_P_E_isoAng* >* vE_P_E_isoAng =
new std::vector < E_P_E_isoAng* >;
191 for (
G4int i = 0 ; i <
n ; i++ )
195 anT_E_P_E_isoAng->insert ( std::pair <
G4double , std::vector < E_P_E_isoAng* >* > ( temp , vE_P_E_isoAng ) );
199 return anT_E_P_E_isoAng;
219 for (
G4int i = 0 ; i < aData->
n ; i++ )
226 anE_isoAng->
n = nl - 2;
227 anE_isoAng->
isoAngle.resize( anE_isoAng->
n );
229 aData->
prob.push_back( prob );
231 for (
G4int j = 0 ; j < anE_isoAng->
n ; j++ )
242 for (
G4int i = 0 ; i < aData->
n - 1 ; i++ )
247 total += ( ( aData->
prob[i] ) * dE );
258 std::map < G4double , std::vector < E_isoAng* >* >* T_E =
new std::map < G4double , std::vector < E_isoAng* >* >;
261 std::istringstream theChannel(std::ios::in);
265 while ( theChannel >> dummy )
270 std::vector < E_isoAng* >* vE_isoAng =
new std::vector < E_isoAng* >;
273 for (
G4int i = 0 ; i <
n ; i++ )
275 T_E->insert ( std::pair <
G4double , std::vector < E_isoAng* >* > ( temp , vE_isoAng ) );
303 for (
G4int i = 0 ; i < aData->
n ; i++ )
327 G4bool findThermalElement =
false;
330 for (
G4int i = 0; i <
n ; i++ )
337 if (
getTS_ID( NULL , theElement ) != -1 )
339 ielement =
getTS_ID( NULL , theElement );
340 findThermalElement =
true;
343 else if (
getTS_ID( theMaterial , theElement ) != -1 )
345 ielement =
getTS_ID( theMaterial , theElement );
346 findThermalElement =
true;
352 if ( findThermalElement ==
true )
364 if ( random <= inelastic/total )
369 std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator it;
370 std::vector<G4double> v_temp;
372 for ( it =
inelasticFSs.find( ielement )->second->begin() ; it !=
inelasticFSs.find( ielement )->second->end() ; it++ )
374 v_temp.push_back( it->first );
378 std::pair < G4double , G4double > tempLH =
find_LH ( aTemp , &v_temp );
382 std::vector< E_P_E_isoAng* >* vNEP_EPM_TL = 0;
383 std::vector< E_P_E_isoAng* >* vNEP_EPM_TH = 0;
385 if ( tempLH.first != 0.0 && tempLH.second != 0.0 )
387 vNEP_EPM_TL =
inelasticFSs.find( ielement )->second->find ( tempLH.first/
kelvin )->second;
388 vNEP_EPM_TH =
inelasticFSs.find( ielement )->second->find ( tempLH.second/
kelvin )->second;
390 else if ( tempLH.first == 0.0 )
392 std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator itm;
394 vNEP_EPM_TL = itm->second;
396 vNEP_EPM_TH = itm->second;
397 tempLH.first = tempLH.second;
398 tempLH.second = itm->first;
400 else if ( tempLH.second == 0.0 )
402 std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator itm;
405 vNEP_EPM_TH = itm->second;
407 vNEP_EPM_TL = itm->second;
408 tempLH.second = tempLH.first;
409 tempLH.first = itm->first;
418 sE =
get_linear_interpolated ( aTemp , std::pair < G4double , G4double > ( tempLH.first , TL.first ) , std::pair < G4double , G4double > ( tempLH.second , TH.first ) );
422 if ( TL.second.n == TH.second.n )
425 anE_isoAng.
n = TL.second.n;
426 for (
G4int i=0 ; i < anE_isoAng.
n ; i++ )
429 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 ] ) );
430 anE_isoAng.
isoAngle.push_back( angle );
432 mu =
getMu( &anE_isoAng );
436 G4HadronicException(__FILE__, __LINE__,
"A problem is found in Thermal Scattering Data! Do not yet supported");
452 std::map < G4double , std::vector< std::pair< G4double , G4double >* >* >::iterator it;
453 std::vector<G4double> v_temp;
455 for ( it =
coherentFSs.find( ielement )->second->begin() ; it !=
coherentFSs.find( ielement )->second->end() ; it++ )
457 v_temp.push_back( it->first );
461 std::pair < G4double , G4double > tempLH =
find_LH ( aTemp , &v_temp );
466 std::vector< std::pair< G4double , G4double >* >* pvE_p_TL = NULL;
467 std::vector< std::pair< G4double , G4double >* >* pvE_p_TH = NULL;
469 if ( tempLH.first != 0.0 && tempLH.second != 0.0 )
471 pvE_p_TL =
coherentFSs.find( ielement )->second->find ( tempLH.first/
kelvin )->second;
472 pvE_p_TH =
coherentFSs.find( ielement )->second->find ( tempLH.first/
kelvin )->second;
474 else if ( tempLH.first == 0.0 )
476 pvE_p_TL =
coherentFSs.find( ielement )->second->find ( v_temp[ 0 ] )->second;
477 pvE_p_TH =
coherentFSs.find( ielement )->second->find ( v_temp[ 1 ] )->second;
478 tempLH.first = tempLH.second;
479 tempLH.second = v_temp[ 1 ];
481 else if ( tempLH.second == 0.0 )
484 std::vector< G4double >::iterator itv;
488 pvE_p_TL =
coherentFSs.find( ielement )->second->find ( *itv )->second;
489 tempLH.second = tempLH.first;
495 G4HadronicException(__FILE__, __LINE__,
"A problem is found in Thermal Scattering Data! Unexpected temperature values in data");
498 std::vector< G4double > vE_T;
499 std::vector< G4double > vp_T;
501 G4int n1 = pvE_p_TL->size();
504 for (
G4int i=1 ; i < n1 ; i++ )
506 if ( (*pvE_p_TL)[i]->first != (*pvE_p_TH)[i]->first )
G4HadronicException(__FILE__, __LINE__,
"A problem is found in Thermal Scattering Data!");
507 vE_T.push_back ( (*pvE_p_TL)[i]->first );
508 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 ) ) );
512 for (
G4int i = 1 ; i <
n ; i++ )
514 if ( E/
eV < vE_T[ i ] )
524 for (
G4int i = 1 ; i < j ; i++ )
526 G4double Pi = vp_T[ i ] / vp_T[ j ];
527 if ( rand_for_mu < Pi )
539 if ( mu < -1.0 ) mu = -1.0;
551 std::map < G4double , std::vector < E_isoAng* >* >::iterator it;
552 std::vector<G4double> v_temp;
556 v_temp.push_back( it->first );
560 std::pair < G4double , G4double > tempLH =
find_LH ( aTemp , &v_temp );
569 if ( tempLH.first != 0.0 && tempLH.second != 0.0 ) {
573 }
else if ( tempLH.first == 0.0 ) {
577 tempLH.first = tempLH.second;
578 tempLH.second = v_temp[ 1 ];
579 }
else if ( tempLH.second == 0.0 ) {
582 std::vector< G4double >::iterator itv;
587 tempLH.second = tempLH.first;
595 if ( anEPM_TL_E.
n == anEPM_TH_E.
n )
597 anEPM_T_E.
n = anEPM_TL_E.
n;
598 for (
G4int i=0 ; i < anEPM_TL_E.
n ; i++ )
601 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 ] ) );
602 anEPM_T_E.
isoAngle.push_back( angle );
604 mu =
getMu ( &anEPM_T_E );
608 G4HadronicException(__FILE__, __LINE__,
"A problem is found in Thermal Scattering Data! Do not yet supported");
638 G4int in = int ( random * ( (*anEPM).n ) );
642 G4double mu_l = (*anEPM).isoAngle[ in-1 ];
643 G4double mu_h = (*anEPM).isoAngle[ in ];
644 result = ( mu_h - mu_l ) * ( random * ( (*anEPM).n ) - in ) + mu_l;
649 G4double D = ( (*anEPM).isoAngle[ 0 ] - ( -1 ) ) + ( 1 - (*anEPM).isoAngle[ (*anEPM).n - 1 ] );
650 G4double ratio = ( (*anEPM).isoAngle[ 0 ] - ( -1 ) ) / D;
654 G4double mu_h = (*anEPM).isoAngle[ 0 ];
655 result = ( mu_h - mu_l ) * x + mu_l;
659 G4double mu_l = (*anEPM).isoAngle[ (*anEPM).n - 1 ];
661 result = ( mu_h - mu_l ) * x + mu_l;
675 if ( aVector->size() == 1 ) {
676 L = aVector->front();
677 H = aVector->front();
682 for ( std::vector< G4double >::iterator
683 it = aVector->begin() ; it != aVector->end() ; it++ ) {
686 if ( it != aVector->begin() ) {
698 if ( H == 0.0 ) L = aVector->back();
701 return std::pair < G4double , G4double > (
L , H );
709 if ( High.first - Low.first != 0 ) {
710 y = ( High.second - Low.second ) / ( High.first - Low.first ) * ( x - Low.first ) + Low.second;
712 if ( High.second == Low.second ) {
715 G4cout <<
"G4NeutronHPThermalScattering liner interpolation err!!" <<
G4endl;
728 std::vector< E_isoAng* >::iterator iv;
730 std::vector< G4double > v_e;
732 for ( iv = vEPM->begin() ; iv != vEPM->end() ; iv++ )
733 v_e.push_back ( (*iv)->energy );
735 std::pair < G4double , G4double > energyLH =
find_LH ( energy , &v_e );
741 if ( energyLH.first != 0.0 && energyLH.second != 0.0 )
743 for ( iv = vEPM->begin() ; iv != vEPM->end() ; iv++ )
745 if ( energyLH.first == (*iv)->energy )
752 else if ( energyLH.first == 0.0 )
754 panEPM_T_EL = (*vEPM)[0];
755 panEPM_T_EH = (*vEPM)[1];
757 else if ( energyLH.second == 0.0 )
759 panEPM_T_EH = (*vEPM).back();
766 if ( panEPM_T_EL->
n == panEPM_T_EH->
n )
769 anEPM_T_E.
n = panEPM_T_EL->
n;
771 for (
G4int i=0 ; i < panEPM_T_EL->
n ; i++ )
774 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 ] ) );
775 anEPM_T_E.
isoAngle.push_back( angle );
780 G4cout <<
"G4NeutronHPThermalScattering Do not Suuport yet." <<
G4endl;
814 for (
G4int i = 0 ; i < n-1 ; i++ )
819 sum_p += ( ( anE_P_E_isoAng->
prob[i] ) * dE );
821 if ( random <= sum_p/total )
823 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 ) );
824 secondary_energy = secondary_energy*
eV;
830 return secondary_energy;
838 std::map< G4double , G4int > map_energy;
840 std::vector< G4double > v_energy;
842 std::vector< E_P_E_isoAng* >::iterator itv;
844 for ( itv = vNEP_EPM->begin(); itv != vNEP_EPM->end(); itv++ )
846 v_energy.push_back( (*itv)->energy );
847 map_energy.insert( std::pair < G4double , G4int > ( (*itv)->energy , i ) );
851 std::pair < G4double , G4double > energyLH =
find_LH ( pE , &v_energy );
856 if ( energyLH.first != 0.0 && energyLH.second != 0.0 )
858 pE_P_E_isoAng_EL = (*vNEP_EPM)[ map_energy.find ( energyLH.first )->second ];
859 pE_P_E_isoAng_EH = (*vNEP_EPM)[ map_energy.find ( energyLH.second )->second ];
861 else if ( energyLH.first == 0.0 )
863 pE_P_E_isoAng_EL = (*vNEP_EPM)[ 0 ];
864 pE_P_E_isoAng_EH = (*vNEP_EPM)[ 1 ];
866 if ( energyLH.second == 0.0 )
868 pE_P_E_isoAng_EH = (*vNEP_EPM).back();
869 itv = vNEP_EPM->end();
872 pE_P_E_isoAng_EL = *itv;
884 sE =
get_linear_interpolated ( pE , std::pair < G4double , G4double > ( energyLH.first , sE_L ) , std::pair < G4double , G4double > ( energyLH.second , sE_H ) );
894 if ( E_isoAng_L.
n == E_isoAng_H.
n )
896 anE_isoAng.
n = E_isoAng_L.
n;
897 for (
G4int j=0 ; j < anE_isoAng.
n ; j++ )
901 anE_isoAng.
isoAngle.push_back( angle );
912 return std::pair< G4double , E_isoAng >( sE , anE_isoAng);
919 std::map < G4String , G4int > co_dic;
924 for (
size_t i = 0 ; i < numberOfMaterials ; i++ )
926 G4Material* material = (*theMaterialTable)[i];
928 for (
size_t j = 0 ; j < numberOfElements ; j++ )
933 G4int ts_ID_of_this_geometry;
935 if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
937 ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) ->
second;
941 ts_ID_of_this_geometry = co_dic.size();
942 co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
949 dic.insert( std::pair < std::pair < G4Material* , const G4Element* > ,
G4int > ( std::pair < G4Material* , const G4Element* > ( material , element ) , ts_ID_of_this_geometry ) );
958 for (
size_t i = 0 ; i < numberOfElements ; i++ )
960 const G4Element* element = (*theElementTable)[i];
965 G4int ts_ID_of_this_geometry;
967 if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
969 ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) ->
second;
973 ts_ID_of_this_geometry = co_dic.size();
974 co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
981 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 ) );
987 G4cout <<
"Neutron HP Thermal Scattering: Following material-element pairs or elements are registered." <<
G4endl;
988 for ( std::map < std::pair < const G4Material* , const G4Element* > ,
G4int >::iterator it =
dic.begin() ; it !=
dic.end() ; it++ )
990 if ( it->first.first != NULL )
992 G4cout <<
"Material " << it->first.first->GetName() <<
" - Element " << it->first.second->GetName() <<
", internal thermal scattering id " << it->second <<
G4endl;
996 G4cout <<
"Element " << it->first.second->GetName() <<
", internal thermal scattering id " << it->second <<
G4endl;
1004 if ( !getenv(
"G4NEUTRONHPDATA" ) )
1005 throw G4HadronicException(__FILE__, __LINE__,
"Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
1006 dirName = getenv(
"G4NEUTRONHPDATA" );
1012 for ( std::map < G4String , G4int >::iterator it = co_dic.begin() ; it != co_dic.end() ; it++ )
1015 G4int ts_ID = it->second;
1018 G4String fsName =
"/ThermalScattering/Coherent/FS/";
1019 G4String fileName = dirName + fsName + tsndlName;
1023 fsName =
"/ThermalScattering/Incoherent/FS/";
1024 fileName = dirName + fsName + tsndlName;
1028 fsName =
"/ThermalScattering/Inelastic/FS/";
1029 fileName = dirName + fsName + tsndlName;
1040 if (
dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) ) !=
dic.end() )
1041 result =
dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) )->second;
G4int getTS_ID(const G4Material *, const G4Element *)
E_isoAng * readAnE_isoAng(std::istream *)
std::map< G4double, std::vector< E_isoAng * > * > * readAnIncoherentFSDATA(G4String)
std::vector< G4double > isoAngle
static G4NeutronHPManager * GetInstance()
G4double GetCoherentCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
const G4String & GetName() const
E_isoAng create_E_isoAng_from_energy(G4double, std::vector< E_isoAng * > *)
static G4MaterialTable * GetMaterialTable()
std::vector< G4Material * > G4MaterialTable
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus)
std::pair< G4double, E_isoAng > create_sE_and_EPM_from_pE_and_vE_P_E_isoAng(G4double, G4double, std::vector< E_P_E_isoAng * > *)
void GetDataStream(G4String, std::istringstream &iss)
G4double getMu(E_isoAng *)
const G4Element * GetElement(G4int iel) const
std::vector< E_isoAng * > vE_isoAngle
std::pair< G4double, G4double > find_LH(G4double, std::vector< G4double > *)
std::map< std::pair< const G4Material *, const G4Element * >, G4int > dic
G4NeutronHPElastic * theHPElastic
void SetMinEnergy(G4double anEnergy)
std::map< G4double, std::vector< E_P_E_isoAng * > * > * readAnInelasticFSDATA(G4String)
void AddThermalElement(G4String, G4String)
G4GLOB_DLL std::ostream G4cout
static size_t GetNumberOfElements()
const G4ParticleDefinition * GetDefinition() const
G4double get_linear_interpolated(G4double, std::pair< G4double, G4double >, std::pair< G4double, G4double >)
G4double GetKineticEnergy() const
std::map< G4int, std::map< G4double, std::vector< E_P_E_isoAng * > * > * > inelasticFSs
static const double second
static const double perCent
static G4Neutron * Neutron()
static size_t GetNumberOfMaterials()
const G4LorentzVector & Get4Momentum() const
~G4NeutronHPThermalScattering()
void AddUserThermalScatteringFile(G4String, G4String)
static const double kelvin
void BuildPhysicsTable(const G4ParticleDefinition &)
G4double GetCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
std::map< G4int, std::map< G4double, std::vector< E_isoAng * > * > * > incoherentFSs
E_P_E_isoAng * readAnE_P_E_isoAng(std::istream *)
void SetEnergyChange(G4double anEnergy)
G4double total(Particle const *const p1, Particle const *const p2)
G4double energy(const ThreeVector &p, const G4double m)
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels() const
std::map< G4double, std::vector< std::pair< G4double, G4double > * > * > * readACoherentFSDATA(G4String)
G4double GetInelasticCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
G4String GetTS_NDL_Name(G4String nameG4Element)
void SetMaxEnergy(const G4double anEnergy)
std::map< G4int, std::map< G4double, std::vector< std::pair< G4double, G4double > * > * > * > coherentFSs
G4HadFinalState theParticleChange
G4NeutronHPThermalScatteringNames names
G4double GetTemperature() const
const G4Material * GetMaterial() const
size_t GetNumberOfElements() const
std::vector< G4Element * > G4ElementTable
const G4String & GetName() const
G4NeutronHPThermalScattering()
static G4ElementTable * GetElementTable()
void AddUserThermalScatteringFile(G4String, G4String)
size_t sizeOfMaterialTable
void SetMomentumChange(const G4ThreeVector &aV)
G4bool IsThisThermalElement(G4String)
G4NeutronHPThermalScatteringData * theXSection
G4double get_secondary_energy_from_E_P_E_isoAng(G4double, E_P_E_isoAng *)
std::vector< G4double > prob