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
00035
00036
00037 #include "G4NeutronHPElasticData.hh"
00038 #include "G4PhysicalConstants.hh"
00039 #include "G4SystemOfUnits.hh"
00040 #include "G4Neutron.hh"
00041 #include "G4ElementTable.hh"
00042 #include "G4NeutronHPData.hh"
00043
00044 G4NeutronHPElasticData::G4NeutronHPElasticData()
00045 :G4VCrossSectionDataSet("NeutronHPElasticXS")
00046 {
00047 SetMinKinEnergy( 0*MeV );
00048 SetMaxKinEnergy( 20*MeV );
00049
00050 ke_cache = 0.0;
00051 xs_cache = 0.0;
00052 element_cache = NULL;
00053 material_cache = NULL;
00054
00055 theCrossSections = 0;
00056 onFlightDB = true;
00057 BuildPhysicsTable( *G4Neutron::Neutron() );
00058 }
00059
00060 G4NeutronHPElasticData::~G4NeutronHPElasticData()
00061 {
00062 if ( theCrossSections != 0 ) theCrossSections->clearAndDestroy();
00063 delete theCrossSections;
00064 }
00065
00066 G4bool G4NeutronHPElasticData::IsIsoApplicable( const G4DynamicParticle* dp ,
00067 G4int , G4int ,
00068 const G4Element* ,
00069 const G4Material* )
00070 {
00071
00072 G4double eKin = dp->GetKineticEnergy();
00073 if ( eKin > GetMaxKinEnergy()
00074 || eKin < GetMinKinEnergy()
00075 || dp->GetDefinition() != G4Neutron::Neutron() ) return false;
00076
00077 return true;
00078 }
00079
00080 G4double G4NeutronHPElasticData::GetIsoCrossSection( const G4DynamicParticle* dp ,
00081 G4int , G4int ,
00082 const G4Isotope* ,
00083 const G4Element* element ,
00084 const G4Material* material )
00085 {
00086 if ( dp->GetKineticEnergy() == ke_cache && element == element_cache && material == material_cache ) return xs_cache;
00087
00088 ke_cache = dp->GetKineticEnergy();
00089 element_cache = element;
00090 material_cache = material;
00091 G4double xs = GetCrossSection( dp , element , material->GetTemperature() );
00092 xs_cache = xs;
00093 return xs;
00094 }
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106 void G4NeutronHPElasticData::BuildPhysicsTable(const G4ParticleDefinition& aP)
00107 {
00108
00109 if(&aP!=G4Neutron::Neutron())
00110 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
00111
00112
00113 if ( getenv( "G4NEUTRONHP_NEGLECT_DOPPLER" ) )
00114 {
00115 G4cout << "Find environment variable of \"G4NEUTRONHP_NEGLECT_DOPPLER\"." << G4endl;
00116 G4cout << "On the fly Doppler broadening will be neglect in the cross section calculation of elastic scattering of neutrons (<20MeV)." << G4endl;
00117 onFlightDB = false;
00118 }
00119
00120 size_t numberOfElements = G4Element::GetNumberOfElements();
00121
00122
00123 if ( theCrossSections == NULL )
00124 theCrossSections = new G4PhysicsTable( numberOfElements );
00125 else
00126 theCrossSections->clearAndDestroy();
00127
00128
00129
00130 static const G4ElementTable *theElementTable = G4Element::GetElementTable();
00131 for( size_t i=0; i<numberOfElements; ++i )
00132 {
00133 G4PhysicsVector* physVec = G4NeutronHPData::
00134 Instance()->MakePhysicsVector((*theElementTable)[i], this);
00135 theCrossSections->push_back(physVec);
00136 }
00137 }
00138
00139 void G4NeutronHPElasticData::DumpPhysicsTable(const G4ParticleDefinition& aP)
00140 {
00141 if(&aP!=G4Neutron::Neutron())
00142 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
00143
00144
00145
00146
00147
00148
00149
00150
00151 G4cout << G4endl;
00152 G4cout << G4endl;
00153 G4cout << "Elastic Cross Section of Neutron HP"<< G4endl;
00154 G4cout << "(Pointwise cross-section at 0 Kelvin.)" << G4endl;
00155 G4cout << G4endl;
00156 G4cout << "Name of Element" << G4endl;
00157 G4cout << "Energy[eV] XS[barn]" << G4endl;
00158 G4cout << G4endl;
00159
00160 size_t numberOfElements = G4Element::GetNumberOfElements();
00161 static const G4ElementTable *theElementTable = G4Element::GetElementTable();
00162
00163 for ( size_t i = 0 ; i < numberOfElements ; ++i )
00164 {
00165
00166 G4cout << (*theElementTable)[i]->GetName() << G4endl;
00167
00168 G4int ie = 0;
00169
00170 for ( ie = 0 ; ie < 130 ; ie++ )
00171 {
00172 G4double eKinetic = 1.0e-5 * std::pow ( 10.0 , ie/10.0 ) *eV;
00173 G4bool outOfRange = false;
00174
00175 if ( eKinetic < 20*MeV )
00176 {
00177 G4cout << eKinetic/eV << " " << (*((*theCrossSections)(i))).GetValue(eKinetic, outOfRange)/barn << G4endl;
00178 }
00179
00180 }
00181
00182 G4cout << G4endl;
00183 }
00184
00185
00186
00187 }
00188
00189 #include "G4Nucleus.hh"
00190 #include "G4NucleiProperties.hh"
00191 #include "G4Neutron.hh"
00192 #include "G4Electron.hh"
00193
00194 G4double G4NeutronHPElasticData::
00195 GetCrossSection(const G4DynamicParticle* aP, const G4Element*anE, G4double aT)
00196 {
00197 G4double result = 0;
00198 G4bool outOfRange;
00199 G4int index = anE->GetIndex();
00200
00201
00202 G4double eKinetic = aP->GetKineticEnergy();
00203
00204
00205
00206
00207 if ( !onFlightDB )
00208 {
00209 G4double factor = 1.0;
00210 if ( eKinetic < aT * k_Boltzmann )
00211 {
00212
00213
00214
00215
00216 }
00217 return ( (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) )* factor;
00218 }
00219
00220 G4ReactionProduct theNeutron( aP->GetDefinition() );
00221 theNeutron.SetMomentum( aP->GetMomentum() );
00222 theNeutron.SetKineticEnergy( eKinetic );
00223
00224
00225 G4Nucleus aNuc;
00226 G4double eps = 0.0001;
00227 G4double theA = anE->GetN();
00228 G4double theZ = anE->GetZ();
00229 G4double eleMass;
00230
00231
00232 eleMass = ( G4NucleiProperties::GetNuclearMass( static_cast<G4int>(theA+eps) , static_cast<G4int>(theZ+eps) )
00233 ) / G4Neutron::Neutron()->GetPDGMass();
00234
00235 G4ReactionProduct boosted;
00236 G4double aXsection;
00237
00238
00239 G4int counter = 0;
00240 G4double buffer = 0;
00241 G4int size = G4int(std::max(10., aT/60*kelvin));
00242 G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutron.GetMomentum();
00243 G4double neutronVMag = neutronVelocity.mag();
00244
00245 while(counter == 0 || std::abs(buffer-result/std::max(1,counter)) > 0.03*buffer)
00246 {
00247 if(counter) buffer = result/counter;
00248 while (counter<size)
00249 {
00250 counter ++;
00251 G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus(eleMass, aT);
00252 boosted.Lorentz(theNeutron, aThermalNuc);
00253 G4double theEkin = boosted.GetKineticEnergy();
00254 aXsection = (*((*theCrossSections)(index))).GetValue(theEkin, outOfRange);
00255
00256 G4ThreeVector targetVelocity = 1./aThermalNuc.GetMass()*aThermalNuc.GetMomentum();
00257 aXsection *= (targetVelocity-neutronVelocity).mag()/neutronVMag;
00258 result += aXsection;
00259 }
00260 size += size;
00261 }
00262 result /= counter;
00263
00264
00265
00266
00267
00268
00269 return result;
00270 }