ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NeutronElasticXS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4NeutronElasticXS.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: G4NeutronElasticXS
32 //
33 // Author Ivantchenko, Geant4, 3-Aug-09
34 //
35 // Modifications:
36 //
37 
38 #include "G4NeutronElasticXS.hh"
39 #include "G4Neutron.hh"
40 #include "G4DynamicParticle.hh"
41 #include "G4ProductionCutsTable.hh"
42 #include "G4Material.hh"
43 #include "G4Element.hh"
44 #include "G4PhysicsLogVector.hh"
45 #include "G4PhysicsVector.hh"
47 #include "G4NistManager.hh"
48 #include "G4SystemOfUnits.hh"
49 
50 #include <fstream>
51 #include <sstream>
52 
53 // factory
54 #include "G4CrossSectionFactory.hh"
55 //
57 
58 using namespace std;
59 
64 
65 #ifdef G4MULTITHREADED
66  G4Mutex G4NeutronElasticXS::neutronElasticXSMutex = G4MUTEX_INITIALIZER;
67 #endif
68 
70  : G4VCrossSectionDataSet(Default_Name()),
71  ggXsection(nullptr),
73  isMaster(false)
74 {
75  // verboseLevel = 0;
76  if(verboseLevel > 0){
77  G4cout << "G4NeutronElasticXS::G4NeutronElasticXS Initialise for Z < "
78  << MAXZEL << G4endl;
79  }
83  temp.resize(13,0.0);
84 }
85 
87 {
88  if(isMaster) {
89  for(G4int i=0; i<MAXZEL; ++i) {
90  delete data[i];
91  data[i] = nullptr;
92  }
93  }
94 }
95 
96 void G4NeutronElasticXS::CrossSectionDescription(std::ostream& outFile) const
97 {
98  outFile << "G4NeutronElasticXS calculates the neutron elastic scattering\n"
99  << "cross section on nuclei using data from the high precision\n"
100  << "neutron database. These data are simplified and smoothed over\n"
101  << "the resonance region in order to reduce CPU time.\n"
102  << "For high energies Glauber-Gribiv cross section is used.\n";
103 }
104 
105 G4bool
107  G4int, const G4Material*)
108 {
109  return true;
110 }
111 
113  G4int, G4int,
114  const G4Element*, const G4Material*)
115 {
116  return true;
117 }
118 
119 G4double
121  G4int ZZ, const G4Material*)
122 {
123  G4double xs = 0.0;
124  G4double ekin = aParticle->GetKineticEnergy();
125 
126  G4int Z = (ZZ >= MAXZEL) ? MAXZEL - 1 : ZZ;
127 
128  auto pv = GetPhysicsVector(Z);
129  if(!pv) { return xs; }
130  // G4cout << "G4NeutronElasticXS::GetCrossSection e= " << ekin
131  // << " Z= " << Z << G4endl;
132 
133  if(ekin <= pv->Energy(0)) {
134  xs = (*pv)[0];
135  } else if(ekin <= pv->GetMaxEnergy()) {
136  xs = pv->LogVectorValue(ekin, aParticle->GetLogKineticEnergy());
137  } else {
139  ekin, Z, aeff[Z]);
140  }
141 
142  if(verboseLevel > 1) {
143  G4cout << "Z= " << Z << " Ekin(MeV)= " << ekin/CLHEP::MeV
144  << ", nElmXSel(b)= " << xs/CLHEP::barn
145  << G4endl;
146  }
147  return xs;
148 }
149 
151  const G4DynamicParticle* aParticle,
152  G4int Z, G4int A,
153  const G4Isotope*, const G4Element*,
154  const G4Material*)
155 {
156  return IsoCrossSection(aParticle->GetKineticEnergy(),
157  aParticle->GetLogKineticEnergy(), Z, A);
158 }
159 
160 G4double
162  G4int ZZ, G4int A)
163 {
164  G4double xs = 0.0;
165  G4int Z = (ZZ >= MAXZEL) ? MAXZEL - 1 : ZZ;
166 
167  // tritium and He3
168  if(3 == A) {
169  return ggXsection->GetElasticElementCrossSection(neutron, ekin, Z, A);
170  }
171  /*
172  G4cout << "IsoCrossSection Z= " << Z << " A= " << A
173  << " Amin= " << amin[Z] << " Amax= " << amax[Z]
174  << " E(MeV)= " << ekin << G4endl;
175  */
176  auto pv = GetPhysicsVector(Z);
177  if(!pv) { return xs; }
178 
179  if(ekin <= pv->Energy(0)) {
180  xs = (*pv)[0];
181  } else if(ekin <= pv->GetMaxEnergy()) {
182  xs = pv->LogVectorValue(ekin, logekin);
183  } else {
185  ekin, Z, aeff[Z]);
186  }
187  xs *= A/aeff[Z];
188  if(verboseLevel > 1) {
189  G4cout << "G4NeutronElasticXS::IsoXS: Z= " << Z << " A= " << A
190  << " Ekin(MeV)= " << ekin/CLHEP::MeV
191  << ", ElmXS(b)= " << xs/CLHEP::barn << G4endl;
192  }
193  return xs;
194 }
195 
197  const G4Element* anElement, G4double kinEnergy, G4double logE)
198 {
199  size_t nIso = anElement->GetNumberOfIsotopes();
200  const G4Isotope* iso = anElement->GetIsotope(0);
201 
202  //G4cout << "SelectIsotope NIso= " << nIso << G4endl;
203  if(1 == nIso) { return iso; }
204 
205  // more than 1 isotope
206  G4int Z = anElement->GetZasInt();
207  //G4cout << "SelectIsotope Z= " << Z << G4endl;
208 
209  const G4double* abundVector = anElement->GetRelativeAbundanceVector();
210  G4double q = G4UniformRand();
211  G4double sum = 0.0;
212  size_t j;
213 
214  // isotope wise cross section not used
215  if(anElement->GetNaturalAbundanceFlag()) {
216  for (j=0; j<nIso; ++j) {
217  sum += abundVector[j];
218  if(q <= sum) {
219  iso = anElement->GetIsotope(j);
220  break;
221  }
222  }
223  return iso;
224  }
225 
226  // use isotope cross sections
227  size_t nn = temp.size();
228  if(nn < nIso) { temp.resize(nIso, 0.); }
229 
230  for (j=0; j<nIso; ++j) {
231  //G4cout << j << "-th isotope " << (*isoVector)[j]->GetN()
232  // << " abund= " << abundVector[j] << G4endl;
233  sum += abundVector[j]*IsoCrossSection(kinEnergy, logE, Z,
234  anElement->GetIsotope(j)->GetN());
235  temp[j] = sum;
236  }
237  sum *= q;
238  for (j = 0; j<nIso; ++j) {
239  if(temp[j] >= sum) {
240  iso = anElement->GetIsotope(j);
241  break;
242  }
243  }
244  return iso;
245 }
246 
247 void
249 {
250  if(verboseLevel > 0){
251  G4cout << "G4NeutronElasticXS::BuildPhysicsTable for "
252  << p.GetParticleName() << G4endl;
253  }
254  if(p.GetParticleName() != "neutron") {
256  ed << p.GetParticleName() << " is a wrong particle type -"
257  << " only neutron is allowed";
258  G4Exception("G4NeutronElasticXS::BuildPhysicsTable(..)","had012",
259  FatalException, ed, "");
260  return;
261  }
262  if(0. == coeff[0]) {
263 #ifdef G4MULTITHREADED
264  G4MUTEXLOCK(&neutronElasticXSMutex);
265  if(0. == coeff[0]) {
266 #endif
267  coeff[0] = 1.0;
268  isMaster = true;
269 #ifdef G4MULTITHREADED
270  }
271  G4MUTEXUNLOCK(&neutronElasticXSMutex);
272 #endif
273  }
274 
275  // it is possible re-initialisation for the second run
276  if(isMaster) {
277 
278  // Access to elements
279  auto theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
280  size_t numOfCouples = theCoupleTable->GetTableSize();
281  for(size_t j=0; j<numOfCouples; ++j) {
282  auto mat = theCoupleTable->GetMaterialCutsCouple(j)->GetMaterial();
283  auto elmVec = mat->GetElementVector();
284  size_t numOfElem = mat->GetNumberOfElements();
285  for (size_t ie = 0; ie < numOfElem; ++ie) {
286  G4int Z = std::max(1,std::min(((*elmVec)[ie])->GetZasInt(), MAXZEL-1));
287  if(!data[Z]) { Initialise(Z); }
288  }
289  }
290  }
291 }
292 
294 {
295  if(!data[Z]) { InitialiseOnFly(Z); }
296  return data[Z];
297 }
298 
300 {
301  // check environment variable
302  // build the complete string identifying the file with the data set
303  if(gDataDirectory.empty()) {
304  char* path = std::getenv("G4PARTICLEXSDATA");
305  if (path) {
306  std::ostringstream ost;
307  ost << path << "/neutron/el";
308  gDataDirectory = ost.str();
309  } else {
310  G4Exception("G4NeutronElasticXS::Initialise(..)","had013",
312  "Environment variable G4PARTICLEXSDATA is not defined");
313  }
314  }
315  return gDataDirectory;
316 }
317 
319 {
320 #ifdef G4MULTITHREADED
321  G4MUTEXLOCK(&neutronElasticXSMutex);
322  if(!data[Z]) {
323 #endif
324  Initialise(Z);
325 #ifdef G4MULTITHREADED
326  }
327  G4MUTEXUNLOCK(&neutronElasticXSMutex);
328 #endif
329 }
330 
332 {
333  if(data[Z]) { return; }
334 
335  // upload data from file
336  data[Z] = new G4PhysicsLogVector();
337 
338  std::ostringstream ost;
339  ost << FindDirectoryPath() << Z ;
340  std::ifstream filein(ost.str().c_str());
341  if (!(filein)) {
343  ed << "Data file <" << ost.str().c_str()
344  << "> is not opened!";
345  G4Exception("G4NeutronElasticXS::Initialise(..)","had014",
346  FatalException, ed, "Check G4PARTICLEXSDATA");
347  return;
348  }
349  if(verboseLevel > 1) {
350  G4cout << "file " << ost.str()
351  << " is opened by G4NeutronElasticXS" << G4endl;
352  }
353 
354  // retrieve data from DB
355  if(!data[Z]->Retrieve(filein, true)) {
357  ed << "Data file <" << ost.str().c_str()
358  << "> is not retrieved!";
359  G4Exception("G4NeutronElasticXS::Initialise(..)","had015",
360  FatalException, ed, "Check G4PARTICLEXSDATA");
361  return;
362  }
363  // smooth transition
364  G4double sig1 = (*(data[Z]))[data[Z]->GetVectorLength()-1];
365  G4double ehigh = data[Z]->GetMaxEnergy();
366  aeff[Z] = nist->GetAtomicMassAmu(Z);
368  ehigh, Z, aeff[Z]);
369  if(sig2 > 0.) { coeff[Z] = sig1/sig2; }
370 }