Geant4  10.01.p02
G4NeutronElasticXS.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4NeutronElasticXS.cc 83697 2014-09-10 07:15:29Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4NeutronElasticXS
34 //
35 // Author Ivantchenko, Geant4, 3-Aug-09
36 //
37 // Modifications:
38 //
39 
40 #include "G4NeutronElasticXS.hh"
41 #include "G4Neutron.hh"
42 #include "G4DynamicParticle.hh"
43 #include "G4Element.hh"
44 #include "G4ElementTable.hh"
45 #include "G4PhysicsLogVector.hh"
46 #include "G4PhysicsVector.hh"
48 #include "G4HadronNucleonXsc.hh"
49 #include "G4NistManager.hh"
50 #include "G4Proton.hh"
51 
52 #include <iostream>
53 #include <fstream>
54 #include <sstream>
55 
56 // factory
57 #include "G4CrossSectionFactory.hh"
58 //
60 
61 using namespace std;
62 
63 std::vector<G4PhysicsVector*>* G4NeutronElasticXS::data = 0;
65 
67  : G4VCrossSectionDataSet(Default_Name()),
69 {
70  // verboseLevel = 0;
71  if(verboseLevel > 0){
72  G4cout << "G4NeutronElasticXS::G4NeutronElasticXS Initialise for Z < "
73  << MAXZEL << G4endl;
74  }
77  isMaster = false;
78 }
79 
81 {
82  delete fNucleon;
83  delete ggXsection;
84  if(isMaster) {
85  for(G4int i=0; i<MAXZEL; ++i) {
86  delete (*data)[i];
87  (*data)[i] = 0;
88  }
89  delete data;
90  data = 0;
91  }
92 }
93 
94 void G4NeutronElasticXS::CrossSectionDescription(std::ostream& outFile) const
95 {
96  outFile << "G4NeutronElasticXS calculates the neutron elastic scattering\n"
97  << "cross section on nuclei using data from the high precision\n"
98  << "neutron database. These data are simplified and smoothed over\n"
99  << "the resonance region in order to reduce CPU time.\n"
100  << "G4NeutronElasticXS is valid for energies up to 20 MeV, for all\n"
101  << "targets through U.\n";
102 }
103 
104 G4bool
106  G4int, const G4Material*)
107 {
108  return true;
109 }
110 
111 G4double
113  G4int Z, const G4Material*)
114 {
115  G4double xs = 0.0;
116  G4double ekin = aParticle->GetKineticEnergy();
117 
118  if(Z < 1 || Z >=MAXZEL) { return xs; }
119 
120  G4int Amean = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
121  G4PhysicsVector* pv = (*data)[Z];
122  // G4cout << "G4NeutronElasticXS::GetCrossSection e= " << ekin
123  // << " Z= " << Z << G4endl;
124 
125  // element was not initialised
126  if(!pv) {
127  Initialise(Z);
128  pv = (*data)[Z];
129  if(!pv) { return xs; }
130  }
131 
132  G4double e1 = pv->Energy(0);
133  if(ekin <= e1) { return (*pv)[0]; }
134 
135  G4int n = pv->GetVectorLength() - 1;
136  G4double e2 = pv->Energy(n);
137 
138  if(ekin <= e2) {
139  xs = pv->Value(ekin);
140  } else if(1 == Z) {
143  } else {
144  ggXsection->GetIsoCrossSection(aParticle, Z, Amean);
146  }
147 
148  if(verboseLevel > 0){
149  G4cout << "ekin= " << ekin << ", XSinel= " << xs << G4endl;
150  }
151  return xs;
152 }
153 
154 void
156 {
157  if(verboseLevel > 0){
158  G4cout << "G4NeutronElasticXS::BuildPhysicsTable for "
159  << p.GetParticleName() << G4endl;
160  }
161  if(p.GetParticleName() != "neutron") {
163  ed << p.GetParticleName() << " is a wrong particle type -"
164  << " only neutron is allowed";
165  G4Exception("G4NeutronElasticXS::BuildPhysicsTable(..)","had012",
166  FatalException, ed, "");
167  return;
168  }
169 
170  if(0 == coeff[0]) {
171  isMaster = true;
172  for(G4int i=0; i<MAXZEL; ++i) { coeff[i] = 1.0; }
173  data = new std::vector<G4PhysicsVector*>;
174  data->resize(MAXZEL, 0);
175  }
176 
177  // it is possible re-initialisation for the second run
178  if(isMaster) {
179 
180  // check environment variable
181  // Build the complete string identifying the file with the data set
182  char* path = getenv("G4NEUTRONXSDATA");
183 
184  G4DynamicParticle* dynParticle =
186 
187  // Access to elements
188  const G4ElementTable* theElmTable = G4Element::GetElementTable();
189  size_t numOfElm = G4Element::GetNumberOfElements();
190  if(numOfElm > 0) {
191  for(size_t i=0; i<numOfElm; ++i) {
192  G4int Z = G4int(((*theElmTable)[i])->GetZ());
193  if(Z < 1) { Z = 1; }
194  else if(Z >= MAXZEL) { Z = MAXZEL-1; }
195  //G4cout << "Z= " << Z << G4endl;
196  // Initialisation
197  if(!(*data)[Z]) { Initialise(Z, dynParticle, path); }
198  }
199  }
200  delete dynParticle;
201  }
202 }
203 
204 void
206  const char* p)
207 {
208  if((*data)[Z]) { return; }
209  const char* path = p;
210  if(!p) {
211  // check environment variable
212  // Build the complete string identifying the file with the data set
213  path = getenv("G4NEUTRONXSDATA");
214  if (!path) {
215  G4Exception("G4NeutronElasticXS::Initialise(..)","had013",
217  "Environment variable G4NEUTRONXSDATA is not defined");
218  return;
219  }
220  }
221  G4DynamicParticle* dynParticle = dp;
222  if(!dp) {
223  dynParticle =
225  }
226 
227  G4int Amean = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
228 
229  // upload data from file
230  (*data)[Z] = new G4PhysicsLogVector();
231 
232  std::ostringstream ost;
233  ost << path << "/elast" << Z ;
234  std::ifstream filein(ost.str().c_str());
235  if (!(filein)) {
237  ed << "Data file <" << ost.str().c_str()
238  << "> is not opened!";
239  G4Exception("G4NeutronElasticXS::Initialise(..)","had014",
240  FatalException, ed, "Check G4NEUTRONXSDATA");
241  return;
242  }else{
243  if(verboseLevel > 1) {
244  G4cout << "file " << ost.str()
245  << " is opened by G4NeutronElasticXS" << G4endl;
246  }
247 
248  // retrieve data from DB
249  if(!((*data)[Z]->Retrieve(filein, true))) {
251  ed << "Data file <" << ost.str().c_str()
252  << "> is not retrieved!";
253  G4Exception("G4NeutronElasticXS::Initialise(..)","had015",
254  FatalException, ed, "Check G4NEUTRONXSDATA");
255  return;
256  }
257 
258  // smooth transition
259  size_t n = (*data)[Z]->GetVectorLength() - 1;
260  G4double emax = (*data)[Z]->Energy(n);
261  G4double sig1 = (*(*data)[Z])[n];
262  dynParticle->SetKineticEnergy(emax);
263  G4double sig2 = 0.0;
264  if(1 == Z) {
265  fNucleon->GetHadronNucleonXscPDG(dynParticle, proton);
267  } else {
268  ggXsection->GetIsoCrossSection(dynParticle, Z, Amean);
270  }
271  if(sig2 > 0.) { coeff[Z] = sig1/sig2; }
272  }
273  if(!dp) { delete dynParticle; }
274 }
G4double GetElasticHadronNucleonXsc()
const G4ParticleDefinition * proton
G4HadronNucleonXsc * fNucleon
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
static G4double coeff[MAXZEL]
virtual G4double GetIsoCrossSection(const G4DynamicParticle *, G4int Z, G4int A, const G4Isotope *iso=0, const G4Element *elm=0, const G4Material *mat=0)
static const G4double e2
size_t GetVectorLength() const
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
const G4String & GetParticleName() const
G4GLOB_DLL std::ostream G4cout
static size_t GetNumberOfElements()
Definition: G4Element.cc:410
bool G4bool
Definition: G4Types.hh:79
G4_DECLARE_XS_FACTORY(G4NeutronElasticXS)
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
const G4int n
G4double Energy(size_t index) const
G4double Value(G4double theEnergy, size_t &lastidx) const
void SetKineticEnergy(G4double aEnergy)
static std::vector< G4PhysicsVector * > * data
static const G4double e1
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
int G4lrint(double ad)
Definition: templates.hh:163
G4GlauberGribovCrossSection * ggXsection
virtual G4double GetElementCrossSection(const G4DynamicParticle *, G4int Z, const G4Material *mat=0)
G4double GetHadronNucleonXscPDG(const G4DynamicParticle *, const G4ParticleDefinition *)
#define G4endl
Definition: G4ios.hh:61
virtual void CrossSectionDescription(std::ostream &) const
std::vector< G4Element * > G4ElementTable
virtual G4bool IsElementApplicable(const G4DynamicParticle *, G4int Z, const G4Material *)
double G4double
Definition: G4Types.hh:76
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:403
const G4int MAXZEL
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
void Initialise(G4int Z, G4DynamicParticle *dp=0, const char *=0)