Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4ParticleInelasticXS.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: G4ParticleInelasticXS.cc 109964 2018-05-11 17:48:23Z vnivanch $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4ParticleInelasticXS
34 //
35 // Author Ivantchenko, Geant4, 3-Aug-09
36 //
37 // Modifications:
38 //
39 
40 #include "G4ParticleInelasticXS.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"
49 #include "G4HadronNucleonXsc.hh"
50 #include "G4NistManager.hh"
51 #include "G4Proton.hh"
52 #include "Randomize.hh"
53 
54 #include <iostream>
55 #include <fstream>
56 #include <sstream>
57 
58 using namespace std;
59 
61  0,
62  1, 4, 6, 9, 10, 12, 14, 16, 19, 20, //1-10
63  23, 24, 27, 28, 31, 32, 35, 36, 39, 40, //11-20
64  45, 46, 50, 50, 55, 54, 59, 58, 63, 64, //21-30
65  69, 70, 75, 0, 0, 0, 0, 0, 0, 90, //31-40
66  0, 92, 0, 0, 0, 102, 107, 106, 113, 112, //41-50
67  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //51-60
68  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //61-70
69  0, 0, 181, 180, 0, 0, 0, 192, 197, 0, //71-80
70  0, 204, 209, 0, 0, 0, 0, 0, 0, 0, //81-90
71  0, 235};
73  0,
74  1, 4, 7, 9, 11, 13, 15, 18, 19, 22, //1-10
75  23, 26, 27, 30, 31, 34, 37, 40, 41, 48, //11-20
76  45, 50, 51, 54, 55, 58, 59, 64, 65, 70, //21-30
77  71, 76, 75, 0, 0, 0, 0, 0, 0, 96, //31-40
78  0, 100, 0, 0, 0, 110, 109, 116, 115, 124, //41-50
79  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //51-60
80  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //61-70
81  0, 0, 181, 186, 0, 0, 0, 198, 197, 0, //71-80
82  0, 208, 209, 0, 0, 0, 0, 0, 0, 0, //81-90
83  0, 238};
84 
86 
88 
89 #ifdef G4MULTITHREADED
90  G4Mutex G4ParticleInelasticXS::particleInelasticXSMutex = G4MUTEX_INITIALIZER;
91 #endif
92 
94  : G4VCrossSectionDataSet("G4ParticleInelasticXS"),
95  ggXsection(nullptr),
96  nnXsection(nullptr),
97  fNucleon(nullptr),
98  particle(part),
100  particleName(""),
101  isMaster(false),
102  emax(20*CLHEP::MeV)
103 {
104  if(!part) {
105  G4Exception("G4ParticleInelasticXS::G4ParticleInelasticXS(..)","had015",
106  FatalException, "NO particle definition in constructor");
107  } else {
108  // verboseLevel = 0;
110  if(verboseLevel > 0){
111  G4cout << "G4ParticleInelasticXS::G4ParticleInelasticXS for "
112  << particleName << " on atoms with Z < " << MAXZINELP << G4endl;
113  }
114  if(particleName == "neutron" || particleName == "proton") {
117  } else {
119  }
120  }
122 }
123 
125 {
126  //G4cout << "G4ParticleInelasticXS::~G4ParticleInelasticXS() "
127  // << " isMaster= " << isMaster << " data: " << data << G4endl;
128  delete fNucleon;
129  delete ggXsection;
130  delete nnXsection;
131  if(isMaster) { delete data; data = nullptr; }
132 }
133 
134 void G4ParticleInelasticXS::CrossSectionDescription(std::ostream& outFile) const
135 {
136  outFile << "G4ParticleInelasticXS calculates n, p, d, t, he3, he4 inelastic \n"
137  << "cross section on nuclei using data from the high precision\n"
138  << "neutron database. These data are simplified and smoothed over\n"
139  << "the resonance region in order to reduce CPU time.\n"
140  << "For high energy Glauber-Gribov cross section model is used\n";
141 }
142 
143 G4bool
145  G4int, const G4Material*)
146 {
147  return true;
148 }
149 
150 G4bool
152  G4int, G4int,
153  const G4Element*, const G4Material*)
154 {
155  return true;
156 }
157 
159  const G4DynamicParticle* aParticle,
160  G4int ZZ, const G4Material*)
161 {
162  G4double xs = 0.0;
163  G4double ekin = aParticle->GetKineticEnergy();
164 
165  G4int Z = (ZZ >= MAXZINELP) ? MAXZINELP - 1 : ZZ;
166 
167  G4PhysicsVector* pv = data->GetElementData(Z);
168  // G4cout << "G4ParticleInelasticXS::GetCrossSection e= " << ekin
169  // << " Z= " << Z << G4endl;
170 
171  // element was not initialised or below threshold
172  if(!pv || ekin <= pv->Energy(0)) { return xs; }
173 
174  if(ekin <= pv->GetMaxEnergy()) {
175  xs = pv->Value(ekin);
176  } else if(1 == Z && fNucleon) {
179  } else if(ggXsection) {
180  G4int Amean =
181  G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
182  ggXsection->GetIsoCrossSection(aParticle, Z, Amean);
184  } else if(nnXsection) {
185  G4int Amean =
186  G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
187  nnXsection->GetZandACrossSection(aParticle, Z, Amean);
189  }
190 
191  if(verboseLevel > 0) {
192  G4cout << "ElmXS: Z= " << Z << " Ekin(MeV)= " << ekin/CLHEP::MeV
193  << " xs(bn)= " << xs/CLHEP::barn << " element data for "
194  << particleName << G4endl;
195  }
196  return xs;
197 }
198 
200  const G4DynamicParticle* aParticle,
201  G4int Z, G4int A,
202  const G4Isotope*, const G4Element*,
203  const G4Material*)
204 {
205  return IsoCrossSection(aParticle->GetKineticEnergy(), Z, A);
206 }
207 
208 G4double
210 {
211  G4double xs = 0.0;
212  G4int Z = (ZZ >= MAXZINELP) ? MAXZINELP - 1 : ZZ;
213  /*
214  G4cout << "IsoCrossSection Z= " << Z << " A= " << A
215  << " Amin= " << amin[Z] << " Amax= " << amax[Z]
216  << " E(MeV)= " << ekin << G4endl;
217  */
218  // first compute isotope cross section
219  if(ekin <=emax && amin[Z]>0 && A >= amin[Z] && A <= amax[Z]) {
220  G4PhysicsVector* pviso = data->GetComponentDataByIndex(Z, A - amin[Z]);
221  if(pviso) {
222  xs = pviso->Value(ekin);
223  if(verboseLevel > 0){
224  G4cout << "IsoXS for " << particleName
225  << " Target Z= " << Z << " A= " << A
226  << " Ekin(MeV)= " << ekin/CLHEP::MeV
227  << " xs(bn)= " << xs/CLHEP::barn << G4endl;
228  }
229  return xs;
230  }
231  }
232  // isotope data are not available or applicable
233  G4PhysicsVector* pv = data->GetElementData(Z);
234  if(pv) { xs = pv->Value(ekin); }
235  if(verboseLevel > 0) {
236  G4cout << "IsoXS for " << particleName
237  << " Target Z= " << Z << " A= " << A
238  << " Ekin(MeV)= " << ekin/CLHEP::MeV
239  << " xs(bn)= " << xs/CLHEP::barn << G4endl;
240  }
241  return xs;
242 }
243 
245  const G4Element* anElement, G4double kinEnergy)
246 {
247  size_t nIso = anElement->GetNumberOfIsotopes();
248  const G4Isotope* iso = anElement->GetIsotope(0);
249 
250  //G4cout << "SelectIsotope NIso= " << nIso << G4endl;
251  if(1 == nIso) { return iso; }
252 
253  // more than 1 isotope
254  G4int Z = anElement->GetZasInt();
255  //G4cout << "SelectIsotope Z= " << Z << G4endl;
256 
257  const G4double* abundVector = anElement->GetRelativeAbundanceVector();
258  G4double q = G4UniformRand();
259  G4double sum = 0.0;
260 
261  // is there isotope wise cross section?
262  size_t j;
263  if(kinEnergy > emax || 0 == amin[Z] || Z >= MAXZINELP) {
264  for (j = 0; j<nIso; ++j) {
265  sum += abundVector[j];
266  if(q <= sum) {
267  iso = anElement->GetIsotope(j);
268  break;
269  }
270  }
271  return iso;
272  }
273 
274  size_t nn = temp.size();
275  if(nn < nIso) { temp.resize(nIso, 0.); }
276 
277  for (j=0; j<nIso; ++j) {
278  //G4cout << j << "-th isotope " << (*isoVector)[j]->GetN()
279  // << " abund= " << abundVector[j] << G4endl;
280  sum += abundVector[j]*IsoCrossSection(kinEnergy, Z,
281  anElement->GetIsotope(j)->GetN());
282  temp[j] = sum;
283  }
284  sum *= q;
285  for (j = 0; j<nIso; ++j) {
286  if(temp[j] >= sum) {
287  iso = anElement->GetIsotope(j);
288  break;
289  }
290  }
291  return iso;
292 }
293 
294 void
296 {
297  if(verboseLevel > 0){
298  G4cout << "G4ParticleInelasticXS::BuildPhysicsTable for "
299  << p.GetParticleName() << G4endl;
300  }
301  if(&p != particle) {
303  ed << p.GetParticleName() << " is a wrong particle type -"
304  << particleName << " is expected";
305  G4Exception("G4ParticleInelasticXS::BuildPhysicsTable(..)","had012",
306  FatalException, ed, "");
307  return;
308  }
309 
310  if(!data) {
311 #ifdef G4MULTITHREADED
312  G4MUTEXLOCK(&particleInelasticXSMutex);
313  if(!data) {
314 #endif
315  isMaster = true;
316  data = new G4ElementData();
317  data->SetName(particleName + "Inelastic");
318  temp.resize(13,0.0);
319 #ifdef G4MULTITHREADED
320  }
321  G4MUTEXUNLOCK(&particleInelasticXSMutex);
322 #endif
323  }
324 
325  // it is possible re-initialisation for the new run
326  if(isMaster) {
327 
328  // check environment variable
329  // Build the complete string identifying the file with the data set
330  char* path = getenv("G4NEUTRONXSDATA");
331 
332  G4DynamicParticle* dynParticle =
334 
335  // Access to elements
336  const G4ElementTable* theElmTable = G4Element::GetElementTable();
337  size_t numOfElm = G4Element::GetNumberOfElements();
338  for(size_t i=0; i<numOfElm; ++i) {
339  G4int Z = ((*theElmTable)[i])->GetZasInt();
340  if(Z >= MAXZINELP) { Z = MAXZINELP-1; }
341  //G4cout << "Z= " << Z << G4endl;
342  // Initialisation
343  if(!(data->GetElementData(Z))) {
344  Initialise(Z, dynParticle, path);
345  }
346  }
347  delete dynParticle;
348  }
349 }
350 
352  const char* p)
353 {
354  if(data->GetElementData(Z)) { return; }
355  const char* path = p;
356  if(!p) {
357  // check environment variable
358  // Build the complete string identifying the file with the data set
359  path = getenv("G4NEUTRONXSDATA");
360  if (!path) {
361  G4Exception("G4ParticleInelasticXS::Initialise(..)","had013",
363  "Environment variable G4NEUTRONXSDATA is not defined");
364  return;
365  }
366  }
367 
368  // upload element data
369  std::ostringstream ost;
370  ost << path << "/" << particleName << "/inel" << Z ;
371  G4PhysicsVector* v = RetrieveVector(ost, true);
372  data->InitialiseForElement(Z, v);
373  /*
374  G4cout << "G4ParticleInelasticXS::Initialise for Z= " << Z
375  << " A= " << Amean << " Amin= " << amin[Z]
376  << " Amax= " << amax[Z] << G4endl;
377  */
378  // upload isotope data
379  if(amin[Z] > 0) {
380  size_t nmax = (size_t)(amax[Z]-amin[Z]+1);
381  data->InitialiseForComponent(Z, nmax);
382 
383  for(G4int A=amin[Z]; A<=amax[Z]; ++A) {
384  std::ostringstream ost1;
385  ost1 << path << "/neutron/inel" << Z << "_" << A;
386  G4PhysicsVector* v1 = RetrieveVector(ost1, false);
387  data->AddComponent(Z, A, v1);
388  }
389  }
390 
391  // smooth transition
392  G4double sig1 = (*v)[v->GetVectorLength()-1];
393  dp->SetKineticEnergy(v->GetMaxEnergy());
394  G4double sig2 = 0.0;
395  G4int Amean =
396  G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
397  if(1 == Z && fNucleon) {
400  } else if(ggXsection) {
401  ggXsection->GetIsoCrossSection(dp, Z, Amean);
403  } else if(nnXsection) {
404  nnXsection->GetZandACrossSection(dp, Z, Amean);
406  }
407  if(sig2 > 0.) { coeff[Z] = sig1/sig2; }
408 }
409 
411 G4ParticleInelasticXS::RetrieveVector(std::ostringstream& ost, G4bool warn)
412 {
413  G4PhysicsLogVector* v = nullptr;
414  std::ifstream filein(ost.str().c_str());
415  if (!(filein)) {
416  if(warn) {
418  ed << "Data file <" << ost.str().c_str()
419  << "> is not opened!";
420  G4Exception("G4ParticleInelasticXS::RetrieveVector(..)","had014",
421  FatalException, ed, "Check G4NEUTRONXSDATA");
422  }
423  } else {
424  if(verboseLevel > 1) {
425  G4cout << "File " << ost.str()
426  << " is opened by G4ParticleInelasticXS" << G4endl;
427  }
428  // retrieve data from DB
429  v = new G4PhysicsLogVector();
430  if(!v->Retrieve(filein, true)) {
432  ed << "Data file <" << ost.str().c_str()
433  << "> is not retrieved!";
434  G4Exception("G4ParticleInelasticXS::RetrieveVector(..)","had015",
435  FatalException, ed, "Check G4NEUTRONXSDATA");
436  }
437  }
438  return v;
439 }
G4double GetIsoCrossSection(const G4DynamicParticle *, G4int Z, G4int A, const G4Isotope *iso=0, const G4Element *elm=0, const G4Material *mat=0)
virtual G4double GetElementCrossSection(const G4DynamicParticle *, G4int Z, const G4Material *mat=nullptr)
G4double GetZandACrossSection(const G4DynamicParticle *, G4int Z, G4int A)
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
G4ComponentGGNuclNuclXsc * nnXsection
void SetKineticEnergy(G4double aEnergy)
static G4double coeff[MAXZINELP]
G4ComponentGGHadronNucleusXsc * ggXsection
#define G4endl
Definition: G4ios.hh:61
const char * p
Definition: xmltok.h:285
void Initialise(G4int Z, G4DynamicParticle *dp, const char *)
const G4String & GetParticleName() const
static const G4double emax
G4double Value(G4double theEnergy, size_t &lastidx) const
G4double IsoCrossSection(G4double ekin, G4int Z, G4int A)
const G4int nmax
const XML_Char const XML_Char * data
Definition: expat.h:268
static const G4int amin[MAXZINELP]
Float_t Z
virtual G4bool Retrieve(std::ifstream &fIn, G4bool ascii) final
G4PhysicsVector * RetrieveVector(std::ostringstream &in, G4bool warn)
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
TString part[npart]
Definition: Style.C:32
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:88
G4double GetHadronNucleonXscPDG(const G4DynamicParticle *, const G4ParticleDefinition *)
std::vector< G4Element * > G4ElementTable
static constexpr double MeV
#define G4UniformRand()
Definition: Randomize.hh:53
double A(double temperature)
G4double GetInelasticHadronNucleonXsc()
void SetForAllAtomsAndEnergies(G4bool val)
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
#define G4MUTEXUNLOCK(mutex)
Definition: G4Threading.hh:234
const G4ParticleDefinition * particle
G4int GetZasInt() const
Definition: G4Element.hh:132
virtual G4bool IsIsoApplicable(const G4DynamicParticle *, G4int Z, G4int A, const G4Element *, 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
G4ParticleInelasticXS(const G4ParticleDefinition *)
virtual const G4Isotope * SelectIsotope(const G4Element *, G4double kinEnergy)
virtual G4bool IsElementApplicable(const G4DynamicParticle *, G4int Z, const G4Material *)
#define G4MUTEXLOCK(mutex)
Definition: G4Threading.hh:233
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
virtual G4double GetIsoCrossSection(const G4DynamicParticle *, G4int Z, G4int A, const G4Isotope *iso, const G4Element *elm, const G4Material *mat)
G4HadronNucleonXsc * fNucleon
const G4Isotope * GetIsotope(G4int iso) const
Definition: G4Element.hh:170
const G4ParticleDefinition * proton
static size_t GetNumberOfElements()
Definition: G4Element.cc:405
virtual void CrossSectionDescription(std::ostream &) const
std::vector< G4double > temp
static G4ElementData * data
static const G4int amax[MAXZINELP]
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
const G4int MAXZINELP