ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NeutronCaptureXS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 // -------------------------------------------------------------------
27 //
28 // GEANT4 Class file
29 //
30 //
31 // File name: G4NeutronCaptureXS
32 //
33 // Author Ivantchenko, Geant4, 3-Aug-09
34 //
35 // Modifications:
36 //
37 
38 #include <fstream>
39 #include <sstream>
40 
41 #include "G4SystemOfUnits.hh"
42 #include "G4NeutronCaptureXS.hh"
43 #include "G4Material.hh"
44 #include "G4Element.hh"
45 #include "G4PhysicsLogVector.hh"
46 #include "G4PhysicsVector.hh"
47 #include "G4DynamicParticle.hh"
48 #include "G4ProductionCutsTable.hh"
49 #include "Randomize.hh"
50 #include "G4Log.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  2, 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 
86 
87 #ifdef G4MULTITHREADED
88  G4Mutex G4NeutronCaptureXS::neutronCaptureXSMutex = G4MUTEX_INITIALIZER;
89 #endif
90 
92  : G4VCrossSectionDataSet(Default_Name()),
93  emax(20*CLHEP::MeV),elimit(1.0e-10*CLHEP::eV)
94 {
95  // verboseLevel = 0;
96  if(verboseLevel > 0){
97  G4cout << "G4NeutronCaptureXS::G4NeutronCaptureXS: Initialise for Z < "
98  << MAXZCAPTURE << G4endl;
99  }
101  isMaster = false;
102  temp.resize(13,0.0);
103 }
104 
106 {
107  if(isMaster) { delete data; data = nullptr; }
108 }
109 
110 void G4NeutronCaptureXS::CrossSectionDescription(std::ostream& outFile) const
111 {
112  outFile << "G4NeutronCaptureXS calculates the neutron capture cross sections\n"
113  << "on nuclei using data from the high precision neutron database.\n"
114  << "These data are simplified and smoothed over the resonance region\n"
115  << "in order to reduce CPU time. G4NeutronCaptureXS is set to zero\n"
116  << "above 20 MeV for all targets. Cross section is zero also for Z>92.\n";
117 }
118 
119 G4bool
121  G4int, const G4Material*)
122 {
123  return true;
124 }
125 
126 G4bool
128  G4int, G4int,
129  const G4Element*, const G4Material*)
130 {
131  return true;
132 }
133 
134 G4double
136  G4int ZZ, const G4Material*)
137 {
138  G4double xs = 0.0;
139  G4double ekin = aParticle->GetKineticEnergy();
140  if(ekin > emax) { return xs; }
141 
142  G4int Z = std::min(ZZ, MAXZCAPTURE-1);
143  G4double logEkin = aParticle->GetLogKineticEnergy();
144  if(ekin < elimit) { ekin = elimit; logEkin = logElimit; }
145 
146  auto pv = GetPhysicsVector(Z);
147  if(!pv) { return xs; }
148 
149  G4double e1 = pv->Energy(0);
150  if(ekin < e1) {
151  xs = (*pv)[0]*std::sqrt(e1/ekin);
152  } else if(ekin <= pv->GetMaxEnergy()) {
153  xs = pv->LogVectorValue(ekin, logEkin);
154  }
155 
156  if(verboseLevel > 1){
157  G4cout << "Ekin= " << ekin/CLHEP::MeV
158  << " ElmXScap(b)= " << xs/CLHEP::barn << G4endl;
159  }
160  return xs;
161 }
162 
163 G4double
165  G4int Z, G4int A,
166  const G4Isotope*, const G4Element*,
167  const G4Material*)
168 {
169  return IsoCrossSection(aParticle->GetKineticEnergy(),
170  aParticle->GetLogKineticEnergy(),
171  Z, A);
172 }
173 
175  G4int ZZ, G4int A)
176 {
177  G4double xs = 0.0;
178  if(eKin > emax) { return xs; }
179 
180  G4int Z = std::min(ZZ, MAXZCAPTURE-1);
181  G4double ekin = eKin;
182  G4double logEkin = logE;
183  if(ekin < elimit) {
184  ekin = elimit;
185  logEkin = logElimit;
186  }
187 
188  auto pv = GetPhysicsVector(Z);
189  if(!pv) { return xs; }
190 
191  if(amin[Z] > 0 && A >= amin[Z] && A <= amax[Z]) {
192  G4PhysicsVector* pviso = data->GetComponentDataByID(Z, A - amin[Z]);
193  if(pviso) {
194  G4double e1 = pviso->Energy(1);
195  if(ekin < e1) {
196  xs = (*pviso)[1]*std::sqrt(e1/ekin);
197  } else if(ekin <= pviso->GetMaxEnergy()) {
198  xs = pviso->LogVectorValue(ekin, logEkin);
199  }
200  if(verboseLevel > 0) {
201  G4cout << "G4NeutronCaptureXS::IsoXS: Ekin(MeV)= " << ekin/MeV
202  << " xs(b)= " << xs/barn
203  << " Z= " << Z << " A= " << A << G4endl;
204  }
205  return xs;
206  }
207  }
208  // isotope data are not available or applicable
209  G4double e1 = pv->Energy(1);
210  if(ekin < e1) {
211  xs = (*pv)[1]*std::sqrt(e1/ekin);
212  } else if(ekin <= pv->GetMaxEnergy()) {
213  xs = pv->LogVectorValue(ekin, logEkin);
214  }
215  if(verboseLevel > 0) {
216  G4cout << "G4NeutronCaptureXS::IsoXS: Ekin(MeV)= " << ekin/MeV
217  << " xs(b)= " << xs/barn
218  << " Z= " << Z << " A= " << A << " no iso XS" << G4endl;
219  }
220  return xs;
221 }
222 
223 const G4Isotope*
225  G4double kinEnergy, G4double logE)
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 
236  const G4double* abundVector = anElement->GetRelativeAbundanceVector();
237  G4double q = G4UniformRand();
238  G4double sum = 0.0;
239 
240  // is there isotope wise cross section?
241  size_t j;
242  if(0 == amin[Z] || Z >= MAXZCAPTURE) {
243  for (j = 0; j<nIso; ++j) {
244  sum += abundVector[j];
245  if(q <= sum) {
246  iso = anElement->GetIsotope(j);
247  break;
248  }
249  }
250  return iso;
251  }
252  size_t nn = temp.size();
253  if(nn < nIso) { temp.resize(nIso, 0.); }
254 
255  for (j=0; j<nIso; ++j) {
256  sum += abundVector[j]*IsoCrossSection(kinEnergy, logE, Z,
257  anElement->GetIsotope(j)->GetN());
258  temp[j] = sum;
259  }
260  sum *= q;
261  for (j = 0; j<nIso; ++j) {
262  if(temp[j] >= sum) {
263  iso = anElement->GetIsotope(j);
264  break;
265  }
266  }
267  return iso;
268 }
269 
270 void
272 {
273  if(verboseLevel > 0){
274  G4cout << "G4NeutronCaptureXS::BuildPhysicsTable for "
275  << p.GetParticleName() << G4endl;
276  }
277  if(p.GetParticleName() != "neutron") {
279  ed << p.GetParticleName() << " is a wrong particle type -"
280  << " only neutron is allowed";
281  G4Exception("G4NeutronCaptureXS::BuildPhysicsTable(..)","had012",
282  FatalException, ed, "");
283  return;
284  }
285 
286  if(!data) {
287 #ifdef G4MULTITHREADED
288  G4MUTEXLOCK(&neutronCaptureXSMutex);
289  if(!data) {
290 #endif
291  isMaster = true;
292  data = new G4ElementData();
293  data->SetName("NeutronCapture");
295 #ifdef G4MULTITHREADED
296  }
297  G4MUTEXUNLOCK(&neutronCaptureXSMutex);
298 #endif
299  }
300 
301  // it is possible re-initialisation for the second run
302  if(isMaster) {
303 
304  // Access to elements
305  auto theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
306  size_t numOfCouples = theCoupleTable->GetTableSize();
307  for(size_t j=0; j<numOfCouples; ++j) {
308  auto mat = theCoupleTable->GetMaterialCutsCouple(j)->GetMaterial();
309  auto elmVec = mat->GetElementVector();
310  size_t numOfElem = mat->GetNumberOfElements();
311  for (size_t ie = 0; ie < numOfElem; ++ie) {
312  G4int Z = std::max(1,std::min(((*elmVec)[ie])->GetZasInt(), MAXZCAPTURE-1));
313  if(!data->GetElementData(Z)) { Initialise(Z); }
314  }
315  }
316  }
317 }
318 
320 {
321  const G4PhysicsVector* pv = data->GetElementData(Z);
322  if(!pv) {
323  InitialiseOnFly(Z);
324  pv = data->GetElementData(Z);
325  }
326  return pv;
327 }
328 
330 {
331  // check environment variable
332  // build the complete string identifying the file with the data set
333  if(gDataDirectory.empty()) {
334  char* path = std::getenv("G4PARTICLEXSDATA");
335  if (path) {
336  std::ostringstream ost;
337  ost << path << "/neutron/cap";
338  gDataDirectory = ost.str();
339  } else {
340  G4Exception("G4NeutronCaptureXS::Initialise(..)","had013",
342  "Environment variable G4PARTICLEXSDATA is not defined");
343  }
344  }
345  return gDataDirectory;
346 }
347 
349 {
350 #ifdef G4MULTITHREADED
351  G4MUTEXLOCK(&neutronCaptureXSMutex);
352  if(!data->GetElementData(Z)) {
353 #endif
354  Initialise(Z);
355 #ifdef G4MULTITHREADED
356  }
357  G4MUTEXUNLOCK(&neutronCaptureXSMutex);
358 #endif
359 }
360 
362 {
363  if(data->GetElementData(Z)) { return; }
364 
365  // upload element data
366  std::ostringstream ost;
367  ost << FindDirectoryPath() << Z ;
368  G4PhysicsVector* v = RetrieveVector(ost, true);
369  data->InitialiseForElement(Z, v);
370 
371  // upload isotope data
372  if(amin[Z] > 0) {
373  size_t nmax = (size_t)(amax[Z]-amin[Z]+1);
374  data->InitialiseForComponent(Z, nmax);
375 
376  for(G4int A=amin[Z]; A<=amax[Z]; ++A) {
377  std::ostringstream ost1;
378  ost1 << gDataDirectory << Z << "_" << A;
379  v = RetrieveVector(ost1, false);
380  data->AddComponent(Z, A, v);
381  }
382  }
383 }
384 
386 G4NeutronCaptureXS::RetrieveVector(std::ostringstream& ost, G4bool warn)
387 {
388  G4PhysicsLogVector* v = nullptr;
389  std::ifstream filein(ost.str().c_str());
390  if (!(filein)) {
391  if(warn) {
393  ed << "Data file <" << ost.str().c_str()
394  << "> is not opened!";
395  G4Exception("G4NeutronCaptureXS::RetrieveVector(..)","had014",
396  FatalException, ed, "Check G4PARTICLEXSDATA");
397  }
398  } else {
399  if(verboseLevel > 1) {
400  G4cout << "File " << ost.str()
401  << " is opened by G4NeutronCaptureXS" << G4endl;
402  }
403  // retrieve data from DB
404  v = new G4PhysicsLogVector();
405  if(!v->Retrieve(filein, true)) {
407  ed << "Data file <" << ost.str().c_str()
408  << "> is not retrieved!";
409  G4Exception("G4NeutronCaptureXS::RetrieveVector(..)","had015",
410  FatalException, ed, "Check G4PARTICLEXSDATA");
411  }
412  }
413  return v;
414 }