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 "G4NeutronHPorLCaptureData.hh"
00035 #include "G4SystemOfUnits.hh"
00036 #include "G4Neutron.hh"
00037 #include "G4ElementTable.hh"
00038 #include "G4NeutronHPData.hh"
00039
00040 #include "G4PhysicsVector.hh"
00041
00042
00043 G4NeutronHPorLCaptureData::G4NeutronHPorLCaptureData()
00044 {
00045 SetMinKinEnergy( 0*MeV );
00046 SetMaxKinEnergy( 20*MeV );
00047
00048 ke_cache = 0.0;
00049 xs_cache = 0.0;
00050 element_cache = NULL;
00051 material_cache = NULL;
00052
00053 }
00054
00055 G4NeutronHPorLCaptureData::~G4NeutronHPorLCaptureData()
00056 {
00057
00058 }
00059
00060 G4bool G4NeutronHPorLCaptureData::IsIsoApplicable( const G4DynamicParticle* dp ,
00061 G4int , G4int ,
00062 const G4Element* element ,
00063 const G4Material* )
00064 {
00065 G4double eKin = dp->GetKineticEnergy();
00066 if ( eKin > GetMaxKinEnergy()
00067 || eKin < GetMinKinEnergy()
00068 || dp->GetDefinition() != G4Neutron::Neutron() ) return false;
00069 if ( unavailable_elements->find( element->GetName() ) != unavailable_elements->end() ) return false;
00070
00071 return true;
00072 }
00073
00074 G4double G4NeutronHPorLCaptureData::GetIsoCrossSection( const G4DynamicParticle* dp ,
00075 G4int , G4int ,
00076 const G4Isotope* ,
00077 const G4Element* element ,
00078 const G4Material* material )
00079 {
00080 if ( dp->GetKineticEnergy() == ke_cache && element == element_cache && material == material_cache ) return xs_cache;
00081
00082 ke_cache = dp->GetKineticEnergy();
00083 element_cache = element;
00084 material_cache = material;
00085 G4double xs = GetCrossSection( dp , element , material->GetTemperature() );
00086 xs_cache = xs;
00087 return xs;
00088
00089 }
00090
00091
00092 G4NeutronHPorLCaptureData::G4NeutronHPorLCaptureData( G4NeutronHPChannel* pChannel , std::set< G4String >* pSet )
00093 :G4VCrossSectionDataSet("NeutronHPorLCaptureXS")
00094 {
00095 theCaptureChannel = pChannel;
00096 unavailable_elements = pSet;
00097
00098 SetMinKinEnergy( 0*MeV );
00099 SetMaxKinEnergy( 20*MeV );
00100
00101 ke_cache = 0.0;
00102 xs_cache = 0.0;
00103 element_cache = NULL;
00104 material_cache = NULL;
00105 }
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118 void G4NeutronHPorLCaptureData::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 G4NeutronHPorLCaptureData::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 G4NeutronHPorLCaptureData::
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 = theCaptureChannel[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 }