Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 93682 2015-10-28 10:09:49Z 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 // factory
58 #include "G4CrossSectionFactory.hh"
59 //
61 
62 using namespace std;
63 
64 const G4int G4NeutronInelasticXS::amin[] = {
65  0,
66  0, 0, 6, 0,10,12,14,16, 0, 0, //1-10
67  0, 0, 0,28, 0, 0, 0,36, 0,40, //11-20
68  0, 0, 0, 0, 0,54, 0,58,63,64, //21-30
69  0,70, 0, 0, 0, 0, 0, 0, 0,90, //31-40
70  0, 0, 0, 0, 0, 0,107,106, 0,112, //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,180, 0, 0, 0, 0, 0, 0, //71-80
74  0,204, 0, 0, 0, 0, 0, 0, 0, 0, //81-90
75  0,235};
76 const G4int G4NeutronInelasticXS::amax[] = {
77  0,
78  0, 0, 7, 0,11,13,15,18, 0, 0, //1-10
79  0, 0, 0,30, 0, 0, 0,40, 0,48, //11-20
80  0, 0, 0, 0, 0,58, 0,64,65,70, //21-30
81  0,76, 0, 0, 0, 0, 0, 0, 0,96, //31-40
82  0, 0, 0, 0, 0, 0,109,116, 0,124, //41-50
83  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //51-60
84  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //61-70
85  0, 0, 0,186, 0, 0, 0, 0, 0, 0, //71-80
86  0,208, 0, 0, 0, 0, 0, 0, 0, 0, //81-90
87  0,238};
88 
89 G4double G4NeutronInelasticXS::coeff[] = {1.0};
90 
91 G4ElementData* G4NeutronInelasticXS::data = 0;
92 
94  : G4VCrossSectionDataSet(Default_Name()),
96 {
97  // verboseLevel = 0;
98  if(verboseLevel > 0){
99  G4cout << "G4NeutronInelasticXS::G4NeutronInelasticXS Initialise for Z < "
100  << MAXZINEL << G4endl;
101  }
102  ggXsection = new G4ComponentGGHadronNucleusXsc();
103  fNucleon = new G4HadronNucleonXsc();
104  isMaster = false;
105 }
106 
108 {
109  //G4cout << "G4NeutronInelasticXS::~G4NeutronInelasticXS() "
110  // << " isMaster= " << isMaster << " data: " << data << G4endl;
111  delete fNucleon;
112  delete ggXsection;
113  if(isMaster) { delete data; data = 0; }
114 }
115 
116 void G4NeutronInelasticXS::CrossSectionDescription(std::ostream& outFile) const
117 {
118  outFile << "G4NeutronInelasticXS calculates the neutron inelastic scattering\n"
119  << "cross section on nuclei using data from the high precision\n"
120  << "neutron database. These data are simplified and smoothed over\n"
121  << "the resonance region in order to reduce CPU time.\n"
122  << "G4NeutronInelasticXS is valid for energies up to 20 MeV, for\n"
123  << "nuclei through U.\n";
124 }
125 
126 G4bool
128  G4int, const G4Material*)
129 {
130  return true;
131 }
132 
133 G4bool
135  G4int /*ZZ*/, G4int /*AA*/,
136  const G4Element*, const G4Material*)
137 {
138  return true;
139 }
140 
142  const G4DynamicParticle* aParticle,
143  G4int Z, const G4Material*)
144 {
145  G4double xs = 0.0;
146  G4double ekin = aParticle->GetKineticEnergy();
147 
148  if(Z >= MAXZINEL) { Z = MAXZINEL - 1; }
149  else if(Z < 1) { Z = 1; }
150 
151  G4int Amean = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
152 
153  G4PhysicsVector* pv = data->GetElementData(Z);
154  // G4cout << "G4NeutronInelasticXS::GetCrossSection e= " << ekin
155  // << " Z= " << Z << G4endl;
156 
157  // element was not initialised
158  if(!pv) {
159  Initialise(Z);
160  pv = data->GetElementData(Z);
161  if(!pv) { return xs; }
162  }
163 
164  G4double e1 = pv->Energy(0);
165  if(ekin <= e1) { return xs; }
166 
167  G4double e2 = pv->GetMaxEnergy();
168 
169  if(ekin <= e2) {
170  xs = pv->Value(ekin);
171  } else if(1 == Z) {
172  fNucleon->GetHadronNucleonXscPDG(aParticle, proton);
173  xs = coeff[1]*fNucleon->GetInelasticHadronNucleonXsc();
174  } else {
175  ggXsection->GetIsoCrossSection(aParticle, Z, Amean);
176  xs = coeff[Z]*ggXsection->GetInelasticGlauberGribovXsc();
177  }
178 
179  if(verboseLevel > 0) {
180  G4cout << "ekin= " << ekin << ", XSinel= " << xs << G4endl;
181  }
182  return xs;
183 }
184 
186  const G4DynamicParticle* aParticle,
187  G4int Z, G4int A,
188  const G4Isotope*, const G4Element*,
189  const G4Material*)
190 {
191  if(Z >= MAXZINEL) { Z = MAXZINEL - 1; }
192  else if(Z < 1) { Z = 1; }
193 
194  return IsoCrossSection(aParticle->GetKineticEnergy(), Z, A);
195 }
196 
197 G4double
198 G4NeutronInelasticXS::IsoCrossSection(G4double ekin, G4int Z, G4int A)
199 {
200  G4double xs = 0.0;
201  /*
202  G4cout << "IsoCrossSection Z= " << Z << " A= " << A
203  << " Amin= " << amin[Z] << " Amax= " << amax[Z]
204  << " E(MeV)= " << ekin << G4endl;
205  */
206  // element was not initialised
207  G4PhysicsVector* pv = data->GetElementData(Z);
208  if(!pv) {
209  Initialise(Z);
210  pv = data->GetElementData(Z);
211  }
212 
213  // isotope cross section exist
214  if(pv && amin[Z] > 0 && A >= amin[Z] && A <= amax[Z]) {
215  pv = data->GetComponentDataByIndex(Z, A - amin[Z]);
216  if(pv && ekin > pv->Energy(0)) { xs = pv->Value(ekin); }
217  }
218  if(verboseLevel > 0){
219  G4cout << "ekin= " << ekin << ", xs= " << xs << G4endl;
220  }
221  return xs;
222 }
223 
225  const G4Element* anElement, G4double kinEnergy)
226 {
227  size_t nIso = anElement->GetNumberOfIsotopes();
228  G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
229  G4Isotope* iso = (*isoVector)[0];
230 
231  //G4cout << "SelectIsotope NIso= " << nIso << G4endl;
232 
233  // more than 1 isotope
234  if(1 < nIso) {
235  G4int Z = G4lrint(anElement->GetZ());
236  //G4cout << "SelectIsotope Z= " << Z << G4endl;
237 
238  G4double* abundVector = anElement->GetRelativeAbundanceVector();
239  G4double q = G4UniformRand();
240  G4double sum = 0.0;
241 
242  // is there isotope wise cross section?
243  size_t j;
244  if(0 == amin[Z] || Z >= MAXZINEL) {
245  for (j = 0; j<nIso; ++j) {
246  sum += abundVector[j];
247  if(q <= sum) {
248  iso = (*isoVector)[j];
249  break;
250  }
251  }
252  } else {
253 
254  // element may be not initialised in unit test
255  if(!data->GetElementData(Z)) { Initialise(Z); }
256  size_t nn = temp.size();
257  if(nn < nIso) { temp.resize(nIso, 0.); }
258 
259  for (j=0; j<nIso; ++j) {
260  //G4cout << j << "-th isotope " << (*isoVector)[j]->GetN()
261  // << " abund= " << abundVector[j] << G4endl;
262  sum += abundVector[j]*IsoCrossSection(kinEnergy, Z,
263  (*isoVector)[j]->GetN());
264  temp[j] = sum;
265  }
266  sum *= q;
267  for (j = 0; j<nIso; ++j) {
268  if(temp[j] >= sum) {
269  iso = (*isoVector)[j];
270  break;
271  }
272  }
273  }
274  }
275  return iso;
276 }
277 
278 void
280 {
281  if(verboseLevel > 0){
282  G4cout << "G4NeutronInelasticXS::BuildPhysicsTable for "
283  << p.GetParticleName() << G4endl;
284  }
285  if(p.GetParticleName() != "neutron") {
287  ed << p.GetParticleName() << " is a wrong particle type -"
288  << " only neutron is allowed";
289  G4Exception("G4NeutronInelasticXS::BuildPhysicsTable(..)","had012",
290  FatalException, ed, "");
291  return;
292  }
293 
294  if(!data) {
295  isMaster = true;
296  data = new G4ElementData();
297  data->SetName("NeutronInelastic");
298  temp.resize(13,0.0);
299  }
300 
301  // it is possible re-initialisation for the new run
302  if(isMaster) {
303 
304  // check environment variable
305  // Build the complete string identifying the file with the data set
306  char* path = getenv("G4NEUTRONXSDATA");
307 
308  G4DynamicParticle* dynParticle =
310 
311  // Access to elements
312  const G4ElementTable* theElmTable = G4Element::GetElementTable();
313  size_t numOfElm = G4Element::GetNumberOfElements();
314  if(numOfElm > 0) {
315  for(size_t i=0; i<numOfElm; ++i) {
316  G4int Z = G4lrint(((*theElmTable)[i])->GetZ());
317  if(Z < 1) { Z = 1; }
318  else if(Z >= MAXZINEL) { Z = MAXZINEL-1; }
319  //G4cout << "Z= " << Z << G4endl;
320  // Initialisation
321  if(!(data->GetElementData(Z))) {
322  Initialise(Z, dynParticle, path);
323  }
324  }
325  }
326  delete dynParticle;
327  }
328 }
329 
330 void
331 G4NeutronInelasticXS::Initialise(G4int Z, G4DynamicParticle* dp,
332  const char* p)
333 {
334  if(data->GetElementData(Z) || Z < 1 || Z >= MAXZINEL) { return; }
335  const char* path = p;
336  if(!p) {
337  // check environment variable
338  // Build the complete string identifying the file with the data set
339  path = getenv("G4NEUTRONXSDATA");
340  if (!path) {
341  G4Exception("G4NeutronInelasticXS::Initialise(..)","had013",
343  "Environment variable G4NEUTRONXSDATA is not defined");
344  return;
345  }
346  }
347  G4DynamicParticle* dynParticle = dp;
348  if(!dp) {
349  dynParticle =
351  }
352 
353  G4int Amean = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
354 
355  // upload element data
356  std::ostringstream ost;
357  ost << path << "/inelast" << Z ;
358  G4PhysicsVector* v = RetrieveVector(ost, true);
359  data->InitialiseForElement(Z, v);
360  /*
361  G4cout << "G4NeutronInelasticXS::Initialise for Z= " << Z
362  << " A= " << Amean << " Amin= " << amin[Z]
363  << " Amax= " << amax[Z] << G4endl;
364  */
365  // upload isotope data
366  if(amin[Z] > 0) {
367  size_t nmax = (size_t)(amax[Z]-amin[Z]+1);
368  data->InitialiseForComponent(Z, nmax);
369 
370  for(G4int A=amin[Z]; A<=amax[Z]; ++A) {
371  std::ostringstream ost1;
372  ost1 << path << "/inelast" << Z << "_" << A;
373  G4PhysicsVector* v1 = RetrieveVector(ost1, false);
374  data->AddComponent(Z, A, v1);
375  }
376  }
377 
378  // smooth transition
379  G4double emax = v->GetMaxEnergy();
380  G4double sig1 = (*v)[v->GetVectorLength() - 1];
381  dynParticle->SetKineticEnergy(emax);
382  G4double sig2 = 0.0;
383  if(1 == Z) {
384  fNucleon->GetHadronNucleonXscPDG(dynParticle, proton);
385  sig2 = fNucleon->GetInelasticHadronNucleonXsc();
386  } else {
387  ggXsection->GetIsoCrossSection(dynParticle, Z, Amean);
388  sig2 = ggXsection->GetInelasticGlauberGribovXsc();
389  }
390  if(sig2 > 0.) { coeff[Z] = sig1/sig2; }
391  if(!dp) { delete dynParticle; }
392 }
393 
395 G4NeutronInelasticXS::RetrieveVector(std::ostringstream& ost, G4bool warn)
396 {
397  G4PhysicsLogVector* v = 0;
398  std::ifstream filein(ost.str().c_str());
399  if (!(filein)) {
400  if(warn) {
402  ed << "Data file <" << ost.str().c_str()
403  << "> is not opened!";
404  G4Exception("G4NeutronInelasticXS::RetrieveVector(..)","had014",
405  FatalException, ed, "Check G4NEUTRONXSDATA");
406  }
407  } else {
408  if(verboseLevel > 1) {
409  G4cout << "File " << ost.str()
410  << " is opened by G4NeutronInelasticXS" << G4endl;
411  }
412  // retrieve data from DB
413  v = new G4PhysicsLogVector();
414  if(!v->Retrieve(filein, true)) {
416  ed << "Data file <" << ost.str().c_str()
417  << "> is not retrieved!";
418  G4Exception("G4NeutronInelasticXS::RetrieveVector(..)","had015",
419  FatalException, ed, "Check G4NEUTRONXSDATA");
420  }
421  }
422  return v;
423 }
size_t GetNumberOfIsotopes() const
Definition: G4Element.hh:159
G4double GetMaxEnergy() const
std::vector< G4Isotope * > G4IsotopeVector
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
const char * p
Definition: xmltok.h:285
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
G4double GetIsoCrossSection(const G4DynamicParticle *, G4int Z, G4int A, const G4Isotope *iso=0, const G4Element *elm=0, const G4Material *mat=0)
const XML_Char const XML_Char * data
Definition: expat.h:268
virtual G4double GetElementCrossSection(const G4DynamicParticle *, G4int Z, const G4Material *mat=0)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
static size_t GetNumberOfElements()
Definition: G4Element.cc:405
bool G4bool
Definition: G4Types.hh:79
virtual G4bool Retrieve(std::ifstream &fIn, G4bool ascii) final
G4double * GetRelativeAbundanceVector() const
Definition: G4Element.hh:167
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
static const G4double emax
G4IsotopeVector * GetIsotopeVector() const
Definition: G4Element.hh:163
#define G4_DECLARE_XS_FACTORY(cross_section)
int G4lrint(double ad)
Definition: templates.hh:163
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:398
G4double GetInelasticHadronNucleonXsc()