Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4NeutronCaptureXS.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: G4NeutronCaptureXS.cc 110787 2018-06-14 06:43:31Z 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 // factory
53 #include "G4CrossSectionFactory.hh"
54 //
56 
57 using namespace std;
58 
60  0,
61  1, 4, 6, 9, 10, 12, 14, 16, 19, 20, //1-10
62  23, 24, 27, 28, 31, 32, 35, 36, 39, 40, //11-20
63  45, 46, 50, 50, 55, 54, 59, 58, 63, 64, //21-30
64  69, 70, 75, 0, 0, 0, 0, 0, 0, 90, //31-40
65  0, 92, 0, 0, 0, 102, 107, 106, 113, 112, //41-50
66  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //51-60
67  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //61-70
68  0, 0, 181, 180, 0, 0, 0, 192, 197, 0, //71-80
69  0, 204, 209, 0, 0, 0, 0, 0, 0, 0, //81-90
70  0, 235};
72  0,
73  1, 4, 7, 9, 11, 13, 15, 18, 19, 22, //1-10
74  23, 26, 27, 30, 31, 34, 37, 40, 41, 48, //11-20
75  45, 50, 51, 54, 55, 58, 59, 64, 65, 70, //21-30
76  71, 76, 75, 0, 0, 0, 0, 0, 0, 96, //31-40
77  0, 100, 0, 0, 0, 110, 109, 116, 115, 124, //41-50
78  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //51-60
79  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //61-70
80  0, 0, 181, 186, 0, 0, 0, 198, 197, 0, //71-80
81  0, 208, 209, 0, 0, 0, 0, 0, 0, 0, //81-90
82  0, 238};
83 
85 
86 #ifdef G4MULTITHREADED
87  G4Mutex G4NeutronCaptureXS::neutronCaptureXSMutex = G4MUTEX_INITIALIZER;
88 #endif
89 
91  : G4VCrossSectionDataSet(Default_Name()),
92  emax(20*CLHEP::MeV),elimit(1.0e-10*CLHEP::eV)
93 {
94  // verboseLevel = 0;
95  if(verboseLevel > 0){
96  G4cout << "G4NeutronCaptureXS::G4NeutronCaptureXS: Initialise for Z < "
97  << MAXZCAPTURE << G4endl;
98  }
99  isMaster = false;
100 }
101 
103 {
104  if(isMaster) { delete data; data = nullptr; }
105 }
106 
107 void G4NeutronCaptureXS::CrossSectionDescription(std::ostream& outFile) const
108 {
109  outFile << "G4NeutronCaptureXS calculates the neutron capture cross sections\n"
110  << "on nuclei using data from the high precision neutron database.\n"
111  << "These data are simplified and smoothed over the resonance region\n"
112  << "in order to reduce CPU time. G4NeutronCaptureXS is set to zero\n"
113  << "above 20 MeV for all targets. Cross section is zero also for Z>92.\n";
114 }
115 
116 G4bool
118  G4int, const G4Material*)
119 {
120  return true;
121 }
122 
123 G4bool
125  G4int, G4int,
126  const G4Element*, const G4Material*)
127 {
128  return true;
129 }
130 
131 G4double
133  G4int Z, const G4Material*)
134 {
135  G4double xs = 0.0;
136  G4double ekin = aParticle->GetKineticEnergy();
137  if(ekin > emax || Z < 1 || Z >= MAXZCAPTURE) { return xs; }
138  if(ekin < elimit) { ekin = elimit; }
139 
140  // element was not initialised
141  G4PhysicsVector* pv = data->GetElementData(Z);
142  if(!pv) { return xs; }
143 
144  G4double e1 = pv->Energy(0);
145  if(ekin < e1) { xs = (*pv)[0]*std::sqrt(e1/ekin); }
146  else if(ekin <= pv->GetMaxEnergy()) { xs = pv->Value(ekin); }
147 
148  if(verboseLevel > 0){
149  G4cout << "ekin= " << ekin << ", xs= " << xs << G4endl;
150  }
151  return xs;
152 }
153 
154 G4double
156  G4int Z, G4int A,
157  const G4Isotope*, const G4Element*,
158  const G4Material*)
159 {
160  return IsoCrossSection(aParticle->GetKineticEnergy(), Z, A);
161 }
162 
164 {
165  G4double xs = 0.0;
166  if(ekin > emax || Z < 1 || Z >= MAXZCAPTURE) { return xs; }
167  if(ekin < elimit) { ekin = elimit; }
168 
169  G4PhysicsVector* pviso = data->GetComponentDataByID(Z, A - amin[Z]);
170  if(pviso) {
171  G4double e1 = pviso->Energy(1);
172  if(ekin < e1) { xs = (*pviso)[1]*std::sqrt(e1/ekin); }
173  else if(ekin <= pviso->GetMaxEnergy()) { xs = pviso->Value(ekin); }
174  if(verboseLevel > 0) {
175  G4cout << "G4NeutronCaptureXS::IsoXS: Ekin(MeV)= " << ekin/MeV
176  << " xs(b)= " << xs/barn
177  << " Z= " << Z << " A= " << A << G4endl;
178  }
179  return xs;
180  }
181  // isotope data are not available or applicable
182  G4PhysicsVector* pv = data->GetElementData(Z);
183  if(pv) {
184  G4double e1 = pv->Energy(1);
185  if(ekin < e1) { xs = (*pv)[1]*std::sqrt(e1/ekin); }
186  else if(ekin <= pv->GetMaxEnergy()) { xs = pv->Value(ekin); }
187  if(verboseLevel > 0) {
188  G4cout << "G4NeutronCaptureXS::IsoXS: Ekin(MeV)= " << ekin/MeV
189  << " xs(b)= " << xs/barn
190  << " Z= " << Z << " A= " << A << G4endl;
191  }
192  }
193  return xs;
194 }
195 
196 const G4Isotope*
198  G4double kinEnergy)
199 {
200  size_t nIso = anElement->GetNumberOfIsotopes();
201  const G4Isotope* iso = anElement->GetIsotope(0);
202 
203  // more than 1 isotope
204  if(1 < nIso) {
205  G4int Z = anElement->GetZasInt();
206 
207  const G4double* abundVector = anElement->GetRelativeAbundanceVector();
208  G4double q = G4UniformRand();
209  G4double sum = 0.0;
210 
211  // is there isotope wise cross section?
212  size_t j;
213  if(0 == amin[Z] || Z >= MAXZCAPTURE) {
214  for (j = 0; j<nIso; ++j) {
215  sum += abundVector[j];
216  if(q <= sum) {
217  iso = anElement->GetIsotope(j);
218  break;
219  }
220  }
221  } else {
222  size_t nn = temp.size();
223  if(nn < nIso) { temp.resize(nIso, 0.); }
224 
225  for (j=0; j<nIso; ++j) {
226  sum += abundVector[j]*IsoCrossSection(kinEnergy, Z,
227  anElement->GetIsotope(j)->GetN());
228  temp[j] = sum;
229  }
230  sum *= q;
231  for (j = 0; j<nIso; ++j) {
232  if(temp[j] >= sum) {
233  iso = anElement->GetIsotope(j);
234  break;
235  }
236  }
237  }
238  }
239  return iso;
240 }
241 
242 void
244 {
245  if(verboseLevel > 0){
246  G4cout << "G4NeutronCaptureXS::BuildPhysicsTable for "
247  << p.GetParticleName() << G4endl;
248  }
249  if(p.GetParticleName() != "neutron") {
251  ed << p.GetParticleName() << " is a wrong particle type -"
252  << " only neutron is allowed";
253  G4Exception("G4NeutronCaptureXS::BuildPhysicsTable(..)","had012",
254  FatalException, ed, "");
255  return;
256  }
257 
258  if(!data) {
259 #ifdef G4MULTITHREADED
260  G4MUTEXLOCK(&neutronCaptureXSMutex);
261  if(!data) {
262 #endif
263  isMaster = true;
264  data = new G4ElementData();
265  data->SetName("NeutronCapture");
266  temp.resize(13,0.0);
267 #ifdef G4MULTITHREADED
268  }
269  G4MUTEXUNLOCK(&neutronCaptureXSMutex);
270 #endif
271  }
272 
273  // it is possible re-initialisation for the second run
274  if(isMaster) {
275 
276  // check environment variable
277  // Build the complete string identifying the file with the data set
278  char* path = getenv("G4NEUTRONXSDATA");
279 
280  // Access to elements
281  const G4ElementTable* theElmTable = G4Element::GetElementTable();
282  size_t numOfElm = G4Element::GetNumberOfElements();
283  for(size_t i=0; i<numOfElm; ++i) {
284  G4int Z = ((*theElmTable)[i])->GetZasInt();
285  if(Z >= MAXZCAPTURE) { Z = MAXZCAPTURE-1; }
286  //G4cout << "Z= " << Z << G4endl;
287  // Initialisation
288  if(!data->GetElementData(Z)) { Initialise(Z, path); }
289  }
290  }
291 }
292 
293 void
295 {
296  if(data->GetElementData(Z) || Z < 1 || Z >= MAXZCAPTURE) { return; }
297  const char* path = p;
298 
299  // check environment variable
300  if(!p) {
301  path = getenv("G4NEUTRONXSDATA");
302  if (!path) {
303  G4Exception("G4NeutronCaptureXS::Initialise(..)","had013",FatalException,
304  "Environment variable G4NEUTRONXSDATA is not defined");
305  return;
306  }
307  }
308 
309  // upload element data
310  std::ostringstream ost;
311  ost << path << "/neutron/cap" << Z ;
312  G4PhysicsVector* v = RetrieveVector(ost, true);
313  data->InitialiseForElement(Z, v);
314 
315  // upload isotope data
316  if(amin[Z] > 0) {
317  size_t nmax = (size_t)(amax[Z]-amin[Z]+1);
318  data->InitialiseForComponent(Z, nmax);
319 
320  for(G4int A=amin[Z]; A<=amax[Z]; ++A) {
321  std::ostringstream ost1;
322  ost1 << path << "/neutron/cap" << Z << "_" << A;
323  v = RetrieveVector(ost1, false);
324  data->AddComponent(Z, A, v);
325  }
326  }
327 }
328 
330 G4NeutronCaptureXS::RetrieveVector(std::ostringstream& ost, G4bool warn)
331 {
332  G4PhysicsLogVector* v = nullptr;
333  std::ifstream filein(ost.str().c_str());
334  if (!(filein)) {
335  if(warn) {
337  ed << "Data file <" << ost.str().c_str()
338  << "> is not opened!";
339  G4Exception("G4NeutronCaptureXS::RetrieveVector(..)","had014",
340  FatalException, ed, "Check G4NEUTRONXSDATA");
341  }
342  } else {
343  if(verboseLevel > 1) {
344  G4cout << "File " << ost.str()
345  << " is opened by G4NeutronCaptureXS" << G4endl;
346  }
347  // retrieve data from DB
348  v = new G4PhysicsLogVector();
349  if(!v->Retrieve(filein, true)) {
351  ed << "Data file <" << ost.str().c_str()
352  << "> is not retrieved!";
353  G4Exception("G4NeutronCaptureXS::RetrieveVector(..)","had015",
354  FatalException, ed, "Check G4NEUTRONXSDATA");
355  }
356  }
357  return v;
358 }
G4double Energy(size_t index) const
virtual G4bool IsElementApplicable(const G4DynamicParticle *, G4int Z, const G4Material *)
std::vector< G4double > temp
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:45
G4double * GetRelativeAbundanceVector() const
Definition: G4Element.hh:167
static constexpr double MeV
Definition: G4SIunits.hh:214
static const G4int amin[MAXZCAPTURE]
#define G4endl
Definition: G4ios.hh:61
const char * p
Definition: xmltok.h:285
virtual G4double GetElementCrossSection(const G4DynamicParticle *, G4int Z, const G4Material *mat=nullptr)
static const G4int amax[MAXZCAPTURE]
const G4String & GetParticleName() const
static const G4double emax
virtual void CrossSectionDescription(std::ostream &) const
G4double Value(G4double theEnergy, size_t &lastidx) const
const G4int nmax
const XML_Char const XML_Char * data
Definition: expat.h:268
virtual G4double GetIsoCrossSection(const G4DynamicParticle *, G4int Z, G4int A, const G4Isotope *iso, const G4Element *elm, const G4Material *mat)
Float_t Z
virtual G4bool Retrieve(std::ifstream &fIn, G4bool ascii) final
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:88
const G4int MAXZCAPTURE
virtual const G4Isotope * SelectIsotope(const G4Element *, G4double kinEnergy)
std::vector< G4Element * > G4ElementTable
#define G4UniformRand()
Definition: Randomize.hh:53
double A(double temperature)
static constexpr double eV
Definition: G4SIunits.hh:215
static G4ElementData * data
#define G4_DECLARE_XS_FACTORY(cross_section)
#define G4MUTEXUNLOCK(mutex)
Definition: G4Threading.hh:234
G4int GetZasInt() const
Definition: G4Element.hh:132
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 G4int
Definition: G4Types.hh:78
void Initialise(G4int Z, const char *)
static constexpr double barn
Definition: G4SIunits.hh:105
#define G4MUTEXLOCK(mutex)
Definition: G4Threading.hh:233
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
const G4Isotope * GetIsotope(G4int iso) const
Definition: G4Element.hh:170
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
static size_t GetNumberOfElements()
Definition: G4Element.cc:405
virtual G4bool IsIsoApplicable(const G4DynamicParticle *, G4int Z, G4int A, const G4Element *, const G4Material *)
G4PhysicsVector * RetrieveVector(std::ostringstream &in, G4bool warn)
G4double IsoCrossSection(G4double ekin, G4int Z, G4int A)
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:398
size_t GetNumberOfIsotopes() const
Definition: G4Element.hh:159
std::mutex G4Mutex
Definition: G4Threading.hh:84