ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleInelasticXS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 // -------------------------------------------------------------------
27 //
28 // GEANT4 Class file
29 //
30 //
31 // File name: G4ParticleInelasticXS
32 //
33 // Author Ivantchenko, Geant4, 3-Aug-09
34 //
35 // Modifications:
36 //
37 
38 #include "G4ParticleInelasticXS.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"
48 #include "G4NistManager.hh"
49 #include "G4Proton.hh"
50 #include "Randomize.hh"
51 #include "G4SystemOfUnits.hh"
52 
53 #include <fstream>
54 #include <sstream>
55 
56 using namespace std;
57 
59  0,
60  1, 4, 6, 9, 10, 12, 14, 16, 19, 20, //1-10
61  23, 24, 27, 28, 31, 32, 35, 36, 39, 40, //11-20
62  45, 46, 50, 50, 55, 54, 59, 58, 63, 64, //21-30
63  69, 70, 75, 0, 0, 0, 0, 0, 0, 90, //31-40
64  0, 92, 0, 0, 0, 102, 107, 106, 113, 112, //41-50
65  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //51-60
66  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //61-70
67  0, 0, 181, 180, 0, 0, 0, 192, 197, 0, //71-80
68  0, 204, 209, 0, 0, 0, 0, 0, 0, 0, //81-90
69  0, 235};
71  0,
72  2, 4, 7, 9, 11, 13, 15, 18, 19, 22, //1-10
73  23, 26, 27, 30, 31, 34, 37, 40, 41, 48, //11-20
74  45, 50, 51, 54, 55, 58, 59, 64, 65, 70, //21-30
75  71, 76, 75, 0, 0, 0, 0, 0, 0, 96, //31-40
76  0, 100, 0, 0, 0, 110, 109, 116, 115, 124, //41-50
77  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //51-60
78  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //61-70
79  0, 0, 181, 186, 0, 0, 0, 198, 197, 0, //71-80
80  0, 208, 209, 0, 0, 0, 0, 0, 0, 0, //81-90
81  0, 238};
82 
87 
88 #ifdef G4MULTITHREADED
89  G4Mutex G4ParticleInelasticXS::particleInelasticXSMutex = G4MUTEX_INITIALIZER;
90 #endif
91 
93  : G4VCrossSectionDataSet("G4ParticleInelasticXS"),
94  ggXsection(nullptr),
95  nnXsection(nullptr),
96  particle(part),
98  isMaster(false)
99 {
100  if(!part) {
101  G4Exception("G4ParticleInelasticXS::G4ParticleInelasticXS(..)","had015",
102  FatalException, "NO particle definition in constructor");
103  } else {
104  verboseLevel = 0;
105  G4String particleName = particle->GetParticleName();
106  if(verboseLevel > 0){
107  G4cout << "G4ParticleInelasticXS::G4ParticleInelasticXS for "
108  << particleName << " on atoms with Z < " << MAXZINELP << G4endl;
109  }
110  if(particleName == "neutron" || particleName == "proton") {
112  } else {
114  }
115  }
118  temp.resize(13,0.0);
119 }
120 
122 {
123  if(isMaster) { delete data; data = nullptr; }
124 }
125 
126 void G4ParticleInelasticXS::CrossSectionDescription(std::ostream& outFile) const
127 {
128  outFile << "G4ParticleInelasticXS calculates n, p, d, t, he3, he4 inelastic \n"
129  << "cross section on nuclei using data from the high precision\n"
130  << "neutron database. These data are simplified and smoothed over\n"
131  << "the resonance region in order to reduce CPU time.\n"
132  << "For high energy Glauber-Gribov cross section model is used\n";
133 }
134 
135 G4bool
137  G4int, const G4Material*)
138 {
139  return true;
140 }
141 
142 G4bool
144  G4int, G4int,
145  const G4Element*, const G4Material*)
146 {
147  return true;
148 }
149 
151  const G4DynamicParticle* aParticle,
152  G4int ZZ, const G4Material*)
153 {
154  G4double xs = 0.0;
155  G4double ekin = aParticle->GetKineticEnergy();
156 
157  G4int Z = (ZZ >= MAXZINELP) ? MAXZINELP - 1 : ZZ;
158 
159  auto pv = GetPhysicsVector(Z);
160  if(!pv) { return xs; }
161  // G4cout << "G4ParticleInelasticXS::GetCrossSection e= " << ekin
162  // << " Z= " << Z << G4endl;
163 
164  // below threshold
165  if(ekin <= pv->Energy(0)) { return xs; }
166 
167  if(ekin <= pv->GetMaxEnergy()) {
168  xs = pv->LogVectorValue(ekin, aParticle->GetLogKineticEnergy());
169  } else {
170  if(ggXsection) {
172  ekin, Z, aeff[Z]);
173  } else {
175  ekin, Z, aeff[Z]);
176  }
177  }
178 
179  if(verboseLevel > 1) {
180  G4cout << "ElmXS: Z= " << Z << " Ekin(MeV)= " << ekin/CLHEP::MeV
181  << " xs(bn)= " << xs/CLHEP::barn << " element data for "
183  }
184  return xs;
185 }
186 
188  const G4DynamicParticle* aParticle,
189  G4int Z, G4int A,
190  const G4Isotope*, const G4Element*,
191  const G4Material*)
192 {
193  return IsoCrossSection(aParticle->GetKineticEnergy(),
194  aParticle->GetLogKineticEnergy(),Z, A);
195 }
196 
197 G4double
199  G4int ZZ, G4int A)
200 {
201  G4double xs = 0.0;
202  G4int Z = (ZZ >= MAXZINELP) ? MAXZINELP - 1 : ZZ;
203 
204  // tritium and He3
205  if(3 == A) {
206  if(ggXsection) {
208  } else {
210  }
211  return xs;
212  }
213  /*
214  G4cout << "IsoCrossSection Z= " << Z << " A= " << A
215  << " Amin= " << amin[Z] << " Amax= " << amax[Z]
216  << " E(MeV)= " << ekin << G4endl;
217  */
218  auto pv = GetPhysicsVector(Z);
219  if(!pv) { return xs; }
220 
221  // below threshold
222  if(ekin <= pv->Energy(0)) { return xs; }
223 
224  // compute isotope cross section if applicable
225  G4double emax = pv->GetMaxEnergy();
226  if(ekin <= emax && amin[Z]>0 && A >= amin[Z] && A <= amax[Z]) {
227  auto pviso = data->GetComponentDataByIndex(Z, A - amin[Z]);
228  if(pviso) {
229  xs = pviso->LogVectorValue(ekin, logE);
230  if(verboseLevel > 1) {
231  G4cout << "G4ParticleInelasticXS::IsoXS: for "
232  << particle->GetParticleName() << " Ekin(MeV)= "
233  << ekin/CLHEP::MeV << " xs(b)= " << xs/CLHEP::barn
234  << " Z= " << Z << " A= " << A << G4endl;
235  }
236  return xs;
237  }
238  }
239  // use element x-section
240  if(ekin <= emax) {
241  xs = pv->LogVectorValue(ekin, logE);
242  } else {
243  if(ggXsection) {
245  ekin, Z, aeff[Z]);
246  } else {
248  ekin, Z, aeff[Z]);
249  }
250  }
251  xs *= A/aeff[Z];
252  if(verboseLevel > 1) {
253  G4cout << "IsoXS for " << particle->GetParticleName()
254  << " Target Z= " << Z << " A= " << A
255  << " Ekin(MeV)= " << ekin/CLHEP::MeV
256  << " xs(bn)= " << xs/CLHEP::barn << G4endl;
257  }
258  return xs;
259 }
260 
262  const G4Element* anElement, G4double kinEnergy, G4double logE)
263 {
264  size_t nIso = anElement->GetNumberOfIsotopes();
265  const G4Isotope* iso = anElement->GetIsotope(0);
266 
267  //G4cout << "SelectIsotope NIso= " << nIso << G4endl;
268  if(1 == nIso) { return iso; }
269 
270  // more than 1 isotope
271  G4int Z = anElement->GetZasInt();
272  //G4cout << "SelectIsotope Z= " << Z << G4endl;
273 
274  const G4double* abundVector = anElement->GetRelativeAbundanceVector();
275  G4double q = G4UniformRand();
276  G4double sum = 0.0;
277  size_t j;
278 
279  // isotope wise cross section not available
280  if(0 == amin[Z] || Z >= MAXZINELP) {
281  for (j=0; j<nIso; ++j) {
282  sum += abundVector[j];
283  if(q <= sum) {
284  iso = anElement->GetIsotope(j);
285  break;
286  }
287  }
288  return iso;
289  }
290 
291  size_t nn = temp.size();
292  if(nn < nIso) { temp.resize(nIso, 0.); }
293 
294  for (j=0; j<nIso; ++j) {
295  //G4cout << j << "-th isotope " << (*isoVector)[j]->GetN()
296  // << " abund= " << abundVector[j] << G4endl;
297  sum += abundVector[j]*IsoCrossSection(kinEnergy, logE, Z,
298  anElement->GetIsotope(j)->GetN());
299  temp[j] = sum;
300  }
301  sum *= q;
302  for (j=0; j<nIso; ++j) {
303  if(temp[j] >= sum) {
304  iso = anElement->GetIsotope(j);
305  break;
306  }
307  }
308  return iso;
309 }
310 
311 void
313 {
314  if(verboseLevel > 0){
315  G4cout << "G4ParticleInelasticXS::BuildPhysicsTable for "
316  << p.GetParticleName() << G4endl;
317  }
318  if(&p != particle) {
320  ed << p.GetParticleName() << " is a wrong particle type -"
321  << particle->GetParticleName() << " is expected";
322  G4Exception("G4ParticleInelasticXS::BuildPhysicsTable(..)","had012",
323  FatalException, ed, "");
324  return;
325  }
326 
327  if(!data) {
328 #ifdef G4MULTITHREADED
329  G4MUTEXLOCK(&particleInelasticXSMutex);
330  if(!data) {
331 #endif
332  isMaster = true;
333  data = new G4ElementData();
334  data->SetName(particle->GetParticleName() + "Inelastic");
335 #ifdef G4MULTITHREADED
336  }
337  G4MUTEXUNLOCK(&particleInelasticXSMutex);
338 #endif
339  }
340 
341  // it is possible re-initialisation for the new run
342  if(isMaster) {
343 
344  // Access to elements
345  auto theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
346  size_t numOfCouples = theCoupleTable->GetTableSize();
347  for(size_t j=0; j<numOfCouples; ++j) {
348  auto mat = theCoupleTable->GetMaterialCutsCouple(j)->GetMaterial();
349  auto elmVec = mat->GetElementVector();
350  size_t numOfElem = mat->GetNumberOfElements();
351  for (size_t ie = 0; ie < numOfElem; ++ie) {
352  G4int Z = std::max(1,std::min(((*elmVec)[ie])->GetZasInt(), MAXZINELP-1));
353  if(!data->GetElementData(Z)) { Initialise(Z); }
354  }
355  }
356  }
357 }
358 
360 {
361  const G4PhysicsVector* pv = data->GetElementData(Z);
362  if(!pv) {
363  InitialiseOnFly(Z);
364  pv = data->GetElementData(Z);
365  }
366  return pv;
367 }
368 
370 {
371  // check environment variable
372  // build the complete string identifying the file with the data set
373  if(gDataDirectory.empty()) {
374  char* path = std::getenv("G4PARTICLEXSDATA");
375  if (path) {
376  std::ostringstream ost;
377  ost << path << "/" << particle->GetParticleName() << "/inel";
378  gDataDirectory = ost.str();
379  } else {
380  G4Exception("G4NeutronInelasticXS::Initialise(..)","had013",
382  "Environment variable G4PARTICLEXSDATA is not defined");
383  }
384  }
385  return gDataDirectory;
386 }
387 
389 {
390 #ifdef G4MULTITHREADED
391  G4MUTEXLOCK(&particleInelasticXSMutex);
392  if(!data->GetElementData(Z)) {
393 #endif
394  Initialise(Z);
395 #ifdef G4MULTITHREADED
396  }
397  G4MUTEXUNLOCK(&particleInelasticXSMutex);
398 #endif
399 }
400 
402 {
403  if(data->GetElementData(Z)) { return; }
404 
405  // upload element data
406  std::ostringstream ost;
407  ost << FindDirectoryPath() << Z ;
408  G4PhysicsVector* v = RetrieveVector(ost, true);
409  data->InitialiseForElement(Z, v);
410  /*
411  G4cout << "G4ParticleInelasticXS::Initialise for Z= " << Z
412  << " A= " << Amean << " Amin= " << amin[Z]
413  << " Amax= " << amax[Z] << G4endl;
414  */
415  // upload isotope data
416  if(amin[Z] > 0) {
417  size_t nmax = (size_t)(amax[Z]-amin[Z]+1);
418  data->InitialiseForComponent(Z, nmax);
419 
420  for(G4int A=amin[Z]; A<=amax[Z]; ++A) {
421  std::ostringstream ost1;
422  ost1 << gDataDirectory << Z << "_" << A;
423  G4PhysicsVector* v1 = RetrieveVector(ost1, false);
424  data->AddComponent(Z, A, v1);
425  }
426  }
427  // smooth transition
428  G4double sig1 = (*v)[v->GetVectorLength()-1];
429  G4double sig2 = 0.0;
430  G4double ehigh = v->GetMaxEnergy();
431  aeff[Z] = fNist->GetAtomicMassAmu(Z);
432  if(ggXsection) {
434  ehigh, Z, aeff[Z]);
435  } else {
437  ehigh, Z, aeff[Z]);
438  }
439  if(sig2 > 0.) { coeff[Z] = sig1/sig2; }
440 }
441 
443 G4ParticleInelasticXS::RetrieveVector(std::ostringstream& ost, G4bool warn)
444 {
445  G4PhysicsLogVector* v = nullptr;
446  std::ifstream filein(ost.str().c_str());
447  if (!(filein)) {
448  if(warn) {
450  ed << "Data file <" << ost.str().c_str()
451  << "> is not opened!";
452  G4Exception("G4ParticleInelasticXS::RetrieveVector(..)","had014",
453  FatalException, ed, "Check G4PARTICLEXSDATA");
454  }
455  } else {
456  if(verboseLevel > 1) {
457  G4cout << "File " << ost.str()
458  << " is opened by G4ParticleInelasticXS" << G4endl;
459  }
460  // retrieve data from DB
461  v = new G4PhysicsLogVector();
462  if(!v->Retrieve(filein, true)) {
464  ed << "Data file <" << ost.str().c_str()
465  << "> is not retrieved!";
466  G4Exception("G4ParticleInelasticXS::RetrieveVector(..)","had015",
467  FatalException, ed, "Check G4PARTICLEXSDATA");
468  }
469  }
470  return v;
471 }