Geant4  10.00.p03
G4NeutronCaptureXS.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: G4NeutronCaptureXS.cc 71324 2013-06-13 16:58:55Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4NeutronCaptureXS
34 //
35 // Author Ivantchenko, Geant4, 3-Aug-09
36 //
37 // Modifications:
38 //
39 
40 #include <fstream>
41 #include <sstream>
42 
43 #include "G4SystemOfUnits.hh"
44 #include "G4NeutronCaptureXS.hh"
45 #include "G4Element.hh"
46 #include "G4ElementTable.hh"
47 #include "G4PhysicsLogVector.hh"
48 #include "G4PhysicsVector.hh"
49 #include "G4DynamicParticle.hh"
50 #include "Randomize.hh"
51 
52 using namespace std;
53 
54 const G4int G4NeutronCaptureXS::amin[] = {0,
55  0, 0, 6, 0,10,12,14,16, 0, 0, //1-10
56  0, 0, 0,28, 0, 0, 0,36, 0,40, //11-20
57  0, 0, 0, 0, 0,54, 0,58,63,64, //21-30
58  0,70, 0, 0, 0, 0, 0, 0, 0,90, //31-40
59  0, 0, 0, 0, 0, 0,107,106, 0,112, //41-50
60  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //51-60
61  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //61-70
62  0, 0, 0,180, 0, 0, 0, 0, 0, 0, //71-80
63  0,204, 0, 0, 0, 0, 0, 0, 0, 0, //81-90
64  0,235};
65 const G4int G4NeutronCaptureXS::amax[] = {0,
66  0, 0, 7, 0,11,13,15,18, 0, 0, //1-10
67  0, 0, 0,30, 0, 0, 0,40, 0,48, //11-20
68  0, 0, 0, 0, 0,58, 0,64,65,70, //21-30
69  0,76, 0, 0, 0, 0, 0, 0, 0,96, //31-40
70  0, 0, 0, 0, 0, 0,109,116, 0,124, //41-50
71  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //51-60
72  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //61-70
73  0, 0, 0,186, 0, 0, 0, 0, 0, 0, //71-80
74  0,208, 0, 0, 0, 0, 0, 0, 0, 0, //81-90
75  0,238};
76 
78  : G4VCrossSectionDataSet("G4NeutronCaptureXS"),
79  emax(20*MeV),elimit(1.0e-10*eV)
80 {
81  // verboseLevel = 0;
82  if(verboseLevel > 0){
83  G4cout << "G4NeutronCaptureXS::G4NeutronCaptureXS: Initialise for Z < "
84  << MAXZCAPTURE << G4endl;
85  }
86  data.SetName("NeutronCapture");
87  work.resize(13,0);
88  temp.resize(13,0.0);
89  isInitialized = false;
90 }
91 
93 {}
94 
96 {
97  outFile << "G4NeutronCaptureXS calculates the neutron capture cross sections\n"
98  << "on nuclei using data from the high precision neutron database.\n"
99  << "These data are simplified and smoothed over the resonance region\n"
100  << "in order to reduce CPU time. G4NeutronCaptureXS is valid up to\n"
101  << "20 MeV for all targets through U.\n";
102 }
103 
104 G4bool
106  G4int, const G4Material*)
107 {
108  return true;
109 }
110 
111 G4bool
113  G4int /*ZZ*/, G4int /*AA*/,
114  const G4Element*, const G4Material*)
115 {
116  return true;
117 }
118 
119 G4double
121  G4int Z, const G4Material*)
122 {
123  G4double xs = 0.0;
124  G4double ekin = aParticle->GetKineticEnergy();
125  if(ekin > emax || Z < 1 || Z >= MAXZCAPTURE) { return xs; }
126  if(ekin < elimit) { ekin = elimit; }
127 
129 
130  // element was not initialised
131  if(!pv) {
132  Initialise(Z);
133  pv = data.GetElementData(Z);
134  if(!pv) { return xs; }
135  }
136 
137  G4double e1 = pv->Energy(0);
138  if(ekin < e1) { xs = (*pv)[0]*std::sqrt(e1/ekin); }
139  else if(ekin <= pv->GetMaxEnergy()) { xs = pv->Value(ekin); }
140 
141  if(verboseLevel > 0){
142  G4cout << "ekin= " << ekin << ", xs= " << xs << G4endl;
143  }
144  return xs;
145 }
146 
147 G4double
149  G4int Z, G4int A,
150  const G4Isotope*, const G4Element*,
151  const G4Material*)
152 {
153  G4double xs = 0.0;
154  G4double ekin = aParticle->GetKineticEnergy();
155  if(ekin <= emax && Z > 0 && Z < MAXZCAPTURE) {
156  xs = IsoCrossSection(ekin, Z, A);
157  }
158  return xs;
159 }
160 
162 {
163  G4double xs = 0.0;
164  if(ekin < elimit) { ekin = elimit; }
165 
167 
168  // element was not initialised
169  if(!pv) {
170  Initialise(Z);
171  pv = data.GetElementData(Z);
172  if(!pv) { return xs; }
173  }
175  if(pviso) { pv = pviso; }
176 
177  G4double e1 = pv->Energy(1);
178  if(ekin < e1) { xs = (*pv)[1]*std::sqrt(e1/ekin); }
179  else if(ekin <= pv->GetMaxEnergy()) { xs = pv->Value(ekin); }
180 
181  if(verboseLevel > 0){
182  G4cout << "G4NeutronCaptureXS::IsoCrossSection: Ekin(MeV)= " << ekin/MeV
183  << " xs(b)= " << xs/barn
184  << " Z= " << Z << " A= " << A << " " << pv->GetVectorLength()
185  << G4endl;
186  }
187  return xs;
188 }
189 
191  G4double kinEnergy)
192 {
193  G4int nIso = anElement->GetNumberOfIsotopes();
194  G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
195  G4Isotope* iso = (*isoVector)[0];
196 
197  // more than 1 isotope
198  if(1 < nIso) {
199  G4int Z = G4lrint(anElement->GetZ());
200  if(Z >= MAXZCAPTURE) { Z = MAXZCAPTURE-1; }
201  G4double* abundVector = anElement->GetRelativeAbundanceVector();
202  G4double q = G4UniformRand();
203  G4double sum = 0.0;
204 
205  // is there isotope wise cross section?
206  if(0 == amin[Z]) {
207  for (G4int j = 0; j<nIso; ++j) {
208  sum += abundVector[j];
209  if(q <= sum) {
210  iso = (*isoVector)[j];
211  break;
212  }
213  }
214  } else {
215  size_t nmax = data.GetNumberOfComponents(Z);
216  if(temp.size() < nmax) { temp.resize(nmax,0.0); }
217  for (size_t i=0; i<nmax; ++i) {
218  G4int A = (*isoVector)[i]->GetN();
219  sum += abundVector[i]*IsoCrossSection(kinEnergy, Z, A);
220  temp[i] = sum;
221  }
222  sum *= q;
223  for (size_t j = 0; j<nmax; ++j) {
224  if(temp[j] >= sum) {
225  iso = (*isoVector)[j];
226  break;
227  }
228  }
229  }
230  }
231  return iso;
232 }
233 
234 void
236 {
237  if(isInitialized) { return; }
238  if(verboseLevel > 0){
239  G4cout << "G4NeutronCaptureXS::BuildPhysicsTable for "
240  << p.GetParticleName() << G4endl;
241  }
242  if(p.GetParticleName() != "neutron") {
244  ed << p.GetParticleName() << " is a wrong particle type -"
245  << " only neutron is allowed";
246  G4Exception("G4NeutronCaptureXS::BuildPhysicsTable(..)","had012",
247  FatalException, ed, "");
248  return;
249  }
250  isInitialized = true;
251 
252  // check environment variable
253  // Build the complete string identifying the file with the data set
254  char* path = getenv("G4NEUTRONXSDATA");
255 
256  // Access to elements
257  const G4ElementTable* theElmTable = G4Element::GetElementTable();
258  size_t numOfElm = G4Element::GetNumberOfElements();
259  if(numOfElm > 0) {
260  for(size_t i=0; i<numOfElm; ++i) {
261  G4int Z = G4int(((*theElmTable)[i])->GetZ());
262  if(Z < 1) { Z = 1; }
263  else if(Z >= MAXZCAPTURE) { Z = MAXZCAPTURE-1; }
264  //G4cout << "Z= " << Z << G4endl;
265  // Initialisation
266  if(!data.GetElementData(Z)) { Initialise(Z, path); }
267  }
268  }
269 }
270 
271 void
273 {
274  if(data.GetElementData(Z)) { return; }
275  const char* path = p;
276 
277  // check environment variable
278  if(!p) {
279  path = getenv("G4NEUTRONXSDATA");
280  if (!path) {
281  G4Exception("G4NeutronCaptureXS::Initialise(..)","had013",FatalException,
282  "Environment variable G4NEUTRONXSDATA is not defined");
283  return;
284  }
285  }
286 
287  // upload element data
288  std::ostringstream ost;
289  ost << path << "/cap" << Z ;
290  G4PhysicsVector* v = RetrieveVector(ost, true);
292 
293  // upload isotope data
294  if(amin[Z] > 0) {
295  size_t n = 0;
296  size_t i = 0;
297  size_t nmax = (size_t)(amax[Z]-amin[Z]+1);
298  if(work.size() < nmax) { work.resize(nmax,0); }
299  for(G4int A=amin[Z]; A<=amax[Z]; ++A) {
300  std::ostringstream ost1;
301  ost1 << path << "/cap" << Z << "_" << A;
302  v = RetrieveVector(ost1, false);
303  if(v) { ++n; }
304  work[i] = v;
305  ++i;
306  }
308  for(size_t j=0; j<i; ++j) {
309  if(work[j]) { data.AddComponent(Z, amin[Z]+j, work[j]); }
310  }
311  }
312 }
313 
315 G4NeutronCaptureXS::RetrieveVector(std::ostringstream& ost, G4bool warn)
316 {
317  G4PhysicsLogVector* v = 0;
318  std::ifstream filein(ost.str().c_str());
319  if (!(filein)) {
320  if(!warn) { return v; }
322  ed << "Data file <" << ost.str().c_str()
323  << "> is not opened!";
324  G4Exception("G4NeutronCaptureXS::RetrieveVector(..)","had014",
325  FatalException, ed, "Check G4NEUTRONXSDATA");
326  }else{
327  if(verboseLevel > 1) {
328  G4cout << "File " << ost.str()
329  << " is opened by G4NeutronCaptureXS" << G4endl;
330  }
331  // retrieve data from DB
332  v = new G4PhysicsLogVector();
333  if(!v->Retrieve(filein, true)) {
335  ed << "Data file <" << ost.str().c_str()
336  << "> is not retrieved!";
337  G4Exception("G4NeutronCaptureXS::RetrieveVector(..)","had015",
338  FatalException, ed, "Check G4NEUTRONXSDATA");
339  }
340  }
341  return v;
342 }
size_t GetNumberOfIsotopes() const
Definition: G4Element.hh:158
static const G4int amin[MAXZCAPTURE]
static const double MeV
Definition: G4SIunits.hh:193
G4PhysicsVector * RetrieveVector(std::ostringstream &in, G4bool warn)
virtual G4bool Retrieve(std::ifstream &fIn, G4bool ascii)
std::vector< G4Isotope * > G4IsotopeVector
void AddComponent(G4int Z, G4int id, G4PhysicsVector *v)
virtual G4double GetElementCrossSection(const G4DynamicParticle *, G4int Z, const G4Material *mat=0)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
virtual G4double GetIsoCrossSection(const G4DynamicParticle *, G4int Z, G4int A, const G4Isotope *iso, const G4Element *elm, const G4Material *mat)
std::ofstream outFile
Definition: GammaRayTel.cc:68
G4PhysicsVector * GetComponentDataByID(G4int Z, G4int id)
std::vector< G4double > temp
G4PhysicsVector * GetElementData(G4int Z)
G4double GetZ() const
Definition: G4Element.hh:131
virtual G4bool IsElementApplicable(const G4DynamicParticle *, G4int Z, const G4Material *)
const G4int MAXZCAPTURE
size_t GetVectorLength() const
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
void InitialiseForElement(G4int Z, G4PhysicsVector *v)
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
static size_t GetNumberOfElements()
Definition: G4Element.cc:410
bool G4bool
Definition: G4Types.hh:79
std::vector< G4PhysicsVector * > work
void Initialise(G4int Z, const char *=0)
G4double * GetRelativeAbundanceVector() const
Definition: G4Element.hh:166
const G4int nmax
const G4int n
G4double Energy(size_t index) const
size_t GetNumberOfComponents(G4int Z)
G4double Value(G4double theEnergy, size_t &lastidx) const
static const G4double A[nN]
G4double IsoCrossSection(G4double ekin, G4int Z, G4int A)
virtual void CrossSectionDescription(std::ostream &) const
static const G4double e1
virtual G4Isotope * SelectIsotope(const G4Element *, G4double kinEnergy)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual G4bool IsIsoApplicable(const G4DynamicParticle *, G4int Z, G4int A, const G4Element *, const G4Material *)
static const double eV
Definition: G4SIunits.hh:194
G4IsotopeVector * GetIsotopeVector() const
Definition: G4Element.hh:162
int G4lrint(double ad)
Definition: templates.hh:163
static const G4int amax[MAXZCAPTURE]
#define G4endl
Definition: G4ios.hh:61
static const double barn
Definition: G4SIunits.hh:95
void SetName(const G4String &nam)
std::vector< G4Element * > G4ElementTable
double G4double
Definition: G4Types.hh:76
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:403
void InitialiseForComponent(G4int Z, G4int nComponents=0)