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 #include "G4NeutronHPFissionData.hh"
00037 #include "G4SystemOfUnits.hh"
00038 #include "G4Neutron.hh"
00039 #include "G4ElementTable.hh"
00040 #include "G4NeutronHPData.hh"
00041
00042 G4NeutronHPFissionData::G4NeutronHPFissionData()
00043 :G4VCrossSectionDataSet("NeutronHPFissionXS")
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 theCrossSections = 0;
00054 BuildPhysicsTable(*G4Neutron::Neutron());
00055 }
00056
00057 G4NeutronHPFissionData::~G4NeutronHPFissionData()
00058 {
00059 if ( theCrossSections != NULL ) theCrossSections->clearAndDestroy();
00060 delete theCrossSections;
00061 }
00062
00063 G4bool G4NeutronHPFissionData::IsIsoApplicable( const G4DynamicParticle* dp ,
00064 G4int , G4int ,
00065 const G4Element* ,
00066 const G4Material* )
00067 {
00068 G4double eKin = dp->GetKineticEnergy();
00069 if ( eKin > GetMaxKinEnergy()
00070 || eKin < GetMinKinEnergy()
00071 || dp->GetDefinition() != G4Neutron::Neutron() ) return false;
00072
00073 return true;
00074 }
00075
00076 G4double G4NeutronHPFissionData::GetIsoCrossSection( const G4DynamicParticle* dp ,
00077 G4int , G4int ,
00078 const G4Isotope* ,
00079 const G4Element* element ,
00080 const G4Material* material )
00081 {
00082 if ( dp->GetKineticEnergy() == ke_cache && element == element_cache && material == material_cache ) return xs_cache;
00083
00084 ke_cache = dp->GetKineticEnergy();
00085 element_cache = element;
00086 material_cache = material;
00087 G4double xs = GetCrossSection( dp , element , material->GetTemperature() );
00088 xs_cache = xs;
00089 return xs;
00090 }
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102 void G4NeutronHPFissionData::BuildPhysicsTable(const G4ParticleDefinition& aP)
00103 {
00104 if(&aP!=G4Neutron::Neutron())
00105 throw G4HadronicException(__FILE__, __LINE__, "Attempt to use NeutronHP data for particles other than neutrons!!!");
00106 size_t numberOfElements = G4Element::GetNumberOfElements();
00107
00108
00109
00110 if ( theCrossSections == NULL )
00111 theCrossSections = new G4PhysicsTable( numberOfElements );
00112 else
00113 theCrossSections->clearAndDestroy();
00114
00115
00116
00117 static const G4ElementTable *theElementTable = G4Element::GetElementTable();
00118 for( size_t i=0; i<numberOfElements; ++i )
00119 {
00120 G4PhysicsVector* physVec = G4NeutronHPData::
00121 Instance()->MakePhysicsVector((*theElementTable)[i], this);
00122 theCrossSections->push_back(physVec);
00123 }
00124 }
00125
00126 void G4NeutronHPFissionData::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
00136
00137
00138 G4cout << G4endl;
00139 G4cout << G4endl;
00140 G4cout << "Fission Cross Section of Neutron HP"<< G4endl;
00141 G4cout << "(Pointwise cross-section at 0 Kelvin.)" << G4endl;
00142 G4cout << G4endl;
00143 G4cout << "Name of Element" << G4endl;
00144 G4cout << "Energy[eV] XS[barn]" << G4endl;
00145 G4cout << G4endl;
00146
00147 size_t numberOfElements = G4Element::GetNumberOfElements();
00148 static const G4ElementTable *theElementTable = G4Element::GetElementTable();
00149
00150 for ( size_t i = 0 ; i < numberOfElements ; ++i )
00151 {
00152
00153 G4cout << (*theElementTable)[i]->GetName() << G4endl;
00154
00155 if ( (*((*theCrossSections)(i))).GetVectorLength() == 0 )
00156 {
00157 G4cout << "The cross-section data of the fission of this element is not available." << G4endl;
00158 G4cout << G4endl;
00159 continue;
00160 }
00161
00162 G4int ie = 0;
00163
00164 for ( ie = 0 ; ie < 130 ; ie++ )
00165 {
00166 G4double eKinetic = 1.0e-5 * std::pow ( 10.0 , ie/10.0 ) *eV;
00167 G4bool outOfRange = false;
00168
00169 if ( eKinetic < 20*MeV )
00170 {
00171 G4cout << eKinetic/eV << " " << (*((*theCrossSections)(i))).GetValue(eKinetic, outOfRange)/barn << G4endl;
00172 }
00173
00174 }
00175
00176 G4cout << G4endl;
00177 }
00178
00179
00180 }
00181
00182 #include "G4NucleiProperties.hh"
00183
00184 G4double G4NeutronHPFissionData::
00185 GetCrossSection(const G4DynamicParticle* aP, const G4Element*anE, G4double aT)
00186 {
00187 G4double result = 0;
00188 if(anE->GetZ()<90) return result;
00189 G4bool outOfRange;
00190 G4int index = anE->GetIndex();
00191
00192
00193 if ( ( ( *theCrossSections )( index ) )->GetVectorLength() == 0 ) return result;
00194
00195
00196 G4double eKinetic = aP->GetKineticEnergy();
00197 G4ReactionProduct theNeutron( aP->GetDefinition() );
00198 theNeutron.SetMomentum( aP->GetMomentum() );
00199 theNeutron.SetKineticEnergy( eKinetic );
00200
00201
00202 G4Nucleus aNuc;
00203 G4double eps = 0.0001;
00204 G4double theA = anE->GetN();
00205 G4double theZ = anE->GetZ();
00206 G4double eleMass;
00207 eleMass = ( G4NucleiProperties::GetNuclearMass( static_cast<G4int>(theA+eps) , static_cast<G4int>(theZ+eps) )
00208 ) / G4Neutron::Neutron()->GetPDGMass();
00209
00210 G4ReactionProduct boosted;
00211 G4double aXsection;
00212
00213
00214 G4int counter = 0;
00215 G4double buffer = 0;
00216 G4int size = G4int(std::max(10., aT/60*kelvin));
00217 G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutron.GetMomentum();
00218 G4double neutronVMag = neutronVelocity.mag();
00219
00220 while(counter == 0 || std::abs(buffer-result/std::max(1,counter)) > 0.01*buffer)
00221 {
00222 if(counter) buffer = result/counter;
00223 while (counter<size)
00224 {
00225 counter ++;
00226 G4ReactionProduct aThermalNuc = aNuc.GetThermalNucleus(eleMass, aT);
00227 boosted.Lorentz(theNeutron, aThermalNuc);
00228 G4double theEkin = boosted.GetKineticEnergy();
00229 aXsection = (*((*theCrossSections)(index))).GetValue(theEkin, outOfRange);
00230
00231 G4ThreeVector targetVelocity = 1./aThermalNuc.GetMass()*aThermalNuc.GetMomentum();
00232 aXsection *= (targetVelocity-neutronVelocity).mag()/neutronVMag;
00233 result += aXsection;
00234 }
00235 size += size;
00236 }
00237 result /= counter;
00238 return result;
00239 }