// // ******************************************************************** // * License and Disclaimer * // * * // * The Geant4 software is copyright of the Copyright Holders of * // * the Geant4 Collaboration. It is provided under the terms and * // * conditions of the Geant4 Software License, included in the file * // * LICENSE and available at http://cern.ch/geant4/license . These * // * include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GEANT4 collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the Geant4 Software license. * // ******************************************************************** // // Thermal Neutron Scattering // Koi, Tatsumi (SLAC/SCCS) // // Class Description: // // Final State Generators for a high precision (based on evaluated data // libraries) description of themal neutron scattering below 4 eV; // Based on Thermal neutron scattering files // from the evaluated nuclear data files ENDF/B-VI, Release2 // To be used in your physics list in case you need this physics. // In this case you want to register an object of this class with // the corresponding process. // 070625 Fix memory leaking at destructor by T. Koi // 081201 Fix memory leaking at destructor by T. Koi // 100729 Add model name in constructor Problem #1116 #include "G4NeutronHPThermalScattering.hh" #include "G4Neutron.hh" #include "G4ElementTable.hh" G4NeutronHPThermalScattering::G4NeutronHPThermalScattering() :G4HadronicInteraction("NeutronHPThermalScattering") { theHPElastic = new G4NeutronHPElastic(); SetMinEnergy( 0.*eV ); SetMaxEnergy( 4*eV ); theXSection = new G4NeutronHPThermalScatteringData(); theXSection->BuildPhysicsTable( *(G4Neutron::Neutron()) ); // Check Elements std::vector< G4int > indexOfThermalElement; static const G4ElementTable* theElementTable = G4Element::GetElementTable(); size_t numberOfElements = G4Element::GetNumberOfElements(); for ( size_t i = 0 ; i < numberOfElements ; i++ ) { if ( names.IsThisThermalElement ( (*theElementTable)[i]->GetName() ) ) { indexOfThermalElement.push_back( i ); } } if ( !getenv("G4NEUTRONHPDATA") ) throw G4HadronicException(__FILE__, __LINE__, "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section files."); dirName = getenv("G4NEUTRONHPDATA"); // Read data // Element (id) -> FS Type -> read file for ( size_t i = 0 ; i < indexOfThermalElement.size() ; i++ ) { //G4cout << "G4NeutronHPThermalScattering " << (*theElementTable)[i]->GetName() << G4endl; G4String tsndlName = names.GetTS_NDL_Name ( (*theElementTable)[ indexOfThermalElement[ i ] ]->GetName() ); //G4cout << "G4NeutronHPThermalScattering " << tsndlName << std::endl; // coherent elastic G4String fsName = "/ThermalScattering/Coherent/FS/"; G4String fileName = dirName + fsName + tsndlName; coherentFSs.insert ( std::pair < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* > ( indexOfThermalElement[ i ] , readACoherentFSDATA( fileName ) ) ); // incoherent elastic fsName = "/ThermalScattering/Incoherent/FS/"; fileName = dirName + fsName + tsndlName; incoherentFSs.insert ( std::pair < G4int , std::map < G4double , std::vector < E_isoAng* >* >* > ( indexOfThermalElement[ i ] , readAnIncoherentFSDATA( fileName ) ) ); // inelastic fsName = "/ThermalScattering/Inelastic/FS/"; fileName = dirName + fsName + tsndlName; inelasticFSs.insert ( std::pair < G4int , std::map < G4double , std::vector < E_P_E_isoAng* >* >* > ( indexOfThermalElement[ i ] , readAnInelasticFSDATA( fileName ) ) ); } } G4NeutronHPThermalScattering::~G4NeutronHPThermalScattering() { for ( std::map < G4int , std::map < G4double , std::vector < E_isoAng* >* >* >::iterator it = incoherentFSs.begin() ; it != incoherentFSs.end() ; it++ ) { std::map < G4double , std::vector < E_isoAng* >* >::iterator itt; for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ ) { std::vector< E_isoAng* >::iterator ittt; for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ ) { delete *ittt; } delete itt->second; } delete it->second; } for ( std::map < G4int , std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* >::iterator it = coherentFSs.begin() ; it != coherentFSs.end() ; it++ ) { std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >::iterator itt; for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ ) { std::vector < std::pair< G4double , G4double >* >::iterator ittt; for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ ) { delete *ittt; } delete itt->second; } delete it->second; } for ( std::map < G4int , std::map < G4double , std::vector < E_P_E_isoAng* >* >* >::iterator it = inelasticFSs.begin() ; it != inelasticFSs.end() ; it++ ) { std::map < G4double , std::vector < E_P_E_isoAng* >* >::iterator itt; for ( itt = it->second->begin() ; itt != it->second->end() ; itt++ ) { std::vector < E_P_E_isoAng* >::iterator ittt; for ( ittt = itt->second->begin(); ittt != itt->second->end() ; ittt++ ) { std::vector < E_isoAng* >::iterator it4; for ( it4 = (*ittt)->vE_isoAngle.begin() ; it4 != (*ittt)->vE_isoAngle.end() ; it4++ ) { delete *it4; } delete *ittt; } delete itt->second; } delete it->second; } delete theHPElastic; delete theXSection; } std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* G4NeutronHPThermalScattering::readACoherentFSDATA( G4String name ) { std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >* aCoherentFSDATA = new std::map < G4double , std::vector < std::pair< G4double , G4double >* >* >; std::ifstream theChannel( name.c_str() ); std::vector< G4double > vBraggE; G4int dummy; while ( theChannel >> dummy ) // MF { theChannel >> dummy; // MT G4double temp; theChannel >> temp; std::vector < std::pair< G4double , G4double >* >* anBragE_P = new std::vector < std::pair< G4double , G4double >* >; G4int n; theChannel >> n; for ( G4int i = 0 ; i < n ; i++ ) { G4double Ei; G4double Pi; if ( aCoherentFSDATA->size() == 0 ) { theChannel >> Ei; vBraggE.push_back( Ei ); } else { Ei = vBraggE[ i ]; } theChannel >> Pi; anBragE_P->push_back ( new std::pair < G4double , G4double > ( Ei , Pi ) ); //G4cout << "Coherent Elastic " << Ei << " " << Pi << G4endl; } aCoherentFSDATA->insert ( std::pair < G4double , std::vector < std::pair< G4double , G4double >* >* > ( temp , anBragE_P ) ); } return aCoherentFSDATA; } std::map < G4double , std::vector < E_P_E_isoAng* >* >* G4NeutronHPThermalScattering::readAnInelasticFSDATA ( G4String name ) { std::map < G4double , std::vector < E_P_E_isoAng* >* >* anT_E_P_E_isoAng = new std::map < G4double , std::vector < E_P_E_isoAng* >* >; std::ifstream theChannel( name.c_str() ); G4int dummy; while ( theChannel >> dummy ) // MF { theChannel >> dummy; // MT G4double temp; theChannel >> temp; std::vector < E_P_E_isoAng* >* vE_P_E_isoAng = new std::vector < E_P_E_isoAng* >; G4int n; theChannel >> n; for ( G4int i = 0 ; i < n ; i++ ) { vE_P_E_isoAng->push_back ( readAnE_P_E_isoAng ( &theChannel ) ); } anT_E_P_E_isoAng->insert ( std::pair < G4double , std::vector < E_P_E_isoAng* >* > ( temp , vE_P_E_isoAng ) ); } theChannel.close(); return anT_E_P_E_isoAng; } E_P_E_isoAng* G4NeutronHPThermalScattering::readAnE_P_E_isoAng( std::ifstream* file ) { E_P_E_isoAng* aData = new E_P_E_isoAng; G4double dummy; G4double energy; G4int nep , nl; *file >> dummy; *file >> energy; aData->energy = energy*eV; *file >> dummy; *file >> dummy; *file >> nep; *file >> nl; aData->n = nep/nl; for ( G4int i = 0 ; i < aData->n ; i++ ) { G4double prob; E_isoAng* anE_isoAng = new E_isoAng; aData->vE_isoAngle.push_back( anE_isoAng ); *file >> energy; anE_isoAng->energy = energy*eV; anE_isoAng->n = nl - 2; anE_isoAng->isoAngle.resize( anE_isoAng->n ); *file >> prob; aData->prob.push_back( prob ); //G4cout << "G4NeutronHPThermalScattering inelastic " << energy/eV << " " << i << " " << prob << " " << aData->prob[ i ] << G4endl; for ( G4int j = 0 ; j < anE_isoAng->n ; j++ ) { G4double x; *file >> x; anE_isoAng->isoAngle[j] = x ; //G4cout << "G4NeutronHPThermalScattering inelastic " << x << anE_isoAng->isoAngle[j] << G4endl; } } // Calcuate sum_of_provXdEs G4double total = 0; for ( G4int i = 0 ; i < aData->n - 1 ; i++ ) { G4double E_L = aData->vE_isoAngle[i]->energy/eV; G4double E_H = aData->vE_isoAngle[i+1]->energy/eV; G4double dE = E_H - E_L; total += ( ( aData->prob[i] ) * dE ); } aData->sum_of_probXdEs = total; return aData; } std::map < G4double , std::vector < E_isoAng* >* >* G4NeutronHPThermalScattering::readAnIncoherentFSDATA ( G4String name ) { std::map < G4double , std::vector < E_isoAng* >* >* T_E = new std::map < G4double , std::vector < E_isoAng* >* >; std::ifstream theChannel( name.c_str() ); G4int dummy; while ( theChannel >> dummy ) // MF { theChannel >> dummy; // MT G4double temp; theChannel >> temp; std::vector < E_isoAng* >* vE_isoAng = new std::vector < E_isoAng* >; G4int n; theChannel >> n; for ( G4int i = 0 ; i < n ; i++ ) vE_isoAng->push_back ( readAnE_isoAng( &theChannel ) ); T_E->insert ( std::pair < G4double , std::vector < E_isoAng* >* > ( temp , vE_isoAng ) ); } theChannel.close(); return T_E; } E_isoAng* G4NeutronHPThermalScattering::readAnE_isoAng( std::ifstream* file ) { E_isoAng* aData = new E_isoAng; G4double dummy; G4double energy; G4int n; *file >> dummy; *file >> energy; *file >> dummy; *file >> dummy; *file >> n; *file >> dummy; aData->energy = energy*eV; aData->n = n-2; aData->isoAngle.resize( n ); *file >> dummy; *file >> dummy; for ( G4int i = 0 ; i < aData->n ; i++ ) *file >> aData->isoAngle[i]; return aData; } G4HadFinalState* G4NeutronHPThermalScattering::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aNucleus ) { // Select Element > Reaction > const G4Material * theMaterial = aTrack.GetMaterial(); G4double aTemp = theMaterial->GetTemperature(); G4int n = theMaterial->GetNumberOfElements(); static const G4ElementTable* theElementTable = G4Element::GetElementTable(); G4bool findThermalElement = false; G4int ielement; for ( G4int i = 0; i < n ; i++ ) { G4int index = theMaterial->GetElement(i)->GetIndex(); if ( aNucleus.GetZ() == (*theElementTable)[index]->GetZ() && ( names.IsThisThermalElement ( (*theElementTable)[index]->GetName() ) ) ) { ielement = index; findThermalElement = true; break; } } if ( findThermalElement == true ) { // Select Reaction (Inelastic, coherent, incoherent) G4ParticleDefinition* pd = const_cast< G4ParticleDefinition* >( aTrack.GetDefinition() ); G4DynamicParticle* dp = new G4DynamicParticle ( pd , aTrack.Get4Momentum() ); G4double total = theXSection->GetCrossSection( dp , (*theElementTable)[ ielement ] , aTemp ); G4double inelastic = theXSection->GetInelasticCrossSection( dp , (*theElementTable)[ ielement ] , aTemp ); G4double random = G4UniformRand(); if ( random <= inelastic/total ) { // Inelastic // T_L and T_H std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator it; std::vector v_temp; v_temp.clear(); for ( it = inelasticFSs.find( ielement )->second->begin() ; it != inelasticFSs.find( ielement )->second->end() ; it++ ) { v_temp.push_back( it->first ); } // T_L T_H std::pair < G4double , G4double > tempLH = find_LH ( aTemp , &v_temp ); // // For T_L aNEP_EPM_TL and T_H aNEP_EPM_TH // std::vector< E_P_E_isoAng* >* vNEP_EPM_TL = 0; std::vector< E_P_E_isoAng* >* vNEP_EPM_TH = 0; if ( tempLH.first != 0.0 && tempLH.second != 0.0 ) { vNEP_EPM_TL = inelasticFSs.find( ielement )->second->find ( tempLH.first/kelvin )->second; vNEP_EPM_TH = inelasticFSs.find( ielement )->second->find ( tempLH.second/kelvin )->second; } else if ( tempLH.first == 0.0 ) { std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator itm; itm = inelasticFSs.find( ielement )->second->begin(); vNEP_EPM_TL = itm->second; itm++; vNEP_EPM_TH = itm->second; } else if ( tempLH.second == 0.0 ) { std::map < G4double , std::vector< E_P_E_isoAng* >* >::iterator itm; itm = inelasticFSs.find( ielement )->second->end(); itm--; vNEP_EPM_TH = itm->second; itm--; vNEP_EPM_TL = itm->second; } // G4double rand_for_sE = G4UniformRand(); 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 ); 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 ); G4double sE; sE = get_linear_interpolated ( aTemp , std::pair < G4double , G4double > ( tempLH.first , TL.first ) , std::pair < G4double , G4double > ( tempLH.second , TH.first ) ); E_isoAng anE_isoAng; if ( TL.second.n == TH.second.n ) { anE_isoAng.energy = sE; anE_isoAng.n = TL.second.n; for ( G4int i=0 ; i < anE_isoAng.n ; i++ ) { G4double angle; 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 ] ) ); anE_isoAng.isoAngle.push_back( angle ); } } else { std::cout << "Do not Suuport yet." << std::endl; } //set theParticleChange.SetEnergyChange( sE ); G4double mu = getMu( &anE_isoAng ); theParticleChange.SetMomentumChange( 0.0 , std::sqrt ( 1 - mu*mu ) , mu ); } else if ( random <= ( inelastic + theXSection->GetCoherentCrossSection( dp , (*theElementTable)[ ielement ] , aTemp ) ) / total ) { // Coherent Elastic G4double E = aTrack.GetKineticEnergy(); // T_L and T_H std::map < G4double , std::vector< std::pair< G4double , G4double >* >* >::iterator it; std::vector v_temp; v_temp.clear(); for ( it = coherentFSs.find( ielement )->second->begin() ; it != coherentFSs.find( ielement )->second->end() ; it++ ) { v_temp.push_back( it->first ); } // T_L T_H std::pair < G4double , G4double > tempLH = find_LH ( aTemp , &v_temp ); // // // For T_L anEPM_TL and T_H anEPM_TH // std::vector< std::pair< G4double , G4double >* >* pvE_p_TL = 0; std::vector< std::pair< G4double , G4double >* >* pvE_p_TH = 0; if ( tempLH.first != 0.0 && tempLH.second != 0.0 ) { pvE_p_TL = coherentFSs.find( ielement )->second->find ( tempLH.first/kelvin )->second; pvE_p_TH = coherentFSs.find( ielement )->second->find ( tempLH.first/kelvin )->second; } else if ( tempLH.first == 0.0 ) { pvE_p_TL = coherentFSs.find( ielement )->second->find ( v_temp[ 0 ] )->second; pvE_p_TH = coherentFSs.find( ielement )->second->find ( v_temp[ 1 ] )->second; } else if ( tempLH.second == 0.0 ) { pvE_p_TL = coherentFSs.find( ielement )->second->find ( v_temp.back() )->second; std::vector< G4double >::iterator itv; itv = v_temp.end(); itv--; itv--; pvE_p_TL = coherentFSs.find( ielement )->second->find ( *itv )->second; } std::vector< G4double > vE_T; std::vector< G4double > vp_T; G4int n1 = pvE_p_TL->size(); //G4int n2 = pvE_p_TH->size(); for ( G4int i=1 ; i < n1 ; i++ ) { if ( (*pvE_p_TL)[i]->first != (*pvE_p_TH)[i]->first ) abort(); vE_T.push_back ( (*pvE_p_TL)[i]->first ); 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 ) ) ); } G4int j = 0; for ( G4int i = 1 ; i < n ; i++ ) { if ( E/eV < vE_T[ i ] ) { j = i-1; break; } } G4double rand_for_mu = G4UniformRand(); G4int k = 0; for ( G4int i = 1 ; i < j ; i++ ) { G4double Pi = vp_T[ i ] / vp_T[ j ]; if ( rand_for_mu < Pi ) { k = i-1; break; } } G4double Ei = vE_T[ j ]; G4double mu = 1 - 2 * Ei / (E/eV) ; theParticleChange.SetEnergyChange( E ); theParticleChange.SetMomentumChange( 0.0 , std::sqrt ( 1 - mu*mu ) , mu ); } else { // InCoherent Elastic // T_L and T_H std::map < G4double , std::vector < E_isoAng* >* >::iterator it; std::vector v_temp; v_temp.clear(); for ( it = incoherentFSs.find( ielement )->second->begin() ; it != incoherentFSs.find( ielement )->second->end() ; it++ ) { v_temp.push_back( it->first ); } // T_L T_H std::pair < G4double , G4double > tempLH = find_LH ( aTemp , &v_temp ); // // For T_L anEPM_TL and T_H anEPM_TH // E_isoAng anEPM_TL_E; E_isoAng anEPM_TH_E; if ( tempLH.first != 0.0 && tempLH.second != 0.0 ) { anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( tempLH.first/kelvin )->second ); anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( tempLH.second/kelvin )->second ); } else if ( tempLH.first == 0.0 ) { anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( v_temp[ 0 ] )->second ); anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( v_temp[ 1 ] )->second ); } else if ( tempLH.second == 0.0 ) { anEPM_TH_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( v_temp.back() )->second ); std::vector< G4double >::iterator itv; itv = v_temp.end(); itv--; itv--; anEPM_TL_E = create_E_isoAng_from_energy ( aTrack.GetKineticEnergy() , incoherentFSs.find( ielement )->second->find ( *itv )->second ); } // E_isoAng for aTemp and aTrack.GetKineticEnergy() E_isoAng anEPM_T_E; if ( anEPM_TL_E.n == anEPM_TH_E.n ) { anEPM_T_E.n = anEPM_TL_E.n; for ( G4int i=0 ; i < anEPM_TL_E.n ; i++ ) { G4double angle; 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 ] ) ); anEPM_T_E.isoAngle.push_back( angle ); } } else { std::cout << "Do not Suuport yet." << std::endl; } // Decide mu G4double mu = getMu ( &anEPM_T_E ); // Set Final State theParticleChange.SetEnergyChange( aTrack.GetKineticEnergy() ); // No energy change in Elastic theParticleChange.SetMomentumChange( 0.0 , std::sqrt ( 1 - mu*mu ) , mu ); } delete dp; return &theParticleChange; } else { // Not thermal element // Neutron HP will handle return theHPElastic -> ApplyYourself( aTrack, aNucleus ); } } G4double G4NeutronHPThermalScattering::getMu( E_isoAng* anEPM ) { G4double random = G4UniformRand(); G4double result = 0.0; G4int in = int ( random * ( (*anEPM).n ) ); if ( in != 0 ) { G4double mu_l = (*anEPM).isoAngle[ in-1 ]; G4double mu_h = (*anEPM).isoAngle[ in ]; result = ( mu_h - mu_l ) * ( random * ( (*anEPM).n ) - in ) + mu_l; } else { G4double x = random * (*anEPM).n; G4double D = ( (*anEPM).isoAngle[ 0 ] - ( -1 ) ) + ( 1 - (*anEPM).isoAngle[ (*anEPM).n - 1 ] ); G4double ratio = ( (*anEPM).isoAngle[ 0 ] - ( -1 ) ) / D; if ( x <= ratio ) { G4double mu_l = -1; G4double mu_h = (*anEPM).isoAngle[ 0 ]; result = ( mu_h - mu_l ) * x + mu_l; } else { G4double mu_l = (*anEPM).isoAngle[ (*anEPM).n - 1 ]; G4double mu_h = 1; result = ( mu_h - mu_l ) * x + mu_l; } } return result; } std::pair < G4double , G4double > G4NeutronHPThermalScattering::find_LH ( G4double x , std::vector< G4double >* aVector ) { G4double L = 0.0; G4double H = 0.0; std::vector< G4double >::iterator it; for ( it = aVector->begin() ; it != aVector->end() ; it++ ) { if ( x <= *it ) { H = *it; if ( it != aVector->begin() ) { it--; L = *it; } else { L = 0.0; } break; } } if ( H == 0.0 ) L = aVector->back(); return std::pair < G4double , G4double > ( L , H ); } G4double G4NeutronHPThermalScattering::get_linear_interpolated ( G4double x , std::pair< G4double , G4double > Low , std::pair< G4double , G4double > High ) { G4double y=0.0; if ( High.first - Low.first != 0 ) y = ( High.second - Low.second ) / ( High.first - Low.first ) * ( x - Low.first ) + Low.second; else std::cout << "G4NeutronHPThermalScattering liner interpolation err!!" << std::endl; return y; } E_isoAng G4NeutronHPThermalScattering::create_E_isoAng_from_energy ( G4double energy , std::vector< E_isoAng* >* vEPM ) { E_isoAng anEPM_T_E; std::vector< E_isoAng* >::iterator iv; std::vector< G4double > v_e; v_e.clear(); for ( iv = vEPM->begin() ; iv != vEPM->end() ; iv++ ) v_e.push_back ( (*iv)->energy ); std::pair < G4double , G4double > energyLH = find_LH ( energy , &v_e ); //std::cout << " " << energy/eV << " " << energyLH.first/eV << " " << energyLH.second/eV << std::endl; E_isoAng* panEPM_T_EL=0; E_isoAng* panEPM_T_EH=0; if ( energyLH.first != 0.0 && energyLH.second != 0.0 ) { for ( iv = vEPM->begin() ; iv != vEPM->end() ; iv++ ) { if ( energyLH.first == (*iv)->energy ) break; } panEPM_T_EL = *iv; iv++; panEPM_T_EH = *iv; } else if ( energyLH.first == 0.0 ) { panEPM_T_EL = (*vEPM)[0]; panEPM_T_EH = (*vEPM)[1]; } else if ( energyLH.second == 0.0 ) { panEPM_T_EH = (*vEPM).back(); iv = vEPM->end(); iv--; iv--; panEPM_T_EL = *iv; } if ( panEPM_T_EL->n == panEPM_T_EH->n ) { anEPM_T_E.energy = energy; anEPM_T_E.n = panEPM_T_EL->n; for ( G4int i=0 ; i < panEPM_T_EL->n ; i++ ) { G4double angle; 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 ] ) ); anEPM_T_E.isoAngle.push_back( angle ); } } else { G4cout << "G4NeutronHPThermalScattering Do not Suuport yet." << G4endl; } return anEPM_T_E; } G4double G4NeutronHPThermalScattering::get_secondary_energy_from_E_P_E_isoAng ( G4double random , E_P_E_isoAng* anE_P_E_isoAng ) { G4double secondary_energy = 0.0; G4int n = anE_P_E_isoAng->n; G4double sum_p = 0.0; // sum_p_H G4double sum_p_L = 0.0; G4double total=0.0; /* delete for speed up for ( G4int i = 0 ; i < n-1 ; i++ ) { G4double E_L = anE_P_E_isoAng->vE_isoAngle[i]->energy/eV; G4double E_H = anE_P_E_isoAng->vE_isoAngle[i+1]->energy/eV; G4double dE = E_H - E_L; total += ( ( anE_P_E_isoAng->prob[i] ) * dE ); } if ( std::abs( total - anE_P_E_isoAng->sum_of_probXdEs ) > 1.0e-14 ) std::cout << total - anE_P_E_isoAng->sum_of_probXdEs << std::endl; */ total = anE_P_E_isoAng->sum_of_probXdEs; for ( G4int i = 0 ; i < n-1 ; i++ ) { G4double E_L = anE_P_E_isoAng->vE_isoAngle[i]->energy/eV; G4double E_H = anE_P_E_isoAng->vE_isoAngle[i+1]->energy/eV; G4double dE = E_H - E_L; sum_p += ( ( anE_P_E_isoAng->prob[i] ) * dE ); if ( random <= sum_p/total ) { 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 ) ); secondary_energy = secondary_energy*eV; //need eV break; } sum_p_L = sum_p; } return secondary_energy; } std::pair< G4double , E_isoAng > G4NeutronHPThermalScattering::create_sE_and_EPM_from_pE_and_vE_P_E_isoAng ( G4double rand_for_sE , G4double pE , std::vector < E_P_E_isoAng* >* vNEP_EPM ) { std::map< G4double , G4int > map_energy; map_energy.clear(); std::vector< G4double > v_energy; v_energy.clear(); std::vector< E_P_E_isoAng* >::iterator itv; G4int i = 0; for ( itv = vNEP_EPM->begin(); itv != vNEP_EPM->end(); itv++ ) { v_energy.push_back( (*itv)->energy ); map_energy.insert( std::pair < G4double , G4int > ( (*itv)->energy , i ) ); i++; } std::pair < G4double , G4double > energyLH = find_LH ( pE , &v_energy ); E_P_E_isoAng* pE_P_E_isoAng_EL = 0; E_P_E_isoAng* pE_P_E_isoAng_EH = 0; if ( energyLH.first != 0.0 && energyLH.second != 0.0 ) { pE_P_E_isoAng_EL = (*vNEP_EPM)[ map_energy.find ( energyLH.first )->second ]; pE_P_E_isoAng_EH = (*vNEP_EPM)[ map_energy.find ( energyLH.second )->second ]; } else if ( energyLH.first == 0.0 ) { pE_P_E_isoAng_EL = (*vNEP_EPM)[ 0 ]; pE_P_E_isoAng_EH = (*vNEP_EPM)[ 1 ]; } if ( energyLH.second == 0.0 ) { pE_P_E_isoAng_EH = (*vNEP_EPM).back(); itv = vNEP_EPM->end(); itv--; itv--; pE_P_E_isoAng_EL = *itv; } G4double sE; G4double sE_L; G4double sE_H; sE_L = get_secondary_energy_from_E_P_E_isoAng ( rand_for_sE , pE_P_E_isoAng_EL ); sE_H = get_secondary_energy_from_E_P_E_isoAng ( rand_for_sE , pE_P_E_isoAng_EH ); sE = get_linear_interpolated ( pE , std::pair < G4double , G4double > ( energyLH.first , sE_L ) , std::pair < G4double , G4double > ( energyLH.second , sE_H ) ); E_isoAng E_isoAng_L = create_E_isoAng_from_energy ( sE , &(pE_P_E_isoAng_EL->vE_isoAngle) ); E_isoAng E_isoAng_H = create_E_isoAng_from_energy ( sE , &(pE_P_E_isoAng_EH->vE_isoAngle) ); E_isoAng anE_isoAng; if ( E_isoAng_L.n == E_isoAng_H.n ) { anE_isoAng.n = E_isoAng_L.n; for ( G4int i=0 ; i < anE_isoAng.n ; i++ ) { G4double angle; angle = get_linear_interpolated ( sE , std::pair< G4double , G4double > ( sE_L , E_isoAng_L.isoAngle[ i ] ) , std::pair< G4double , G4double > ( sE_H , E_isoAng_H.isoAngle[ i ] ) ); anE_isoAng.isoAngle.push_back( angle ); } } else { std::cout << "Do not Suuport yet." << std::endl; } return std::pair< G4double , E_isoAng >( sE , anE_isoAng); }