00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 #include "G4NeutronHPorLElasticData.hh"
00035 #include "G4SystemOfUnits.hh"
00036 #include "G4Neutron.hh"
00037 #include "G4ElementTable.hh"
00038 #include "G4NeutronHPData.hh"
00039
00040 #include "G4PhysicsVector.hh"
00041
00042 G4NeutronHPorLElasticData::G4NeutronHPorLElasticData()
00043 {
00044 SetMinKinEnergy( 0*MeV );
00045 SetMaxKinEnergy( 20*MeV );
00046
00047 ke_cache = 0.0;
00048 xs_cache = 0.0;
00049 element_cache = NULL;
00050 material_cache = NULL;
00051
00052 }
00053
00054 G4NeutronHPorLElasticData::~G4NeutronHPorLElasticData()
00055 {
00056
00057 }
00058
00059 G4bool G4NeutronHPorLElasticData::IsIsoApplicable( const G4DynamicParticle* dp ,
00060 G4int , G4int ,
00061 const G4Element* element ,
00062 const G4Material* )
00063 {
00064 G4double eKin = dp->GetKineticEnergy();
00065 if ( eKin > GetMaxKinEnergy()
00066 || eKin < GetMinKinEnergy()
00067 || dp->GetDefinition() != G4Neutron::Neutron() ) return false;
00068 if ( unavailable_elements->find( element->GetName() ) != unavailable_elements->end() ) return false;
00069
00070 return true;
00071 }
00072
00073 G4double G4NeutronHPorLElasticData::GetIsoCrossSection( const G4DynamicParticle* dp ,
00074 G4int , G4int ,
00075 const G4Isotope* ,
00076 const G4Element* element ,
00077 const G4Material* material )
00078 {
00079 if ( dp->GetKineticEnergy() == ke_cache && element == element_cache && material == material_cache ) return xs_cache;
00080
00081 ke_cache = dp->GetKineticEnergy();
00082 element_cache = element;
00083 material_cache = material;
00084 G4double xs = GetCrossSection( dp , element , material->GetTemperature() );
00085 xs_cache = xs;
00086 return xs;
00087
00088 }
00089
00090
00091 G4NeutronHPorLElasticData::G4NeutronHPorLElasticData( G4NeutronHPChannel* pChannel , std::set< G4String >* pSet )
00092 :G4VCrossSectionDataSet("NeutronHPorLElasticXS")
00093 {
00094 theElasticChannel = pChannel;
00095 unavailable_elements = pSet;
00096
00097 SetMinKinEnergy( 0*MeV );
00098 SetMaxKinEnergy( 20*MeV );
00099
00100 ke_cache = 0.0;
00101 xs_cache = 0.0;
00102 element_cache = NULL;
00103 material_cache = NULL;
00104 }
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118 void G4NeutronHPorLElasticData::BuildPhysicsTable( const G4ParticleDefinition& aP )
00119 {
00120 if( &aP!=G4Neutron::Neutron() )
00121 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
00122 }
00123
00124
00125
00126 void G4NeutronHPorLElasticData::DumpPhysicsTable(const G4ParticleDefinition& aP)
00127 {
00128 if(&aP!=G4Neutron::Neutron())
00129 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
00130
00131 }
00132
00133
00134
00135 #include "G4Nucleus.hh"
00136 #include "G4NucleiProperties.hh"
00137 #include "G4Neutron.hh"
00138 #include "G4Electron.hh"
00139
00140 G4double G4NeutronHPorLElasticData::
00141 GetCrossSection(const G4DynamicParticle* aP, const G4Element*anE, G4double aT)
00142 {
00143
00144
00145 G4double result = 0;
00146
00147 G4int index = anE->GetIndex();
00148
00149
00150 G4double eKinetic = aP->GetKineticEnergy();
00151 G4ReactionProduct theNeutron( aP->GetDefinition() );
00152 theNeutron.SetMomentum( aP->GetMomentum() );
00153 theNeutron.SetKineticEnergy( eKinetic );
00154
00155
00156 G4Nucleus aNuc;
00157 G4double eps = 0.0001;
00158 G4double theA = anE->GetN();
00159 G4double theZ = anE->GetZ();
00160 G4double eleMass;
00161 eleMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theA+eps), static_cast<G4int>(theZ+eps))
00162 ) / G4Neutron::Neutron()->GetPDGMass();
00163
00164 G4ReactionProduct boosted;
00165 G4double aXsection;
00166
00167
00168 G4int counter = 0;
00169 G4double buffer = 0;
00170 G4int size = G4int(std::max(10., aT/60*kelvin));
00171 G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutron.GetMomentum();
00172 G4double neutronVMag = neutronVelocity.mag();
00173 while(counter == 0 || std::abs(buffer-result/std::max(1,counter)) > 0.03*buffer)
00174 {
00175 if(counter) buffer = result/counter;
00176 while (counter<size)
00177 {
00178 counter ++;
00179 G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus(eleMass, aT);
00180 boosted.Lorentz(theNeutron, aThermalNuc);
00181 G4double theEkin = boosted.GetKineticEnergy();
00182
00183 aXsection = theElasticChannel[index].GetXsec( theEkin );
00184
00185 G4ThreeVector targetVelocity = 1./aThermalNuc.GetMass()*aThermalNuc.GetMomentum();
00186 aXsection *= (targetVelocity-neutronVelocity).mag()/neutronVMag;
00187 result += aXsection;
00188 }
00189 size += size;
00190 }
00191 result /= counter;
00192 return result;
00193 }