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 "G4NeutronHPInelasticData.hh"
00038 #include "G4PhysicalConstants.hh"
00039 #include "G4SystemOfUnits.hh"
00040 #include "G4Neutron.hh"
00041 #include "G4ElementTable.hh"
00042 #include "G4NeutronHPData.hh"
00043
00044 G4NeutronHPInelasticData::G4NeutronHPInelasticData()
00045 :G4VCrossSectionDataSet("NeutronHPInelasticXS")
00046 {
00047
00048 SetMinKinEnergy( 0*MeV );
00049 SetMaxKinEnergy( 20*MeV );
00050
00051 ke_cache = 0.0;
00052 xs_cache = 0.0;
00053 element_cache = NULL;
00054 material_cache = NULL;
00055
00056 onFlightDB = true;
00057 theCrossSections = 0;
00058 BuildPhysicsTable(*G4Neutron::Neutron());
00059 }
00060
00061 G4NeutronHPInelasticData::~G4NeutronHPInelasticData()
00062 {
00063 if ( theCrossSections != 0 ) theCrossSections->clearAndDestroy();
00064 delete theCrossSections;
00065 }
00066
00067 G4bool G4NeutronHPInelasticData::IsIsoApplicable( const G4DynamicParticle* dp ,
00068 G4int , G4int ,
00069 const G4Element* ,
00070 const G4Material* )
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 G4NeutronHPInelasticData::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
00107 void G4NeutronHPInelasticData::BuildPhysicsTable(const G4ParticleDefinition& aP)
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 inelastic scattering of neutrons (<20MeV)." << G4endl;
00117 onFlightDB = false;
00118 }
00119
00120 size_t numberOfElements = G4Element::GetNumberOfElements();
00121
00122
00123
00124
00125 if ( theCrossSections == NULL )
00126 theCrossSections = new G4PhysicsTable( numberOfElements );
00127 else
00128 theCrossSections->clearAndDestroy();
00129
00130
00131
00132 static const G4ElementTable *theElementTable = G4Element::GetElementTable();
00133 for( size_t i=0; i<numberOfElements; ++i )
00134 {
00135 G4PhysicsVector* physVec = G4NeutronHPData::
00136 Instance()->MakePhysicsVector((*theElementTable)[i], this);
00137 theCrossSections->push_back(physVec);
00138 }
00139 }
00140
00141 void G4NeutronHPInelasticData::DumpPhysicsTable(const G4ParticleDefinition& aP)
00142 {
00143 if(&aP!=G4Neutron::Neutron())
00144 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
00145
00146
00147
00148
00149
00150
00151
00152
00153 G4cout << G4endl;
00154 G4cout << G4endl;
00155 G4cout << "Inelastic Cross Section of Neutron HP"<< G4endl;
00156 G4cout << "(Pointwise cross-section at 0 Kelvin.)" << G4endl;
00157 G4cout << G4endl;
00158 G4cout << "Name of Element" << G4endl;
00159 G4cout << "Energy[eV] XS[barn]" << G4endl;
00160 G4cout << G4endl;
00161
00162 size_t numberOfElements = G4Element::GetNumberOfElements();
00163 static const G4ElementTable *theElementTable = G4Element::GetElementTable();
00164
00165 for ( size_t i = 0 ; i < numberOfElements ; ++i )
00166 {
00167
00168 G4cout << (*theElementTable)[i]->GetName() << G4endl;
00169
00170 G4int ie = 0;
00171
00172 for ( ie = 0 ; ie < 130 ; ie++ )
00173 {
00174 G4double eKinetic = 1.0e-5 * std::pow ( 10.0 , ie/10.0 ) *eV;
00175 G4bool outOfRange = false;
00176
00177 if ( eKinetic < 20*MeV )
00178 {
00179 G4cout << eKinetic/eV << " " << (*((*theCrossSections)(i))).GetValue(eKinetic, outOfRange)/barn << G4endl;
00180 }
00181
00182 }
00183
00184 G4cout << G4endl;
00185 }
00186
00187
00188 }
00189
00190 #include "G4NucleiProperties.hh"
00191
00192 G4double G4NeutronHPInelasticData::
00193 GetCrossSection(const G4DynamicParticle* aP, const G4Element*anE, G4double aT)
00194 {
00195 G4double result = 0;
00196 G4bool outOfRange;
00197 G4int index = anE->GetIndex();
00198
00199
00200 G4double eKinetic = aP->GetKineticEnergy();
00201
00202
00203
00204
00205 if ( !onFlightDB )
00206 {
00207 G4double factor = 1.0;
00208 if ( eKinetic < aT * k_Boltzmann )
00209 {
00210
00211
00212
00213
00214 }
00215 return ( (*((*theCrossSections)(index))).GetValue(eKinetic, outOfRange) )* factor;
00216 }
00217
00218 G4ReactionProduct theNeutron( aP->GetDefinition() );
00219 theNeutron.SetMomentum( aP->GetMomentum() );
00220 theNeutron.SetKineticEnergy( eKinetic );
00221
00222
00223 G4Nucleus aNuc;
00224 G4double eps = 0.0001;
00225 G4double theA = anE->GetN();
00226 G4double theZ = anE->GetZ();
00227 G4double eleMass;
00228 eleMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theA+eps), static_cast<G4int>(theZ+eps))
00229 ) / G4Neutron::Neutron()->GetPDGMass();
00230
00231 G4ReactionProduct boosted;
00232 G4double aXsection;
00233
00234
00235 G4int counter = 0;
00236 G4int failCount = 0;
00237 G4double buffer = 0;
00238 G4int size = G4int(std::max(10., aT/60*kelvin));
00239 G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutron.GetMomentum();
00240 G4double neutronVMag = neutronVelocity.mag();
00241
00242 while(counter == 0 || std::abs(buffer-result/std::max(1,counter)) > 0.01*buffer)
00243 {
00244 if(counter) buffer = result/counter;
00245 while (counter<size)
00246 {
00247 counter ++;
00248 G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus(eleMass, aT);
00249 boosted.Lorentz(theNeutron, aThermalNuc);
00250 G4double theEkin = boosted.GetKineticEnergy();
00251 aXsection = (*((*theCrossSections)(index))).GetValue(theEkin, outOfRange);
00252 if(aXsection <0)
00253 {
00254 if(failCount<1000)
00255 {
00256 failCount++;
00257 counter--;
00258 continue;
00259 }
00260 else
00261 {
00262 aXsection = 0;
00263 }
00264 }
00265
00266 G4ThreeVector targetVelocity = 1./aThermalNuc.GetMass()*aThermalNuc.GetMomentum();
00267 aXsection *= (targetVelocity-neutronVelocity).mag()/neutronVMag;
00268 result += aXsection;
00269 }
00270 size += size;
00271 }
00272 result /= counter;
00273
00274
00275
00276
00277
00278
00279 return result;
00280 }