Geant4  10.01.p02
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 83697 2014-09-10 07:15:29Z 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 
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};
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 
90 
92 
94  : G4VCrossSectionDataSet(Default_Name()),
96 {
97  // verboseLevel = 0;
98  if(verboseLevel > 0){
99  G4cout << "G4NeutronInelasticXS::G4NeutronInelasticXS Initialise for Z < "
100  << MAXZINEL << G4endl;
101  }
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 < 1 || Z >= MAXZINEL) { return xs; }
149  G4int Amean = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
150 
152  // G4cout << "G4NeutronInelasticXS::GetCrossSection e= " << ekin
153  // << " Z= " << Z << G4endl;
154 
155  // element was not initialised
156  if(!pv) {
157  Initialise(Z);
158  pv = data->GetElementData(Z);
159  if(!pv) { return xs; }
160  }
161 
162  G4double e1 = pv->Energy(0);
163  if(ekin <= e1) { return xs; }
164 
165  G4double e2 = pv->GetMaxEnergy();
166 
167  if(ekin <= e2) {
168  xs = pv->Value(ekin);
169  } else if(1 == Z) {
172  } else {
173  ggXsection->GetIsoCrossSection(aParticle, Z, Amean);
175  }
176 
177  if(verboseLevel > 0) {
178  G4cout << "ekin= " << ekin << ", XSinel= " << xs << G4endl;
179  }
180  return xs;
181 }
182 
184  const G4DynamicParticle* aParticle,
185  G4int Z, G4int A,
186  const G4Isotope*, const G4Element*,
187  const G4Material*)
188 {
189  G4double xs = 0.0;
190  G4double ekin = aParticle->GetKineticEnergy();
191  if(Z > 0 && Z < MAXZINEL) { xs = IsoCrossSection(ekin, Z, A); }
192  return xs;
193 }
194 
195 G4double
197 {
198  G4double xs = 0.0;
199 
201 
202  // element was not initialised
203  if(!pv) {
204  Initialise(Z);
205  pv = data->GetElementData(Z);
206  if(!pv) { return xs; }
207  }
208  G4PhysicsVector* pviso = data->GetComponentDataByID(Z, A);
209  if(pviso) { pv = pviso; }
210 
211  xs = pv->Value(ekin);
212 
213  if(verboseLevel > 0){
214  G4cout << "ekin= " << ekin << ", xs= " << xs << G4endl;
215  }
216  return xs;
217 }
218 
220  const G4Element* anElement, G4double kinEnergy)
221 {
222  G4int nIso = anElement->GetNumberOfIsotopes();
223  G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
224  G4Isotope* iso = (*isoVector)[0];
225 
226  // more than 1 isotope
227  if(1 < nIso) {
228  G4int Z = G4lrint(anElement->GetZ());
229  if(Z >= MAXZINEL) { Z = MAXZINEL - 1; }
230  G4double* abundVector = anElement->GetRelativeAbundanceVector();
231  G4double q = G4UniformRand();
232  G4double sum = 0.0;
233 
234  // is there isotope wise cross section?
235  if(0 == amin[Z]) {
236  for (G4int j = 0; j<nIso; ++j) {
237  sum += abundVector[j];
238  if(q <= sum) {
239  iso = (*isoVector)[j];
240  break;
241  }
242  }
243  } else {
244  size_t nmax = data->GetNumberOfComponents(Z);
245  if(temp.size() < nmax) { temp.resize(nmax,0.0); }
246  for (size_t i=0; i<nmax; ++i) {
247  G4int A = (*isoVector)[i]->GetN();
248  sum += abundVector[i]*IsoCrossSection(kinEnergy, Z, A);
249  temp[i] = sum;
250  }
251  sum *= q;
252  for (size_t j = 0; j<nmax; ++j) {
253  if(temp[j] >= sum) {
254  iso = (*isoVector)[j];
255  break;
256  }
257  }
258  }
259  }
260  return iso;
261 }
262 
263 void
265 {
266  if(verboseLevel > 0){
267  G4cout << "G4NeutronInelasticXS::BuildPhysicsTable for "
268  << p.GetParticleName() << G4endl;
269  }
270  if(p.GetParticleName() != "neutron") {
272  ed << p.GetParticleName() << " is a wrong particle type -"
273  << " only neutron is allowed";
274  G4Exception("G4NeutronInelasticXS::BuildPhysicsTable(..)","had012",
275  FatalException, ed, "");
276  return;
277  }
278 
279  if(!data) {
280  isMaster = true;
281  data = new G4ElementData();
282  data->SetName("NeutronInelastic");
283  work.resize(13,0);
284  temp.resize(13,0.0);
285  }
286 
287  // it is possible re-initialisation for the new run
288  if(isMaster) {
289 
290  // check environment variable
291  // Build the complete string identifying the file with the data set
292  char* path = getenv("G4NEUTRONXSDATA");
293 
294  G4DynamicParticle* dynParticle =
296 
297  // Access to elements
298  const G4ElementTable* theElmTable = G4Element::GetElementTable();
299  size_t numOfElm = G4Element::GetNumberOfElements();
300  if(numOfElm > 0) {
301  for(size_t i=0; i<numOfElm; ++i) {
302  G4int Z = G4lrint(((*theElmTable)[i])->GetZ());
303  if(Z < 1) { Z = 1; }
304  else if(Z >= MAXZINEL) { Z = MAXZINEL-1; }
305  //G4cout << "Z= " << Z << G4endl;
306  // Initialisation
307  if(!(data->GetElementData(Z))) {
308  Initialise(Z, dynParticle, path);
309  }
310  }
311  }
312  delete dynParticle;
313  }
314 }
315 
316 void
318  const char* p)
319 {
320  if(data->GetElementData(Z)) { return; }
321  const char* path = p;
322  if(!p) {
323  // check environment variable
324  // Build the complete string identifying the file with the data set
325  path = getenv("G4NEUTRONXSDATA");
326  if (!path) {
327  G4Exception("G4NeutronInelasticXS::Initialise(..)","had013",
329  "Environment variable G4NEUTRONXSDATA is not defined");
330  return;
331  }
332  }
333  G4DynamicParticle* dynParticle = dp;
334  if(!dp) {
335  dynParticle =
337  }
338 
339  G4int Amean = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
340 
341  // upload element data
342  std::ostringstream ost;
343  ost << path << "/inelast" << Z ;
344  G4PhysicsVector* v = RetrieveVector(ost, true);
345  data->InitialiseForElement(Z, v);
346  /*
347  G4cout << "G4NeutronInelasticXS::Initialise for Z= " << Z
348  << " A= " << Amean << " Amin= " << amin[Z]
349  << " Amax= " << amax[Z] << G4endl;
350  */
351  // upload isotope data
352  if(amin[Z] > 0) {
353  size_t n = 0;
354  size_t i = 0;
355  size_t nmax = (size_t)(amax[Z]-amin[Z]+1);
356  if(work.size() < nmax) { work.resize(nmax,0); }
357  for(G4int A=amin[Z]; A<=amax[Z]; ++A) {
358  std::ostringstream ost1;
359  ost1 << path << "/inelast" << Z << "_" << A;
360  G4PhysicsVector* v1 = RetrieveVector(ost1, false);
361  if(v1) { ++n; }
362  work[i] = v1;
363  ++i;
364  }
365 
366  //G4cout << " n= " << n << G4endl;
367 
369  for(size_t j=0; j<i; ++j) {
370  if(work[j]) { data->AddComponent(Z, amin[Z]+j, work[j]); }
371  }
372  }
373 
374  // smooth transition
375  G4double emax = v->GetMaxEnergy();
376  G4double sig1 = (*v)[v->GetVectorLength() - 1];
377  dynParticle->SetKineticEnergy(emax);
378  G4double sig2 = 0.0;
379  if(1 == Z) {
380  fNucleon->GetHadronNucleonXscPDG(dynParticle, proton);
382  } else {
383  ggXsection->GetIsoCrossSection(dynParticle, Z, Amean);
385  }
386  if(sig2 > 0.) { coeff[Z] = sig1/sig2; }
387  if(!dp) { delete dynParticle; }
388 }
389 
391 G4NeutronInelasticXS::RetrieveVector(std::ostringstream& ost, G4bool warn)
392 {
393  G4PhysicsLogVector* v = 0;
394  std::ifstream filein(ost.str().c_str());
395  if (!(filein)) {
396  if(!warn) { return v; }
398  ed << "Data file <" << ost.str().c_str()
399  << "> is not opened!";
400  G4Exception("G4NeutronInelasticXS::RetrieveVector(..)","had014",
401  FatalException, ed, "Check G4NEUTRONXSDATA");
402  }else{
403  if(verboseLevel > 1) {
404  G4cout << "File " << ost.str()
405  << " is opened by G4NeutronInelasticXS" << G4endl;
406  }
407  // retrieve data from DB
408  v = new G4PhysicsLogVector();
409  if(!v->Retrieve(filein, true)) {
411  ed << "Data file <" << ost.str().c_str()
412  << "> is not retrieved!";
413  G4Exception("G4NeutronInelasticXS::RetrieveVector(..)","had015",
414  FatalException, ed, "Check G4NEUTRONXSDATA");
415  }
416  }
417  return v;
418 }
size_t GetNumberOfIsotopes() const
Definition: G4Element.hh:158
void Initialise(G4int Z, G4DynamicParticle *dp=0, const char *=0)
G4_DECLARE_XS_FACTORY(G4NeutronInelasticXS)
virtual G4bool Retrieve(std::ifstream &fIn, G4bool ascii)
G4double GetMaxEnergy() const
std::vector< G4Isotope * > G4IsotopeVector
void AddComponent(G4int Z, G4int id, G4PhysicsVector *v)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4PhysicsVector * GetComponentDataByID(G4int Z, G4int id)
G4PhysicsVector * GetElementData(G4int Z)
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
static G4double coeff[MAXZINEL]
static const G4double e2
const G4int MAXZINEL
static const G4int amax[MAXZINEL]
const G4ParticleDefinition * proton
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)
void InitialiseForElement(G4int Z, G4PhysicsVector *v)
#define G4UniformRand()
Definition: Randomize.hh:93
G4GLOB_DLL std::ostream G4cout
static size_t GetNumberOfElements()
Definition: G4Element.cc:410
std::vector< G4PhysicsVector * > work
bool G4bool
Definition: G4Types.hh:79
static const G4int amin[MAXZINEL]
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
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]
void SetKineticEnergy(G4double aEnergy)
static const G4double e1
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double IsoCrossSection(G4double ekin, G4int Z, G4int A)
G4IsotopeVector * GetIsotopeVector() const
Definition: G4Element.hh:162
int G4lrint(double ad)
Definition: templates.hh:163
std::vector< G4double > temp
virtual void CrossSectionDescription(std::ostream &) const
G4HadronNucleonXsc * fNucleon
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 &)
void SetName(const G4String &nam)
std::vector< G4Element * > G4ElementTable
double G4double
Definition: G4Types.hh:76
static G4ElementData * data
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:403
void InitialiseForComponent(G4int Z, G4int nComponents=0)
G4PhysicsVector * RetrieveVector(std::ostringstream &in, G4bool warn)
G4double GetInelasticHadronNucleonXsc()
G4GlauberGribovCrossSection * ggXsection