ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NeutronInelasticXS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 // -------------------------------------------------------------------
27 //
28 // GEANT4 Class file
29 //
30 //
31 // File name: G4NeutronInelasticXS
32 //
33 // Author Ivantchenko, Geant4, 3-Aug-09
34 //
35 
36 
37 #include "G4NeutronInelasticXS.hh"
38 #include "G4Neutron.hh"
39 #include "G4DynamicParticle.hh"
40 #include "G4ProductionCutsTable.hh"
41 #include "G4Material.hh"
42 #include "G4Element.hh"
43 #include "G4PhysicsLogVector.hh"
44 #include "G4PhysicsVector.hh"
46 #include "G4NistManager.hh"
47 #include "Randomize.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 
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  2, 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 
89 
90 #ifdef G4MULTITHREADED
91  G4Mutex G4NeutronInelasticXS::neutronInelasticXSMutex = G4MUTEX_INITIALIZER;
92 #endif
93 
95  : G4VCrossSectionDataSet(Default_Name()),
97 {
98  // verboseLevel = 0;
99  if(verboseLevel > 0){
100  G4cout << "G4NeutronInelasticXS::G4NeutronInelasticXS Initialise for Z < "
101  << MAXZINEL << G4endl;
102  }
106  isMaster = false;
107  temp.resize(13,0.0);
108 }
109 
111 {
112  if(isMaster) { delete data; data = nullptr; }
113 }
114 
115 void G4NeutronInelasticXS::CrossSectionDescription(std::ostream& outFile) const
116 {
117  outFile << "G4NeutronInelasticXS calculates the neutron inelastic scattering\n"
118  << "cross section on nuclei using data from the high precision\n"
119  << "neutron database. These data are simplified and smoothed over\n"
120  << "the resonance region in order to reduce CPU time.\n"
121  << "For high energy Glauber-Gribov cross section model is used\n";
122 }
123 
124 G4bool
126  G4int, const G4Material*)
127 {
128  return true;
129 }
130 
131 G4bool
133  G4int, G4int,
134  const G4Element*, const G4Material*)
135 {
136  return true;
137 }
138 
140  const G4DynamicParticle* aParticle,
141  G4int ZZ, const G4Material*)
142 {
143  G4double xs = 0.0;
144  G4double ekin = aParticle->GetKineticEnergy();
145 
146  G4int Z = (ZZ >= MAXZINEL) ? MAXZINEL - 1 : ZZ;
147 
148  auto pv = GetPhysicsVector(Z);
149  if(!pv) { return xs; }
150  // G4cout << "G4NeutronInelasticXS::GetCrossSection e= " << ekin
151  // << " Z= " << Z << G4endl;
152 
153  // below threshold
154  if(ekin <= pv->Energy(0)) { return xs; }
155 
156  if(ekin <= pv->GetMaxEnergy()) {
157  xs = pv->LogVectorValue(ekin, aParticle->GetLogKineticEnergy());
158  } else {
160  ekin, Z, aeff[Z]);
161  }
162 
163  if(verboseLevel > 1) {
164  G4cout << "Z= " << Z << " Ekin(MeV)= " << ekin/CLHEP::MeV
165  << ", ElmXSinel(b)= " << xs/CLHEP::barn
166  << G4endl;
167  }
168  return xs;
169 }
170 
172  const G4DynamicParticle* aParticle,
173  G4int Z, G4int A,
174  const G4Isotope*, const G4Element*,
175  const G4Material*)
176 {
177  return IsoCrossSection(aParticle->GetKineticEnergy(),
178  aParticle->GetLogKineticEnergy(), Z, A);
179 }
180 
181 G4double
183  G4int ZZ, G4int A)
184 {
185  G4double xs = 0.0;
186  G4int Z = (ZZ >= MAXZINEL) ? MAXZINEL - 1 : ZZ;
187 
188  // tritium and He3
189  if(3 == A) {
191  }
192 
193  /*
194  G4cout << "IsoCrossSection Z= " << Z << " A= " << A
195  << " Amin= " << amin[Z] << " Amax= " << amax[Z]
196  << " E(MeV)= " << ekin << G4endl;
197  */
198  auto pv = GetPhysicsVector(Z);
199  if(!pv) { return xs; }
200 
201  // below threshold
202  if(ekin <= pv->Energy(0)) { return xs; }
203 
204  // compute isotope cross section if applicable
205  G4double emax = pv->GetMaxEnergy();
206  if(ekin <= emax && amin[Z] > 0 && A >= amin[Z] && A <= amax[Z]) {
207  auto pviso = data->GetComponentDataByIndex(Z, A - amin[Z]);
208  if(pviso) {
209  xs = pviso->LogVectorValue(ekin, logekin);
210  if(verboseLevel > 1) {
211  G4cout << "G4NeutronInelasticXS::IsoXS: Ekin(MeV)= "
212  << ekin/CLHEP::MeV
213  << " xs(b)= " << xs/CLHEP::barn
214  << " Z= " << Z << " A= " << A << G4endl;
215  }
216  return xs;
217  }
218  }
219 
220  // use element x-section
221  if(ekin <= emax) {
222  xs = pv->LogVectorValue(ekin, logekin);
223  } else {
225  ekin, Z, aeff[Z]);
226  }
227  xs *= A/aeff[Z];
228  if(verboseLevel > 1) {
229  G4cout << "G4NeutronInelasticXS::IsoXS: Z= " << Z << " A= " << A
230  << " Ekin(MeV)= " << ekin/CLHEP::MeV
231  << ", ElmXS(b)= " << xs/CLHEP::barn << G4endl;
232  }
233  return xs;
234 }
235 
237  const G4Element* anElement, G4double kinEnergy, G4double logE)
238 {
239  size_t nIso = anElement->GetNumberOfIsotopes();
240  const G4Isotope* iso = anElement->GetIsotope(0);
241 
242  //G4cout << "SelectIsotope NIso= " << nIso << G4endl;
243  if(1 == nIso) { return iso; }
244 
245  // more than 1 isotope
246  G4int Z = anElement->GetZasInt();
247  //G4cout << "SelectIsotope Z= " << Z << G4endl;
248 
249  const G4double* abundVector = anElement->GetRelativeAbundanceVector();
250  G4double q = G4UniformRand();
251  G4double sum = 0.0;
252  size_t j;
253 
254  // isotope wise cross section not available
255  if(0 == amin[Z] || Z >= MAXZINEL) {
256  for (j=0; j<nIso; ++j) {
257  sum += abundVector[j];
258  if(q <= sum) {
259  iso = anElement->GetIsotope(j);
260  break;
261  }
262  }
263  return iso;
264  }
265 
266  // use isotope cross sections
267  size_t nn = temp.size();
268  if(nn < nIso) { temp.resize(nIso, 0.); }
269 
270  for (j=0; j<nIso; ++j) {
271  //G4cout << j << "-th isotope " << (*isoVector)[j]->GetN()
272  // << " abund= " << abundVector[j] << G4endl;
273  sum += abundVector[j]*IsoCrossSection(kinEnergy, logE, Z,
274  anElement->GetIsotope(j)->GetN());
275  temp[j] = sum;
276  }
277  sum *= q;
278  for (j = 0; j<nIso; ++j) {
279  if(temp[j] >= sum) {
280  iso = anElement->GetIsotope(j);
281  break;
282  }
283  }
284  return iso;
285 }
286 
287 void
289 {
290  if(verboseLevel > 0) {
291  G4cout << "G4NeutronInelasticXS::BuildPhysicsTable for "
292  << p.GetParticleName() << G4endl;
293  }
294  if(p.GetParticleName() != "neutron") {
296  ed << p.GetParticleName() << " is a wrong particle type -"
297  << " only neutron is allowed";
298  G4Exception("G4NeutronInelasticXS::BuildPhysicsTable(..)","had012",
299  FatalException, ed, "");
300  return;
301  }
302 
303  if(!data) {
304 #ifdef G4MULTITHREADED
305  G4MUTEXLOCK(&neutronInelasticXSMutex);
306  if(!data) {
307 #endif
308  isMaster = true;
309  data = new G4ElementData();
310  data->SetName("NeutronInelastic");
312 #ifdef G4MULTITHREADED
313  }
314  G4MUTEXUNLOCK(&neutronInelasticXSMutex);
315 #endif
316  }
317 
318  // it is possible re-initialisation for the new run
319  if(isMaster) {
320 
321  // Access to elements
322  auto theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
323  size_t numOfCouples = theCoupleTable->GetTableSize();
324  for(size_t j=0; j<numOfCouples; ++j) {
325  auto mat = theCoupleTable->GetMaterialCutsCouple(j)->GetMaterial();
326  auto elmVec = mat->GetElementVector();
327  size_t numOfElem = mat->GetNumberOfElements();
328  for (size_t ie = 0; ie < numOfElem; ++ie) {
329  G4int Z = std::max(1,std::min(((*elmVec)[ie])->GetZasInt(), MAXZINEL-1));
330  if(!data->GetElementData(Z)) { Initialise(Z); }
331  }
332  }
333  }
334 }
335 
337 {
338  const G4PhysicsVector* pv = data->GetElementData(Z);
339  if(!pv) {
340  InitialiseOnFly(Z);
341  pv = data->GetElementData(Z);
342  }
343  return pv;
344 }
345 
347 {
348  // check environment variable
349  // build the complete string identifying the file with the data set
350  if(gDataDirectory.empty()) {
351  char* path = std::getenv("G4PARTICLEXSDATA");
352  if (path) {
353  std::ostringstream ost;
354  ost << path << "/neutron/inel";
355  gDataDirectory = ost.str();
356  } else {
357  G4Exception("G4NeutronInelasticXS::Initialise(..)","had013",
359  "Environment variable G4PARTICLEXSDATA is not defined");
360  }
361  }
362  return gDataDirectory;
363 }
364 
366 {
367 #ifdef G4MULTITHREADED
368  G4MUTEXLOCK(&neutronInelasticXSMutex);
369  if(!data->GetElementData(Z)) {
370 #endif
371  Initialise(Z);
372 #ifdef G4MULTITHREADED
373  }
374  G4MUTEXUNLOCK(&neutronInelasticXSMutex);
375 #endif
376 }
377 
379 {
380  if(data->GetElementData(Z)) { return; }
381 
382  // upload element data
383  std::ostringstream ost;
384  ost << FindDirectoryPath() << Z;
385  G4PhysicsVector* v = RetrieveVector(ost, true);
386  data->InitialiseForElement(Z, v);
387  /*
388  G4cout << "G4NeutronInelasticXS::Initialise for Z= " << Z
389  << " A= " << Amean << " Amin= " << amin[Z]
390  << " Amax= " << amax[Z] << G4endl;
391  */
392  // upload isotope data
393  if(amin[Z] > 0) {
394  size_t nmax = (size_t)(amax[Z]-amin[Z]+1);
395  data->InitialiseForComponent(Z, nmax);
396 
397  for(G4int A=amin[Z]; A<=amax[Z]; ++A) {
398  std::ostringstream ost1;
399  ost1 << gDataDirectory << Z << "_" << A;
400  G4PhysicsVector* v1 = RetrieveVector(ost1, false);
401  data->AddComponent(Z, A, v1);
402  }
403  }
404 
405  // smooth transition
406  G4double sig1 = (*v)[v->GetVectorLength()-1];
407  G4double ehigh= v->GetMaxEnergy();
408  aeff[Z] = nist->GetAtomicMassAmu(Z);
410  ehigh, Z, aeff[Z]);
411  if(sig2 > 0.) { coeff[Z] = sig1/sig2; }
412 }
413 
415 G4NeutronInelasticXS::RetrieveVector(std::ostringstream& ost, G4bool warn)
416 {
417  G4PhysicsLogVector* v = nullptr;
418  std::ifstream filein(ost.str().c_str());
419  if (!(filein)) {
420  if(warn) {
422  ed << "Data file <" << ost.str().c_str()
423  << "> is not opened!";
424  G4Exception("G4NeutronInelasticXS::RetrieveVector(..)","had014",
425  FatalException, ed, "Check G4PARTICLEXSDATA");
426  }
427  } else {
428  if(verboseLevel > 1) {
429  G4cout << "File " << ost.str()
430  << " is opened by G4NeutronInelasticXS" << G4endl;
431  }
432  // retrieve data from DB
433  v = new G4PhysicsLogVector();
434  if(!v->Retrieve(filein, true)) {
436  ed << "Data file <" << ost.str().c_str()
437  << "> is not retrieved!";
438  G4Exception("G4NeutronInelasticXS::RetrieveVector(..)","had015",
439  FatalException, ed, "Check G4PARTICLEXSDATA");
440  }
441  }
442  return v;
443 }