Geant4_10
G4NeutronInelasticXS.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: G4NeutronInelasticXS.cc 71324 2013-06-13 16:58:55Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4NeutronInelasticXS
34 //
35 // Author Ivantchenko, Geant4, 3-Aug-09
36 //
37 // Modifications:
38 //
39 
40 #include "G4NeutronInelasticXS.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 #include "Randomize.hh"
52 
53 #include <iostream>
54 #include <fstream>
55 #include <sstream>
56 
57 using namespace std;
58 
59 const G4int G4NeutronInelasticXS::amin[] = {0,
60  0, 0, 6, 0,10,12,14,16, 0, 0, //1-10
61  0, 0, 0,28, 0, 0, 0,36, 0,40, //11-20
62  0, 0, 0, 0, 0,54, 0,58,63,64, //21-30
63  0,70, 0, 0, 0, 0, 0, 0, 0,90, //31-40
64  0, 0, 0, 0, 0, 0,107,106, 0,112, //41-50
65  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //51-60
66  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //61-70
67  0, 0, 0,180, 0, 0, 0, 0, 0, 0, //71-80
68  0,204, 0, 0, 0, 0, 0, 0, 0, 0, //81-90
69  0,235};
70 const G4int G4NeutronInelasticXS::amax[] = {0,
71  0, 0, 7, 0,11,13,15,18, 0, 0, //1-10
72  0, 0, 0,30, 0, 0, 0,40, 0,48, //11-20
73  0, 0, 0, 0, 0,58, 0,64,65,70, //21-30
74  0,76, 0, 0, 0, 0, 0, 0, 0,96, //31-40
75  0, 0, 0, 0, 0, 0,109,116, 0,124, //41-50
76  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //51-60
77  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //61-70
78  0, 0, 0,186, 0, 0, 0, 0, 0, 0, //71-80
79  0,208, 0, 0, 0, 0, 0, 0, 0, 0, //81-90
80  0,238};
81 
83  : G4VCrossSectionDataSet("G4NeutronInelasticXS"),
85 {
86  // verboseLevel = 0;
87  if(verboseLevel > 0){
88  G4cout << "G4NeutronInelasticXS::G4NeutronInelasticXS Initialise for Z < "
89  << MAXZINEL << G4endl;
90  }
91  data.SetName("NeutronInelastic");
92  work.resize(13,0);
93  temp.resize(13,0.0);
94  coeff.resize(MAXZINEL, 1.0);
95  ggXsection = new G4GlauberGribovCrossSection();
96  fNucleon = new G4HadronNucleonXsc();
97  isInitialized = false;
98 }
99 
101 {
102  delete fNucleon;
103 }
104 
106 {
107  outFile << "G4NeutronInelasticXS calculates the neutron inelastic scattering\n"
108  << "cross section on nuclei using data from the high precision\n"
109  << "neutron database. These data are simplified and smoothed over\n"
110  << "the resonance region in order to reduce CPU time.\n"
111  << "G4NeutronInelasticXS is valid for energies up to 20 MeV, for\n"
112  << "nuclei through U.\n";
113 }
114 
115 G4bool
117  G4int, const G4Material*)
118 {
119  return true;
120 }
121 
122 G4bool
124  G4int /*ZZ*/, G4int /*AA*/,
125  const G4Element*, const G4Material*)
126 {
127  return true;
128 }
129 
130 G4double
132  G4int Z, const G4Material*)
133 {
134  G4double xs = 0.0;
135  G4double ekin = aParticle->GetKineticEnergy();
136 
137  if(Z < 1 || Z >= MAXZINEL) { return xs; }
138  G4int Amean = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
139 
140  G4PhysicsVector* pv = data.GetElementData(Z);
141  // G4cout << "G4NeutronInelasticXS::GetCrossSection e= " << ekin << " Z= "
142  // << Z << G4endl;
143 
144  // element was not initialised
145  if(!pv) {
146  Initialise(Z);
147  pv = data.GetElementData(Z);
148  if(!pv) { return xs; }
149  }
150 
151  G4double e1 = pv->Energy(0);
152  if(ekin <= e1) { return xs; }
153 
154  G4double e2 = pv->GetMaxEnergy();
155 
156  if(ekin <= e2) {
157  xs = pv->Value(ekin);
158  } else if(1 == Z) {
159  fNucleon->GetHadronNucleonXscPDG(aParticle, proton);
160  xs = coeff[1]*fNucleon->GetInelasticHadronNucleonXsc();
161  } else {
162  ggXsection->GetIsoCrossSection(aParticle, Z, Amean);
163  xs = coeff[Z]*ggXsection->GetInelasticGlauberGribovXsc();
164  }
165 
166  if(verboseLevel > 0) {
167  G4cout << "ekin= " << ekin << ", XSinel= " << xs << G4endl;
168  }
169  return xs;
170 }
171 
172 G4double
174  G4int Z, G4int A,
175  const G4Isotope*, const G4Element*,
176  const G4Material*)
177 {
178  G4double xs = 0.0;
179  G4double ekin = aParticle->GetKineticEnergy();
180  if(Z > 0 && Z < MAXZINEL) { xs = IsoCrossSection(ekin, Z, A); }
181  return xs;
182 }
183 
184 G4double
185 G4NeutronInelasticXS::IsoCrossSection(G4double ekin, G4int Z, G4int A)
186 {
187  G4double xs = 0.0;
188 
189  G4PhysicsVector* pv = data.GetElementData(Z);
190 
191  // element was not initialised
192  if(!pv) {
193  Initialise(Z);
194  pv = data.GetElementData(Z);
195  if(!pv) { return xs; }
196  }
197  G4PhysicsVector* pviso = data.GetComponentDataByID(Z, A);
198  if(pviso) { pv = pviso; }
199 
200  xs = pv->Value(ekin);
201 
202  if(verboseLevel > 0){
203  G4cout << "ekin= " << ekin << ", xs= " << xs << G4endl;
204  }
205  return xs;
206 }
207 
209  G4double kinEnergy)
210 {
211  G4int nIso = anElement->GetNumberOfIsotopes();
212  G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
213  G4Isotope* iso = (*isoVector)[0];
214 
215  // more than 1 isotope
216  if(1 < nIso) {
217  G4int Z = G4lrint(anElement->GetZ());
218  if(Z >= MAXZINEL) { Z = MAXZINEL - 1; }
219  G4double* abundVector = anElement->GetRelativeAbundanceVector();
220  G4double q = G4UniformRand();
221  G4double sum = 0.0;
222 
223  // is there isotope wise cross section?
224  if(0 == amin[Z]) {
225  for (G4int j = 0; j<nIso; ++j) {
226  sum += abundVector[j];
227  if(q <= sum) {
228  iso = (*isoVector)[j];
229  break;
230  }
231  }
232  } else {
233  size_t nmax = data.GetNumberOfComponents(Z);
234  if(temp.size() < nmax) { temp.resize(nmax,0.0); }
235  for (size_t i=0; i<nmax; ++i) {
236  G4int A = (*isoVector)[i]->GetN();
237  sum += abundVector[i]*IsoCrossSection(kinEnergy, Z, A);
238  temp[i] = sum;
239  }
240  sum *= q;
241  for (size_t j = 0; j<nmax; ++j) {
242  if(temp[j] >= sum) {
243  iso = (*isoVector)[j];
244  break;
245  }
246  }
247  }
248  }
249  return iso;
250 }
251 
252 void
254 {
255  if(isInitialized) { return; }
256  if(verboseLevel > 0){
257  G4cout << "G4NeutronCaptureXS::BuildPhysicsTable for "
258  << p.GetParticleName() << G4endl;
259  }
260  if(p.GetParticleName() != "neutron") {
262  ed << p.GetParticleName() << " is a wrong particle type -"
263  << " only neutron is allowed";
264  G4Exception("G4NeutronInelasticXS::BuildPhysicsTable(..)","had012",
265  FatalException, ed, "");
266  return;
267  }
268  isInitialized = true;
269 
270  // check environment variable
271  // Build the complete string identifying the file with the data set
272  char* path = getenv("G4NEUTRONXSDATA");
273 
274  G4DynamicParticle* dynParticle =
276 
277  // Access to elements
278  const G4ElementTable* theElmTable = G4Element::GetElementTable();
279  size_t numOfElm = G4Element::GetNumberOfElements();
280  if(numOfElm > 0) {
281  for(size_t i=0; i<numOfElm; ++i) {
282  G4int Z = G4lrint(((*theElmTable)[i])->GetZ());
283  if(Z < 1) { Z = 1; }
284  else if(Z >= MAXZINEL) { Z = MAXZINEL-1; }
285  //G4cout << "Z= " << Z << G4endl;
286  // Initialisation
287  if(!data.GetElementData(Z)) { Initialise(Z, dynParticle, path); }
288  }
289  }
290  delete dynParticle;
291 }
292 
293 void
294 G4NeutronInelasticXS::Initialise(G4int Z, G4DynamicParticle* dp,
295  const char* p)
296 {
297  if(data.GetElementData(Z)) { return; }
298  const char* path = p;
299  if(!p) {
300  // check environment variable
301  // Build the complete string identifying the file with the data set
302  path = getenv("G4NEUTRONXSDATA");
303  if (!path) {
304  G4Exception("G4NeutronInelasticXS::Initialise(..)","had013",
306  "Environment variable G4NEUTRONXSDATA is not defined");
307  return;
308  }
309  }
310  G4DynamicParticle* dynParticle = dp;
311  if(!dp) {
312  dynParticle =
314  }
315 
316  G4int Amean = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
317 
318  // upload element data
319  std::ostringstream ost;
320  ost << path << "/inelast" << Z ;
321  G4PhysicsVector* v = RetrieveVector(ost, true);
322  data.InitialiseForElement(Z, v);
323 
324  // upload isotope data
325  if(amin[Z] > 0) {
326  size_t n = 0;
327  size_t i = 0;
328  size_t nmax = (size_t)(amax[Z]-amin[Z]+1);
329  if(work.size() < nmax) { work.resize(nmax,0); }
330  for(G4int A=amin[Z]; A<=amax[Z]; ++A) {
331  std::ostringstream ost1;
332  ost1 << path << "/cap" << Z << "_" << A;
333  G4PhysicsVector* v1 = RetrieveVector(ost1, false);
334  if(v1) { ++n; }
335  work[i] = v1;
336  ++i;
337  }
338  data.InitialiseForComponent(Z, n);
339  for(size_t j=0; j<i; ++j) {
340  if(work[j]) { data.AddComponent(Z, amin[Z]+j, work[j]); }
341  }
342  }
343 
344  // smooth transition
345  G4double emax = v->GetMaxEnergy();
346  G4double sig1 = (*v)[v->GetVectorLength() - 1];
347  dynParticle->SetKineticEnergy(emax);
348  G4double sig2 = 0.0;
349  if(1 == Z) {
350  fNucleon->GetHadronNucleonXscPDG(dynParticle, proton);
351  sig2 = fNucleon->GetInelasticHadronNucleonXsc();
352  } else {
353  ggXsection->GetIsoCrossSection(dynParticle, Z, Amean);
354  sig2 = ggXsection->GetInelasticGlauberGribovXsc();
355  }
356  if(sig2 > 0.) { coeff[Z] = sig1/sig2; }
357  if(!dp) { delete dynParticle; }
358 }
359 
361 G4NeutronInelasticXS::RetrieveVector(std::ostringstream& ost, G4bool warn)
362 {
363  G4PhysicsLogVector* v = 0;
364  std::ifstream filein(ost.str().c_str());
365  if (!(filein)) {
366  if(!warn) { return v; }
368  ed << "Data file <" << ost.str().c_str()
369  << "> is not opened!";
370  G4Exception("G4NeutronElasticXS::RetrieveVector(..)","had014",
371  FatalException, ed, "Check G4NEUTRONXSDATA");
372  }else{
373  if(verboseLevel > 1) {
374  G4cout << "File " << ost.str()
375  << " is opened by G4NeutronCaptureXS" << G4endl;
376  }
377  // retrieve data from DB
378  v = new G4PhysicsLogVector();
379  if(!v->Retrieve(filein, true)) {
381  ed << "Data file <" << ost.str().c_str()
382  << "> is not retrieved!";
383  G4Exception("G4NeutronInelasticXS::RetrieveVector(..)","had015",
384  FatalException, ed, "Check G4NEUTRONXSDATA");
385  }
386  }
387  return v;
388 }
size_t GetNumberOfIsotopes() const
Definition: G4Element.hh:158
virtual G4bool Retrieve(std::ifstream &fIn, G4bool ascii)
G4double GetMaxEnergy() const
std::vector< G4Isotope * > G4IsotopeVector
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)
G4double GetZ() const
Definition: G4Element.hh:131
const G4int MAXZINEL
virtual G4bool IsIsoApplicable(const G4DynamicParticle *, G4int Z, G4int A, const G4Element *, const G4Material *)
size_t GetVectorLength() const
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
const G4String & GetParticleName() const
virtual G4double GetElementCrossSection(const G4DynamicParticle *, G4int Z, const G4Material *mat=0)
Char_t n[5]
#define G4UniformRand()
Definition: Randomize.hh:87
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
G4double * GetRelativeAbundanceVector() const
Definition: G4Element.hh:166
const G4int nmax
virtual G4double GetIsoCrossSection(const G4DynamicParticle *, G4int Z, G4int A, const G4Isotope *iso, const G4Element *elm, const G4Material *mat)
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)
tuple v
Definition: test.py:18
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4IsotopeVector * GetIsotopeVector() const
Definition: G4Element.hh:162
int G4lrint(double ad)
Definition: templates.hh:163
G4int nIso
virtual void CrossSectionDescription(std::ostream &) const
virtual G4bool IsElementApplicable(const G4DynamicParticle *, G4int Z, const G4Material *)
virtual G4Isotope * SelectIsotope(const G4Element *, G4double kinEnergy)
G4double GetHadronNucleonXscPDG(const G4DynamicParticle *, const G4ParticleDefinition *)
#define G4endl
Definition: G4ios.hh:61
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
std::vector< G4Element * > G4ElementTable
double G4double
Definition: G4Types.hh:76
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:395
const XML_Char const XML_Char * data
Definition: expat.h:268
G4double GetInelasticHadronNucleonXsc()