Geant4_10
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 76889 2013-11-18 13:01:55Z 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 using namespace std;
57 
59  : G4VCrossSectionDataSet("G4NeutronElasticXS"),
60  proton(G4Proton::Proton()), maxZ(92)
61 {
62  // verboseLevel = 0;
63  if(verboseLevel > 0){
64  G4cout << "G4NeutronElasticXS::G4NeutronElasticXS Initialise for Z < "
65  << maxZ + 1 << G4endl;
66  }
67  data.resize(maxZ+1, 0);
68  coeff.resize(maxZ+1, 1.0);
69  ggXsection = new G4GlauberGribovCrossSection();
70  fNucleon = new G4HadronNucleonXsc();
71  isInitialized = false;
72 }
73 
75 {
76  delete fNucleon;
77  /*
78  for(G4int i=0; i<=maxZ; ++i) {
79  delete data[i];
80  }
81  */
82 }
83 
85 {
86  outFile << "G4NeutronElasticXS calculates the neutron elastic scattering\n"
87  << "cross section on nuclei using data from the high precision\n"
88  << "neutron database. These data are simplified and smoothed over\n"
89  << "the resonance region in order to reduce CPU time.\n"
90  << "G4NeutronElasticXS is valid for energies up to 20 MeV, for all\n"
91  << "targets through U.\n";
92 }
93 
94 G4bool
96  G4int, const G4Material*)
97 {
98  return true;
99 }
100 
101 G4double
103  G4int Z, const G4Material*)
104 {
105  G4double xs = 0.0;
106  G4double ekin = aParticle->GetKineticEnergy();
107 
108  if(Z < 1 || Z > maxZ) { return xs; }
109 
110  G4int Amean = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
111  G4PhysicsVector* pv = data[Z];
112  // G4cout << "G4NeutronElasticXS::GetCrossSection e= " << ekin
113  // << " Z= " << Z << G4endl;
114 
115  // element was not initialised
116  if(!pv) {
117  Initialise(Z);
118  pv = data[Z];
119  if(!pv) { return xs; }
120  }
121 
122  G4double e1 = pv->Energy(0);
123  if(ekin <= e1) { return (*pv)[0]; }
124 
125  G4int n = pv->GetVectorLength() - 1;
126  G4double e2 = pv->Energy(n);
127 
128  if(ekin <= e2) {
129  xs = pv->Value(ekin);
130  } else if(1 == Z) {
131  fNucleon->GetHadronNucleonXscPDG(aParticle, proton);
132  xs = coeff[1]*fNucleon->GetElasticHadronNucleonXsc();
133  } else {
134  ggXsection->GetIsoCrossSection(aParticle, Z, Amean);
135  xs = coeff[Z]*ggXsection->GetElasticGlauberGribovXsc();
136  }
137 
138  if(verboseLevel > 0){
139  G4cout << "ekin= " << ekin << ", XSinel= " << xs << G4endl;
140  }
141  return xs;
142 }
143 
144 void
146 {
147  if(isInitialized) { return; }
148  if(verboseLevel > 0){
149  G4cout << "G4NeutronElasticXS::BuildPhysicsTable for "
150  << p.GetParticleName() << G4endl;
151  }
152  if(p.GetParticleName() != "neutron") {
154  ed << p.GetParticleName() << " is a wrong particle type -"
155  << " only neutron is allowed";
156  G4Exception("G4NeutronElasticXS::BuildPhysicsTable(..)","had012",
157  FatalException, ed, "");
158  return;
159  }
160  isInitialized = true;
161 
162  // check environment variable
163  // Build the complete string identifying the file with the data set
164  char* path = getenv("G4NEUTRONXSDATA");
165 
166  G4DynamicParticle* dynParticle =
168 
169  // Access to elements
170  const G4ElementTable* theElmTable = G4Element::GetElementTable();
171  size_t numOfElm = G4Element::GetNumberOfElements();
172  if(numOfElm > 0) {
173  for(size_t i=0; i<numOfElm; ++i) {
174  G4int Z = G4int(((*theElmTable)[i])->GetZ());
175  if(Z < 1) { Z = 1; }
176  else if(Z > maxZ) { Z = maxZ; }
177  //G4cout << "Z= " << Z << G4endl;
178  // Initialisation
179  if(!data[Z]) { Initialise(Z, dynParticle, path); }
180  }
181  }
182  delete dynParticle;
183 }
184 
185 void
186 G4NeutronElasticXS::Initialise(G4int Z, G4DynamicParticle* dp,
187  const char* p)
188 {
189  if(data[Z]) { return; }
190  const char* path = p;
191  if(!p) {
192  // check environment variable
193  // Build the complete string identifying the file with the data set
194  path = getenv("G4NEUTRONXSDATA");
195  if (!path) {
196  G4Exception("G4NeutronElasticXS::Initialise(..)","had013",
198  "Environment variable G4NEUTRONXSDATA is not defined");
199  return;
200  }
201  }
202  G4DynamicParticle* dynParticle = dp;
203  if(!dp) {
204  dynParticle =
206  }
207 
208  G4int Amean = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
209 
210  // upload data from file
211  data[Z] = new G4PhysicsLogVector();
212 
213  std::ostringstream ost;
214  ost << path << "/elast" << Z ;
215  std::ifstream filein(ost.str().c_str());
216  if (!(filein)) {
218  ed << "Data file <" << ost.str().c_str()
219  << "> is not opened!";
220  G4Exception("G4NeutronElasticXS::Initialise(..)","had014",
221  FatalException, ed, "Check G4NEUTRONXSDATA");
222  return;
223  }else{
224  if(verboseLevel > 1) {
225  G4cout << "file " << ost.str()
226  << " is opened by G4NeutronElasticXS" << G4endl;
227  }
228 
229  // retrieve data from DB
230  if(!data[Z]->Retrieve(filein, true)) {
232  ed << "Data file <" << ost.str().c_str()
233  << "> is not retrieved!";
234  G4Exception("G4NeutronElasticXS::Initialise(..)","had015",
235  FatalException, ed, "Check G4NEUTRONXSDATA");
236  return;
237  }
238 
239  // smooth transition
240  size_t n = data[Z]->GetVectorLength() - 1;
241  G4double emax = data[Z]->Energy(n);
242  G4double sig1 = (*data[Z])[n];
243  dynParticle->SetKineticEnergy(emax);
244  G4double sig2 = 0.0;
245  if(1 == Z) {
246  fNucleon->GetHadronNucleonXscPDG(dynParticle, proton);
247  sig2 = fNucleon->GetElasticHadronNucleonXsc();
248  } else {
249  ggXsection->GetIsoCrossSection(dynParticle, Z, Amean);
250  sig2 = ggXsection->GetElasticGlauberGribovXsc();
251  }
252  if(sig2 > 0.) { coeff[Z] = sig1/sig2; }
253  }
254  if(!dp) { delete dynParticle; }
255 }
G4double GetElasticHadronNucleonXsc()
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
std::ofstream outFile
Definition: GammaRayTel.cc:68
const char * p
Definition: xmltok.h:285
virtual G4double GetIsoCrossSection(const G4DynamicParticle *, G4int Z, G4int A, const G4Isotope *iso=0, const G4Element *elm=0, const G4Material *mat=0)
size_t GetVectorLength() const
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
const G4String & GetParticleName() const
Char_t n[5]
G4GLOB_DLL std::ostream G4cout
static size_t GetNumberOfElements()
Definition: G4Element.cc:402
Float_t Z
Definition: plot.C:39
bool G4bool
Definition: G4Types.hh:79
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4double Energy(size_t index) const
G4double Value(G4double theEnergy, size_t &lastidx) const
void SetKineticEnergy(G4double aEnergy)
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
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:395
const XML_Char const XML_Char * data
Definition: expat.h:268
virtual void BuildPhysicsTable(const G4ParticleDefinition &)