Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4NeutronInelasticXS.cc
이 파일의 문서화 페이지로 가기
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 110787 2018-06-14 06:43:31Z 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  1, 4, 6, 9, 10, 12, 14, 16, 19, 20, //1-10
67  23, 24, 27, 28, 31, 32, 35, 36, 39, 40, //11-20
68  45, 46, 50, 50, 55, 54, 59, 58, 63, 64, //21-30
69  69, 70, 75, 0, 0, 0, 0, 0, 0, 90, //31-40
70  0, 92, 0, 0, 0, 102, 107, 106, 113, 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, 181, 180, 0, 0, 0, 192, 197, 0, //71-80
74  0, 204, 209, 0, 0, 0, 0, 0, 0, 0, //81-90
75  0, 235};
77  0,
78  1, 4, 7, 9, 11, 13, 15, 18, 19, 22, //1-10
79  23, 26, 27, 30, 31, 34, 37, 40, 41, 48, //11-20
80  45, 50, 51, 54, 55, 58, 59, 64, 65, 70, //21-30
81  71, 76, 75, 0, 0, 0, 0, 0, 0, 96, //31-40
82  0, 100, 0, 0, 0, 110, 109, 116, 115, 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, 181, 186, 0, 0, 0, 198, 197, 0, //71-80
86  0, 208, 209, 0, 0, 0, 0, 0, 0, 0, //81-90
87  0, 238};
88 
90 
92 
93 #ifdef G4MULTITHREADED
94  G4Mutex G4NeutronInelasticXS::neutronInelasticXSMutex = G4MUTEX_INITIALIZER;
95 #endif
96 
98  : G4VCrossSectionDataSet(Default_Name()),
99  proton(G4Proton::Proton()), emax(20*CLHEP::MeV)
100 {
101  // verboseLevel = 0;
102  if(verboseLevel > 0){
103  G4cout << "G4NeutronInelasticXS::G4NeutronInelasticXS Initialise for Z < "
104  << MAXZINEL << G4endl;
105  }
108  isMaster = false;
109 }
110 
112 {
113  //G4cout << "G4NeutronInelasticXS::~G4NeutronInelasticXS() "
114  // << " isMaster= " << isMaster << " data: " << data << G4endl;
115  delete fNucleon;
116  delete ggXsection;
117  if(isMaster) { delete data; data = nullptr; }
118 }
119 
120 void G4NeutronInelasticXS::CrossSectionDescription(std::ostream& outFile) const
121 {
122  outFile << "G4NeutronInelasticXS calculates the neutron inelastic scattering\n"
123  << "cross section on nuclei using data from the high precision\n"
124  << "neutron database. These data are simplified and smoothed over\n"
125  << "the resonance region in order to reduce CPU time.\n"
126  << "For high energy Glauber-Gribov cross section model is used\n";
127 }
128 
129 G4bool
131  G4int, const G4Material*)
132 {
133  return true;
134 }
135 
136 G4bool
138  G4int, G4int,
139  const G4Element*, const G4Material*)
140 {
141  return true;
142 }
143 
145  const G4DynamicParticle* aParticle,
146  G4int ZZ, const G4Material*)
147 {
148  G4double xs = 0.0;
149  G4double ekin = aParticle->GetKineticEnergy();
150 
151  G4int Z = (ZZ >= MAXZINEL) ? MAXZINEL - 1 : ZZ;
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 || ekin <= pv->Energy(0)) { return xs; }
159 
160  if(ekin <= pv->GetMaxEnergy()) {
161  xs = pv->Value(ekin);
162  } else if(1 == Z) {
165  } else {
166  G4int Amean =
167  G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
168  ggXsection->GetIsoCrossSection(aParticle, Z, Amean);
170  }
171 
172  if(verboseLevel > 0) {
173  G4cout << "Z= " << Z << " Ekin(MeV)= " << ekin/CLHEP::MeV
174  << ", nElmXSinel(bn)= " << xs/CLHEP::barn
175  << G4endl;
176  }
177  return xs;
178 }
179 
181  const G4DynamicParticle* aParticle,
182  G4int Z, G4int A,
183  const G4Isotope*, const G4Element*,
184  const G4Material*)
185 {
186  return IsoCrossSection(aParticle->GetKineticEnergy(), Z, A);
187 }
188 
189 G4double
191 {
192  G4double xs = 0.0;
193  G4int Z = (ZZ >= MAXZINEL) ? MAXZINEL - 1 : ZZ;
194 
195  /*
196  G4cout << "IsoCrossSection Z= " << Z << " A= " << A
197  << " Amin= " << amin[Z] << " Amax= " << amax[Z]
198  << " E(MeV)= " << ekin << G4endl;
199  */
200  // first compute isotope cross section
201  if(ekin <=emax && amin[Z]>0 && A >= amin[Z] && A <= amax[Z]) {
202  G4PhysicsVector* pviso = data->GetComponentDataByIndex(Z, A - amin[Z]);
203  if(pviso) {
204  xs = pviso->Value(ekin);
205  if(verboseLevel > 0) {
206  G4cout << "IsoXS: Z= " << Z << " A= " << A
207  << " Ekin(MeV)= " << ekin/CLHEP::MeV
208  << ", nElmXSinel(bn)= " << xs/CLHEP::barn << G4endl;
209  }
210  return xs;
211  }
212  }
213  // isotope data are not available or applicable
214  G4PhysicsVector* pv = data->GetElementData(Z);
215  if(pv) { xs = pv->Value(ekin); }
216  if(verboseLevel > 0) {
217  G4cout << "IsoXS: Z= " << Z << " A= " << A
218  << " Ekin(MeV)= " << ekin/CLHEP::MeV
219  << ", nElmXSinel(bn)= " << xs/CLHEP::barn << G4endl;
220  }
221  return xs;
222 }
223 
225  const G4Element* anElement, G4double kinEnergy)
226 {
227  size_t nIso = anElement->GetNumberOfIsotopes();
228  const G4Isotope* iso = anElement->GetIsotope(0);
229 
230  //G4cout << "SelectIsotope NIso= " << nIso << G4endl;
231  if(1 == nIso) { return iso; }
232 
233  // more than 1 isotope
234  G4int Z = anElement->GetZasInt();
235  //G4cout << "SelectIsotope Z= " << Z << G4endl;
236 
237  const G4double* abundVector = anElement->GetRelativeAbundanceVector();
238  G4double q = G4UniformRand();
239  G4double sum = 0.0;
240 
241  // isotope wise cross section not available
242  size_t j;
243  if(kinEnergy > emax || 0 >= amin[Z] || Z >= MAXZINEL) {
244  for (j = 0; j<nIso; ++j) {
245  sum += abundVector[j];
246  if(q <= sum) {
247  iso = anElement->GetIsotope(j);
248  break;
249  }
250  }
251  return iso;
252  }
253 
254  // use isotope cross sections
255  size_t nn = temp.size();
256  if(nn < nIso) { temp.resize(nIso, 0.); }
257 
258  for (j=0; j<nIso; ++j) {
259  //G4cout << j << "-th isotope " << (*isoVector)[j]->GetN()
260  // << " abund= " << abundVector[j] << G4endl;
261  sum += abundVector[j]*IsoCrossSection(kinEnergy, Z,
262  anElement->GetIsotope(j)->GetN());
263  temp[j] = sum;
264  }
265  sum *= q;
266  for (j = 0; j<nIso; ++j) {
267  if(temp[j] >= sum) {
268  iso = anElement->GetIsotope(j);
269  break;
270  }
271  }
272  return iso;
273 }
274 
275 void
277 {
278  if(verboseLevel > 0){
279  G4cout << "G4NeutronInelasticXS::BuildPhysicsTable for "
280  << p.GetParticleName() << G4endl;
281  }
282  if(p.GetParticleName() != "neutron") {
284  ed << p.GetParticleName() << " is a wrong particle type -"
285  << " only neutron is allowed";
286  G4Exception("G4NeutronInelasticXS::BuildPhysicsTable(..)","had012",
287  FatalException, ed, "");
288  return;
289  }
290 
291  if(!data) {
292 #ifdef G4MULTITHREADED
293  G4MUTEXLOCK(&neutronInelasticXSMutex);
294  if(!data) {
295 #endif
296  isMaster = true;
297  data = new G4ElementData();
298  data->SetName("NeutronInelastic");
299  temp.resize(13,0.0);
300 #ifdef G4MULTITHREADED
301  }
302  G4MUTEXUNLOCK(&neutronInelasticXSMutex);
303 #endif
304  }
305 
306  // it is possible re-initialisation for the new run
307  if(isMaster) {
308 
309  // check environment variable
310  // Build the complete string identifying the file with the data set
311  char* path = getenv("G4NEUTRONXSDATA");
312 
313  G4DynamicParticle* dynParticle =
315 
316  // Access to elements
317  const G4ElementTable* theElmTable = G4Element::GetElementTable();
318  size_t numOfElm = G4Element::GetNumberOfElements();
319  for(size_t i=0; i<numOfElm; ++i) {
320  G4int Z = ((*theElmTable)[i])->GetZasInt();
321  if(Z >= MAXZINEL) { Z = MAXZINEL-1; }
322  //G4cout << "Z= " << Z << G4endl;
323  // Initialisation
324  if(!(data->GetElementData(Z))) {
325  Initialise(Z, dynParticle, path);
326  }
327  }
328  delete dynParticle;
329  }
330 }
331 
332 void
334  const char* p)
335 {
336  if(data->GetElementData(Z) || Z < 1 || Z >= MAXZINEL) { return; }
337  const char* path = p;
338  if(!p) {
339  // check environment variable
340  // Build the complete string identifying the file with the data set
341  path = getenv("G4NEUTRONXSDATA");
342  if (!path) {
343  G4Exception("G4NeutronInelasticXS::Initialise(..)","had013",
345  "Environment variable G4NEUTRONXSDATA is not defined");
346  return;
347  }
348  }
349  // upload element data
350  std::ostringstream ost;
351  ost << path << "/neutron/inel" << Z ;
352  G4PhysicsVector* v = RetrieveVector(ost, true);
353  data->InitialiseForElement(Z, v);
354  /*
355  G4cout << "G4NeutronInelasticXS::Initialise for Z= " << Z
356  << " A= " << Amean << " Amin= " << amin[Z]
357  << " Amax= " << amax[Z] << G4endl;
358  */
359  // upload isotope data
360  if(amin[Z] > 0) {
361  size_t nmax = (size_t)(amax[Z]-amin[Z]+1);
362  data->InitialiseForComponent(Z, nmax);
363 
364  for(G4int A=amin[Z]; A<=amax[Z]; ++A) {
365  std::ostringstream ost1;
366  ost1 << path << "/neutron/inel" << Z << "_" << A;
367  G4PhysicsVector* v1 = RetrieveVector(ost1, false);
368  data->AddComponent(Z, A, v1);
369  }
370  }
371 
372  // smooth transition
373  G4double sig1 = (*v)[v->GetVectorLength()-1];
374  dp->SetKineticEnergy(v->GetMaxEnergy());
375  G4double sig2 = 0.0;
376  if(1 == Z) {
379  } else {
380  G4int Amean =
381  G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
382  ggXsection->GetIsoCrossSection(dp, Z, Amean);
384  }
385  if(sig2 > 0.) { coeff[Z] = sig1/sig2; }
386 }
387 
389 G4NeutronInelasticXS::RetrieveVector(std::ostringstream& ost, G4bool warn)
390 {
391  G4PhysicsLogVector* v = nullptr;
392  std::ifstream filein(ost.str().c_str());
393  if (!(filein)) {
394  if(warn) {
396  ed << "Data file <" << ost.str().c_str()
397  << "> is not opened!";
398  G4Exception("G4NeutronInelasticXS::RetrieveVector(..)","had014",
399  FatalException, ed, "Check G4NEUTRONXSDATA");
400  }
401  } else {
402  if(verboseLevel > 1) {
403  G4cout << "File " << ost.str()
404  << " is opened by G4NeutronInelasticXS" << G4endl;
405  }
406  // retrieve data from DB
407  v = new G4PhysicsLogVector();
408  if(!v->Retrieve(filein, true)) {
410  ed << "Data file <" << ost.str().c_str()
411  << "> is not retrieved!";
412  G4Exception("G4NeutronInelasticXS::RetrieveVector(..)","had015",
413  FatalException, ed, "Check G4NEUTRONXSDATA");
414  }
415  }
416  return v;
417 }
static G4double coeff[MAXZINEL]
G4HadronNucleonXsc * fNucleon
G4double GetIsoCrossSection(const G4DynamicParticle *, G4int Z, G4int A, const G4Isotope *iso=0, const G4Element *elm=0, const G4Material *mat=0)
CLHEP::Hep3Vector G4ThreeVector
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:45
G4double * GetRelativeAbundanceVector() const
Definition: G4Element.hh:167
static constexpr double MeV
Definition: G4SIunits.hh:214
void SetKineticEnergy(G4double aEnergy)
#define G4endl
Definition: G4ios.hh:61
const char * p
Definition: xmltok.h:285
static const G4int amax[MAXZINEL]
const G4String & GetParticleName() const
virtual G4double GetIsoCrossSection(const G4DynamicParticle *, G4int Z, G4int A, const G4Isotope *iso, const G4Element *elm, const G4Material *mat)
void Initialise(G4int Z, G4DynamicParticle *dp, const char *)
std::vector< G4double > temp
const G4int MAXZINEL
static const G4double emax
G4double Value(G4double theEnergy, size_t &lastidx) const
G4ComponentGGHadronNucleusXsc * ggXsection
virtual void CrossSectionDescription(std::ostream &) const
const G4int nmax
const XML_Char const XML_Char * data
Definition: expat.h:268
Float_t Z
virtual G4bool Retrieve(std::ifstream &fIn, G4bool ascii) final
static const G4int amin[MAXZINEL]
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:88
G4double GetHadronNucleonXscPDG(const G4DynamicParticle *, const G4ParticleDefinition *)
G4PhysicsVector * RetrieveVector(std::ostringstream &in, G4bool warn)
std::vector< G4Element * > G4ElementTable
static constexpr double MeV
#define G4UniformRand()
Definition: Randomize.hh:53
double A(double temperature)
G4double IsoCrossSection(G4double ekin, G4int Z, G4int A)
G4double GetInelasticHadronNucleonXsc()
virtual G4bool IsIsoApplicable(const G4DynamicParticle *, G4int Z, G4int A, const G4Element *, const G4Material *)
#define G4_DECLARE_XS_FACTORY(cross_section)
#define G4MUTEXUNLOCK(mutex)
Definition: G4Threading.hh:234
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
G4int GetZasInt() const
Definition: G4Element.hh:132
virtual G4bool IsElementApplicable(const G4DynamicParticle *, G4int Z, const G4Material *)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
G4int GetN() const
Definition: G4Isotope.hh:94
int G4lrint(double ad)
Definition: templates.hh:151
G4double GetMaxEnergy() const
int G4int
Definition: G4Types.hh:78
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
virtual G4double GetElementCrossSection(const G4DynamicParticle *, G4int Z, const G4Material *mat=nullptr)
static G4ElementData * data
#define G4MUTEXLOCK(mutex)
Definition: G4Threading.hh:233
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
virtual const G4Isotope * SelectIsotope(const G4Element *, G4double kinEnergy)
const G4Isotope * GetIsotope(G4int iso) const
Definition: G4Element.hh:170
static size_t GetNumberOfElements()
Definition: G4Element.cc:405
const G4ParticleDefinition * proton
size_t GetVectorLength() const
static constexpr double barn
Definition: SystemOfUnits.h:85
static G4NistManager * Instance()
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:398
size_t GetNumberOfIsotopes() const
Definition: G4Element.hh:159
std::mutex G4Mutex
Definition: G4Threading.hh:84