75 std::map < G4double , std::vector < E_isoAng* >* >::iterator itt;
76 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
78 std::vector< E_isoAng* >::iterator ittt;
79 for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
88 for ( std::map <
G4int , std::map <
G4double , std::vector < std::pair< G4double , G4double >* >* >* >::iterator it =
coherentFSs.begin() ; it !=
coherentFSs.end() ; it++ )
90 std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >::iterator itt;
91 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
93 std::vector < std::pair< G4double , G4double >* >::iterator ittt;
94 for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
105 std::map < G4double , std::vector < E_P_E_isoAng* >* >::iterator itt;
106 for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ )
108 std::vector < E_P_E_isoAng* >::iterator ittt;
109 for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ )
111 std::vector < E_isoAng* >::iterator it4;
112 for ( it4 = (*ittt)->vE_isoAngle.begin() ; it4 != (*ittt)->vE_isoAngle.end() ; it4++ )
132 std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* aCoherentFSDATA =
new std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >;
135 std::istringstream theChannel(std::ios::in);
138 std::vector< G4double > vBraggE;
141 while ( theChannel >> dummy )
146 std::vector < std::pair< G4double , G4double >* >* anBragE_P =
new std::vector < std::pair< G4double , G4double >* >;
150 for (
G4int i = 0 ; i <
n ; i++ )
154 if ( aCoherentFSDATA->size() == 0 )
157 vBraggE.push_back( Ei );
164 anBragE_P->push_back (
new std::pair < G4double , G4double > ( Ei , Pi ) );
167 aCoherentFSDATA->insert ( std::pair <
G4double , std::vector < std::pair< G4double , G4double >* >* > ( temp , anBragE_P ) );
170 return aCoherentFSDATA;
177 std::map < G4double , std::vector < E_P_E_isoAng* >* >* anT_E_P_E_isoAng =
new std::map < G4double , std::vector < E_P_E_isoAng* >* >;
180 std::istringstream theChannel(std::ios::in);
184 while ( theChannel >> dummy )
189 std::vector < E_P_E_isoAng* >* vE_P_E_isoAng =
new std::vector < E_P_E_isoAng* >;
192 for (
G4int i = 0 ; i <
n ; i++ )
196 anT_E_P_E_isoAng->insert ( std::pair <
G4double , std::vector < E_P_E_isoAng* >* > ( temp , vE_P_E_isoAng ) );
200 return anT_E_P_E_isoAng;
220 for (
G4int i = 0 ; i < aData->
n ; i++ )
227 anE_isoAng->
n = nl - 2;
228 anE_isoAng->
isoAngle.resize( anE_isoAng->
n );
230 aData->
prob.push_back( prob );
232 for (
G4int j = 0 ; j < anE_isoAng->
n ; j++ )
243 for (
G4int i = 0 ; i < aData->
n - 1 ; i++ )
248 total += ( ( aData->
prob[i] ) * dE );
259 std::map < G4double , std::vector < E_isoAng* >* >* T_E =
new std::map < G4double , std::vector < E_isoAng* >* >;
262 std::istringstream theChannel(std::ios::in);
266 while ( theChannel >> dummy )
271 std::vector < E_isoAng* >* vE_isoAng =
new std::vector < E_isoAng* >;
274 for (
G4int i = 0 ; i <
n ; i++ )
276 T_E->insert ( std::pair <
G4double , std::vector < E_isoAng* >* > ( temp , vE_isoAng ) );
304 for (
G4int i = 0 ; i < aData->
n ; i++ )
328 G4bool findThermalElement =
false;
331 for (
G4int i = 0; i <
n ; i++ )
338 if (
getTS_ID( NULL , theElement ) != -1 )
340 ielement =
getTS_ID( NULL , theElement );
341 findThermalElement =
true;
344 else if (
getTS_ID( theMaterial , theElement ) != -1 )
346 ielement =
getTS_ID( theMaterial , theElement );
347 findThermalElement =
true;
353 if ( findThermalElement ==
true )
365 if ( random <= inelastic/total )
370 std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator it;
371 std::vector<G4double> v_temp;
373 for ( it =
inelasticFSs.find( ielement )->second->begin() ; it !=
inelasticFSs.find( ielement )->second->end() ; it++ )
375 v_temp.push_back( it->first );
379 std::pair < G4double , G4double > tempLH =
find_LH ( aTemp , &v_temp );
383 std::vector< E_P_E_isoAng* >* vNEP_EPM_TL = 0;
384 std::vector< E_P_E_isoAng* >* vNEP_EPM_TH = 0;
386 if ( tempLH.first != 0.0 && tempLH.second != 0.0 )
388 vNEP_EPM_TL =
inelasticFSs.find( ielement )->second->find ( tempLH.first/
kelvin )->second;
389 vNEP_EPM_TH =
inelasticFSs.find( ielement )->second->find ( tempLH.second/
kelvin )->second;
391 else if ( tempLH.first == 0.0 )
393 std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator itm;
395 vNEP_EPM_TL = itm->second;
397 vNEP_EPM_TH = itm->second;
398 tempLH.first = tempLH.second;
399 tempLH.second = itm->first;
401 else if ( tempLH.second == 0.0 )
403 std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator itm;
406 vNEP_EPM_TH = itm->second;
408 vNEP_EPM_TL = itm->second;
409 tempLH.second = tempLH.first;
410 tempLH.first = itm->first;
419 sE =
get_linear_interpolated ( aTemp , std::pair < G4double , G4double > ( tempLH.first , TL.first ) , std::pair < G4double , G4double > ( tempLH.second , TH.first ) );
423 if ( TL.second.n == TH.second.n )
426 anE_isoAng.
n = TL.second.n;
427 for (
G4int i=0 ; i < anE_isoAng.
n ; i++ )
430 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 ] ) );
431 anE_isoAng.
isoAngle.push_back( angle );
433 mu =
getMu( &anE_isoAng );
437 G4HadronicException(__FILE__, __LINE__,
"A problem is found in Thermal Scattering Data! Do not yet supported");
453 std::map < G4double , std::vector< std::pair< G4double , G4double >* >* >::iterator it;
454 std::vector<G4double> v_temp;
456 for ( it =
coherentFSs.find( ielement )->second->begin() ; it !=
coherentFSs.find( ielement )->second->end() ; it++ )
458 v_temp.push_back( it->first );
462 std::pair < G4double , G4double > tempLH =
find_LH ( aTemp , &v_temp );
467 std::vector< std::pair< G4double , G4double >* >* pvE_p_TL = NULL;
468 std::vector< std::pair< G4double , G4double >* >* pvE_p_TH = NULL;
470 if ( tempLH.first != 0.0 && tempLH.second != 0.0 )
472 pvE_p_TL =
coherentFSs.find( ielement )->second->find ( tempLH.first/
kelvin )->second;
473 pvE_p_TH =
coherentFSs.find( ielement )->second->find ( tempLH.first/
kelvin )->second;
475 else if ( tempLH.first == 0.0 )
477 pvE_p_TL =
coherentFSs.find( ielement )->second->find ( v_temp[ 0 ] )->second;
478 pvE_p_TH =
coherentFSs.find( ielement )->second->find ( v_temp[ 1 ] )->second;
479 tempLH.first = tempLH.second;
480 tempLH.second = v_temp[ 1 ];
482 else if ( tempLH.second == 0.0 )
485 std::vector< G4double >::iterator itv;
489 pvE_p_TL =
coherentFSs.find( ielement )->second->find ( *itv )->second;
490 tempLH.second = tempLH.first;
496 G4HadronicException(__FILE__, __LINE__,
"A problem is found in Thermal Scattering Data! Unexpected temperature values in data");
499 std::vector< G4double > vE_T;
500 std::vector< G4double > vp_T;
502 G4int n1 = pvE_p_TL->size();
505 for (
G4int i=1 ; i < n1 ; i++ )
507 if ( (*pvE_p_TL)[i]->first != (*pvE_p_TH)[i]->first )
G4HadronicException(__FILE__, __LINE__,
"A problem is found in Thermal Scattering Data!");
508 vE_T.push_back ( (*pvE_p_TL)[i]->first );
509 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 ) ) );
513 for (
G4int i = 1 ; i <
n ; i++ )
515 if ( E/
eV < vE_T[ i ] )
525 for (
G4int i = 1 ; i < j ; i++ )
527 G4double Pi = vp_T[ i ] / vp_T[ j ];
528 if ( rand_for_mu < Pi )
540 if ( mu < -1.0 ) mu = -1.0;
552 std::map < G4double , std::vector < E_isoAng* >* >::iterator it;
553 std::vector<G4double> v_temp;
557 v_temp.push_back( it->first );
561 std::pair < G4double , G4double > tempLH =
find_LH ( aTemp , &v_temp );
570 if ( tempLH.first != 0.0 && tempLH.second != 0.0 ) {
574 }
else if ( tempLH.first == 0.0 ) {
578 tempLH.first = tempLH.second;
579 tempLH.second = v_temp[ 1 ];
580 }
else if ( tempLH.second == 0.0 ) {
583 std::vector< G4double >::iterator itv;
588 tempLH.second = tempLH.first;
596 if ( anEPM_TL_E.
n == anEPM_TH_E.
n )
598 anEPM_T_E.
n = anEPM_TL_E.
n;
599 for (
G4int i=0 ; i < anEPM_TL_E.
n ; i++ )
602 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 ] ) );
603 anEPM_T_E.
isoAngle.push_back( angle );
605 mu =
getMu ( &anEPM_T_E );
609 G4HadronicException(__FILE__, __LINE__,
"A problem is found in Thermal Scattering Data! Do not yet supported");
639 G4int in = int ( random * ( (*anEPM).n ) );
643 G4double mu_l = (*anEPM).isoAngle[ in-1 ];
644 G4double mu_h = (*anEPM).isoAngle[ in ];
645 result = ( mu_h - mu_l ) * ( random * ( (*anEPM).n ) - in ) + mu_l;
650 G4double D = ( (*anEPM).isoAngle[ 0 ] - ( -1 ) ) + ( 1 - (*anEPM).isoAngle[ (*anEPM).n - 1 ] );
651 G4double ratio = ( (*anEPM).isoAngle[ 0 ] - ( -1 ) ) / D;
655 G4double mu_h = (*anEPM).isoAngle[ 0 ];
656 result = ( mu_h - mu_l ) * x + mu_l;
660 G4double mu_l = (*anEPM).isoAngle[ (*anEPM).n - 1 ];
662 result = ( mu_h - mu_l ) * x + mu_l;
676 if ( aVector->size() == 1 ) {
677 L = aVector->front();
678 H = aVector->front();
683 for ( std::vector< G4double >::iterator
684 it = aVector->begin() ; it != aVector->end() ; it++ ) {
687 if ( it != aVector->begin() ) {
699 if ( H == 0.0 ) L = aVector->back();
702 return std::pair < G4double , G4double > (
L , H );
710 if ( High.first - Low.first != 0 ) {
711 y = ( High.second - Low.second ) / ( High.first - Low.first ) * ( x - Low.first ) + Low.second;
713 if ( High.second == Low.second ) {
716 G4cout <<
"G4ParticleHPThermalScattering liner interpolation err!!" <<
G4endl;
729 std::vector< E_isoAng* >::iterator iv;
731 std::vector< G4double > v_e;
733 for ( iv = vEPM->begin() ; iv != vEPM->end() ; iv++ )
734 v_e.push_back ( (*iv)->energy );
736 std::pair < G4double , G4double > energyLH =
find_LH ( energy , &v_e );
742 if ( energyLH.first != 0.0 && energyLH.second != 0.0 )
744 for ( iv = vEPM->begin() ; iv != vEPM->end() ; iv++ )
746 if ( energyLH.first == (*iv)->energy )
753 else if ( energyLH.first == 0.0 )
755 panEPM_T_EL = (*vEPM)[0];
756 panEPM_T_EH = (*vEPM)[1];
758 else if ( energyLH.second == 0.0 )
760 panEPM_T_EH = (*vEPM).back();
767 if ( panEPM_T_EL->
n == panEPM_T_EH->
n )
770 anEPM_T_E.
n = panEPM_T_EL->
n;
772 for (
G4int i=0 ; i < panEPM_T_EL->
n ; i++ )
775 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 ] ) );
776 anEPM_T_E.
isoAngle.push_back( angle );
781 G4cout <<
"G4ParticleHPThermalScattering Do not Suuport yet." <<
G4endl;
815 for (
G4int i = 0 ; i < n-1 ; i++ )
820 sum_p += ( ( anE_P_E_isoAng->
prob[i] ) * dE );
822 if ( random <= sum_p/total )
824 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 ) );
825 secondary_energy = secondary_energy*
eV;
831 return secondary_energy;
839 std::map< G4double , G4int > map_energy;
841 std::vector< G4double > v_energy;
843 std::vector< E_P_E_isoAng* >::iterator itv;
845 for ( itv = vNEP_EPM->begin(); itv != vNEP_EPM->end(); itv++ )
847 v_energy.push_back( (*itv)->energy );
848 map_energy.insert( std::pair < G4double , G4int > ( (*itv)->energy , i ) );
852 std::pair < G4double , G4double > energyLH =
find_LH ( pE , &v_energy );
857 if ( energyLH.first != 0.0 && energyLH.second != 0.0 )
859 pE_P_E_isoAng_EL = (*vNEP_EPM)[ map_energy.find ( energyLH.first )->second ];
860 pE_P_E_isoAng_EH = (*vNEP_EPM)[ map_energy.find ( energyLH.second )->second ];
862 else if ( energyLH.first == 0.0 )
864 pE_P_E_isoAng_EL = (*vNEP_EPM)[ 0 ];
865 pE_P_E_isoAng_EH = (*vNEP_EPM)[ 1 ];
867 if ( energyLH.second == 0.0 )
869 pE_P_E_isoAng_EH = (*vNEP_EPM).back();
870 itv = vNEP_EPM->end();
873 pE_P_E_isoAng_EL = *itv;
885 sE =
get_linear_interpolated ( pE , std::pair < G4double , G4double > ( energyLH.first , sE_L ) , std::pair < G4double , G4double > ( energyLH.second , sE_H ) );
895 if ( E_isoAng_L.
n == E_isoAng_H.
n )
897 anE_isoAng.
n = E_isoAng_L.
n;
898 for (
G4int j=0 ; j < anE_isoAng.
n ; j++ )
902 anE_isoAng.
isoAngle.push_back( angle );
913 return std::pair< G4double , E_isoAng >( sE , anE_isoAng);
920 std::map < G4String , G4int > co_dic;
925 for (
size_t i = 0 ; i < numberOfMaterials ; i++ )
927 G4Material* material = (*theMaterialTable)[i];
929 for (
size_t j = 0 ; j < numberOfElements ; j++ )
934 G4int ts_ID_of_this_geometry;
936 if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
938 ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) ->
second;
942 ts_ID_of_this_geometry = co_dic.size();
943 co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
950 dic.insert( std::pair < std::pair < G4Material* , const G4Element* > ,
G4int > ( std::pair < G4Material* , const G4Element* > ( material , element ) , ts_ID_of_this_geometry ) );
959 for (
size_t i = 0 ; i < numberOfElements ; i++ )
961 const G4Element* element = (*theElementTable)[i];
966 G4int ts_ID_of_this_geometry;
968 if ( co_dic.find ( ts_ndl_name ) != co_dic.end() )
970 ts_ID_of_this_geometry = co_dic.find ( ts_ndl_name ) ->
second;
974 ts_ID_of_this_geometry = co_dic.size();
975 co_dic.insert ( std::pair< G4String , G4int >( ts_ndl_name , ts_ID_of_this_geometry ) );
982 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 ) );
988 G4cout <<
"Neutron HP Thermal Scattering: Following material-element pairs or elements are registered." <<
G4endl;
989 for ( std::map < std::pair < const G4Material* , const G4Element* > ,
G4int >::iterator it =
dic.begin() ; it !=
dic.end() ; it++ )
991 if ( it->first.first != NULL )
993 G4cout <<
"Material " << it->first.first->GetName() <<
" - Element " << it->first.second->GetName() <<
", internal thermal scattering id " << it->second <<
G4endl;
997 G4cout <<
"Element " << it->first.second->GetName() <<
", internal thermal scattering id " << it->second <<
G4endl;
1005 if ( !getenv(
"G4NEUTRONHPDATA" ) )
1006 throw G4HadronicException(__FILE__, __LINE__,
"Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files.");
1007 dirName = getenv(
"G4NEUTRONHPDATA" );
1013 for ( std::map < G4String , G4int >::iterator it = co_dic.begin() ; it != co_dic.end() ; it++ )
1016 G4int ts_ID = it->second;
1019 G4String fsName =
"/ThermalScattering/Coherent/FS/";
1020 G4String fileName = dirName + fsName + tsndlName;
1024 fsName =
"/ThermalScattering/Incoherent/FS/";
1025 fileName = dirName + fsName + tsndlName;
1029 fsName =
"/ThermalScattering/Inelastic/FS/";
1030 fileName = dirName + fsName + tsndlName;
1041 if (
dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) ) !=
dic.end() )
1042 result =
dic.find( std::pair < const G4Material* , const G4Element* > ( material , element ) )->second;
static G4ParticleHPManager * GetInstance()
G4ParticleHPThermalScatteringNames names
~G4ParticleHPThermalScattering()
std::vector< G4double > isoAngle
std::map< G4double, std::vector< std::pair< G4double, G4double > * > * > * readACoherentFSDATA(G4String)
G4double get_secondary_energy_from_E_P_E_isoAng(G4double, E_P_E_isoAng *)
size_t sizeOfMaterialTable
std::map< G4int, std::map< G4double, std::vector< E_isoAng * > * > * > incoherentFSs
const G4String & GetName() const
void AddThermalElement(G4String, G4String)
static G4MaterialTable * GetMaterialTable()
std::vector< G4Material * > G4MaterialTable
void GetDataStream(G4String, std::istringstream &iss)
E_isoAng create_E_isoAng_from_energy(G4double, std::vector< E_isoAng * > *)
std::pair< G4double, E_isoAng > create_sE_and_EPM_from_pE_and_vE_P_E_isoAng(G4double, G4double, std::vector< E_P_E_isoAng * > *)
const G4Element * GetElement(G4int iel) const
std::map< G4double, std::vector< E_P_E_isoAng * > * > * readAnInelasticFSDATA(G4String)
G4double GetCoherentCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
void AddUserThermalScatteringFile(G4String, G4String)
void SetMinEnergy(G4double anEnergy)
std::map< G4double, std::vector< E_isoAng * > * > * readAnIncoherentFSDATA(G4String)
G4GLOB_DLL std::ostream G4cout
static size_t GetNumberOfElements()
const G4ParticleDefinition * GetDefinition() const
std::map< G4int, std::map< G4double, std::vector< std::pair< G4double, G4double > * > * > * > coherentFSs
E_P_E_isoAng * readAnE_P_E_isoAng(std::istream *)
G4ParticleHPThermalScatteringData * theXSection
G4double GetKineticEnergy() const
G4int getTS_ID(const G4Material *, const G4Element *)
static const double second
static const double perCent
static G4Neutron * Neutron()
E_isoAng * readAnE_isoAng(std::istream *)
static size_t GetNumberOfMaterials()
const G4LorentzVector & Get4Momentum() const
static const double kelvin
G4double getMu(E_isoAng *)
void SetEnergyChange(G4double anEnergy)
G4double total(Particle const *const p1, Particle const *const p2)
G4double get_linear_interpolated(G4double, std::pair< G4double, G4double >, std::pair< G4double, G4double >)
void AddUserThermalScatteringFile(G4String, G4String)
G4double energy(const ThreeVector &p, const G4double m)
G4bool IsThisThermalElement(G4String)
G4ParticleHPElastic * theHPElastic
std::map< G4int, std::map< G4double, std::vector< E_P_E_isoAng * > * > * > inelasticFSs
std::pair< G4double, G4double > find_LH(G4double, std::vector< G4double > *)
void SetMaxEnergy(const G4double anEnergy)
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &aTargetNucleus)
G4HadFinalState theParticleChange
G4double GetTemperature() const
const G4Material * GetMaterial() const
size_t GetNumberOfElements() const
G4double GetCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
G4String GetTS_NDL_Name(G4String nameG4Element)
std::vector< G4Element * > G4ElementTable
G4double GetInelasticCrossSection(const G4DynamicParticle *, const G4Element *, const G4Material *)
const G4String & GetName() const
static G4ElementTable * GetElementTable()
std::map< std::pair< const G4Material *, const G4Element * >, G4int > dic
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels() const
G4ParticleHPThermalScattering()
void BuildPhysicsTable(const G4ParticleDefinition &)
void SetMomentumChange(const G4ThreeVector &aV)
std::vector< E_isoAng * > vE_isoAngle
std::vector< G4double > prob