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
00038
00039
00040 #include "G4HadronicException.hh"
00041 #include "G4NeutronInelasticXS.hh"
00042 #include "G4Neutron.hh"
00043 #include "G4DynamicParticle.hh"
00044 #include "G4Element.hh"
00045 #include "G4ElementTable.hh"
00046 #include "G4PhysicsLogVector.hh"
00047 #include "G4PhysicsVector.hh"
00048 #include "G4GlauberGribovCrossSection.hh"
00049 #include "G4HadronNucleonXsc.hh"
00050 #include "G4NistManager.hh"
00051 #include "G4Proton.hh"
00052 #include "Randomize.hh"
00053
00054 #include <iostream>
00055 #include <fstream>
00056 #include <sstream>
00057
00058 using namespace std;
00059
00060 G4NeutronInelasticXS::G4NeutronInelasticXS()
00061 : G4VCrossSectionDataSet("G4NeutronInelasticXS"),
00062 proton(G4Proton::Proton()), maxZ(92)
00063 {
00064
00065 if(verboseLevel > 0){
00066 G4cout << "G4NeutronInelasticXS::G4NeutronInelasticXS Initialise for Z < "
00067 << maxZ + 1 << G4endl;
00068 }
00069 data.resize(maxZ+1, 0);
00070 coeff.resize(maxZ+1, 1.0);
00071 ggXsection = new G4GlauberGribovCrossSection();
00072 fNucleon = new G4HadronNucleonXsc();
00073 isInitialized = false;
00074 }
00075
00076 G4NeutronInelasticXS::~G4NeutronInelasticXS()
00077 {
00078 delete fNucleon;
00079 for(G4int i=0; i<=maxZ; ++i) {
00080 delete data[i];
00081 }
00082 }
00083
00084 void G4NeutronInelasticXS::CrossSectionDescription(std::ostream& outFile) const
00085 {
00086 outFile << "G4NeutronInelasticXS calculates the neutron inelastic scattering\n"
00087 << "cross section on nuclei using data from the high precision\n"
00088 << "neutron database. These data are simplified and smoothed over\n"
00089 << "the resonance region in order to reduce CPU time.\n"
00090 << "G4NeutronInelasticXS is valid for energies up to 20 MeV, for\n"
00091 << "nuclei through U.\n";
00092 }
00093
00094 G4bool
00095 G4NeutronInelasticXS::IsElementApplicable(const G4DynamicParticle*,
00096 G4int, const G4Material*)
00097 {
00098 return true;
00099 }
00100
00101 G4bool
00102 G4NeutronInelasticXS::IsIsoApplicable(const G4DynamicParticle*,
00103 G4int , G4int ,
00104 const G4Element*, const G4Material*)
00105 {
00106 return false;
00107 }
00108
00109 G4double
00110 G4NeutronInelasticXS::GetElementCrossSection(const G4DynamicParticle* aParticle,
00111 G4int Z, const G4Material*)
00112 {
00113 G4double xs = 0.0;
00114 G4double ekin = aParticle->GetKineticEnergy();
00115
00116 if(Z < 1 || Z > maxZ) { return xs; }
00117 G4int Amean = G4int(G4NistManager::Instance()->GetAtomicMassAmu(Z)+0.5);
00118 G4PhysicsVector* pv = data[Z];
00119
00120
00121
00122 if(!pv) {
00123 Initialise(Z);
00124 pv = data[Z];
00125 if(!pv) { return xs; }
00126 }
00127
00128 G4double e1 = pv->Energy(0);
00129 if(ekin <= e1) { return xs; }
00130
00131 G4int n = pv->GetVectorLength() - 1;
00132 G4double e2 = pv->Energy(n);
00133 if(ekin <= e2) {
00134 xs = pv->Value(ekin);
00135 } else if(1 == Z) {
00136 fNucleon->GetHadronNucleonXscPDG(aParticle, proton);
00137 xs = coeff[1]*fNucleon->GetInelasticHadronNucleonXsc();
00138 } else {
00139 ggXsection->GetIsoCrossSection(aParticle, Z, Amean);
00140 xs = coeff[Z]*ggXsection->GetInelasticGlauberGribovXsc();
00141 }
00142
00143 if(verboseLevel > 0) {
00144 G4cout << "ekin= " << ekin << ", XSinel= " << xs << G4endl;
00145 }
00146 return xs;
00147 }
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228 void
00229 G4NeutronInelasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
00230 {
00231 if(isInitialized) { return; }
00232 if(verboseLevel > 0){
00233 G4cout << "G4NeutronCaptureXS::BuildPhysicsTable for "
00234 << p.GetParticleName() << G4endl;
00235 }
00236 if(p.GetParticleName() != "neutron") {
00237 throw G4HadronicException(__FILE__, __LINE__,"Wrong particle type");
00238 return;
00239 }
00240 isInitialized = true;
00241
00242
00243
00244 char* path = getenv("G4NEUTRONXSDATA");
00245 if (!path){
00246 throw G4HadronicException(__FILE__, __LINE__,
00247 "G4NEUTRONXSDATA environment variable not defined");
00248 return;
00249 }
00250
00251 G4DynamicParticle* dynParticle =
00252 new G4DynamicParticle(G4Neutron::Neutron(),G4ThreeVector(1,0,0),1);
00253
00254
00255 const G4ElementTable* theElmTable = G4Element::GetElementTable();
00256 size_t numOfElm = G4Element::GetNumberOfElements();
00257 if(numOfElm > 0) {
00258 for(size_t i=0; i<numOfElm; ++i) {
00259 G4int Z = G4int(((*theElmTable)[i])->GetZ());
00260 if(Z < 1) { Z = 1; }
00261 else if(Z > maxZ) { Z = maxZ; }
00262
00263
00264 if(!data[Z]) { Initialise(Z, dynParticle, path); }
00265 }
00266 }
00267 delete dynParticle;
00268 }
00269
00270 void
00271 G4NeutronInelasticXS::Initialise(G4int Z, G4DynamicParticle* dp,
00272 const char* p)
00273 {
00274 if(data[Z]) { return; }
00275 const char* path = p;
00276 if(!p) {
00277
00278
00279 path = getenv("G4NEUTRONXSDATA");
00280 if (!path) {
00281 throw G4HadronicException(__FILE__, __LINE__,
00282 "G4NEUTRONXSDATA environment variable not defined");
00283 return;
00284 }
00285 }
00286 G4DynamicParticle* dynParticle = dp;
00287 if(!dp) {
00288 dynParticle =
00289 new G4DynamicParticle(G4Neutron::Neutron(),G4ThreeVector(1,0,0),1);
00290 }
00291
00292 G4int Amean = G4int(G4NistManager::Instance()->GetAtomicMassAmu(Z)+0.5);
00293
00294
00295 data[Z] = new G4PhysicsLogVector();
00296
00297 std::ostringstream ost;
00298 ost << path << "/inelast" << Z ;
00299 std::ifstream filein(ost.str().c_str());
00300
00301 if (!(filein)) {
00302 throw G4HadronicException(__FILE__, __LINE__,"NO data sets opened");
00303 return;
00304 }else{
00305 if(verboseLevel > 1) {
00306 G4cout << "file " << ost.str()
00307 << " is opened by G4NeutronInelasticXS" << G4endl;
00308 }
00309
00310
00311 data[Z]->Retrieve(filein, true);
00312
00313
00314 size_t n = data[Z]->GetVectorLength() - 1;
00315 G4double emax = data[Z]->Energy(n);
00316 G4double sig1 = (*data[Z])[n];
00317 dynParticle->SetKineticEnergy(emax);
00318 G4double sig2 = 0.0;
00319 if(1 == Z) {
00320 fNucleon->GetHadronNucleonXscPDG(dynParticle, proton);
00321 sig2 = fNucleon->GetInelasticHadronNucleonXsc();
00322 } else {
00323 ggXsection->GetIsoCrossSection(dynParticle, Z, Amean);
00324 sig2 = ggXsection->GetInelasticGlauberGribovXsc();
00325 }
00326 if(sig2 > 0.) { coeff[Z] = sig1/sig2; }
00327 }
00328 if(!dp) { delete dynParticle; }
00329 }