| Geant4 Cross Reference |
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 // V. Ivanchenko September 2023
27 //
28 // G4CrossSectionHP is a generic class implementing
29 // cross sections for neutrons, protons and light ions
30 // It is an alternative to code developed by J.P. Wellisch & T.Koi
31 //
32
33 #include <fstream>
34 #include <sstream>
35 #include <thread>
36
37 #include "G4CrossSectionHP.hh"
38 #include "G4Material.hh"
39 #include "G4Element.hh"
40 #include "G4ElementDataRegistry.hh"
41 #include "G4ElementTable.hh"
42 #include "G4IsotopeList.hh"
43 #include "G4HadronicParameters.hh"
44 #include "G4ParticleHPManager.hh"
45 #include "G4ParticleDefinition.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4Pow.hh"
48 #include "G4Log.hh"
49 #include "Randomize.hh"
50 #include "G4RandomDirection.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4Neutron.hh"
53 #include "G4NucleiProperties.hh"
54 #include "G4AutoLock.hh"
55
56 namespace
57 {
58 G4Mutex theHPXS = G4MUTEX_INITIALIZER;
59 }
60
61 G4CrossSectionHP::G4CrossSectionHP(const G4ParticleDefinition* p,
62 const G4String& nameData,
63 const G4String& nameDir, G4double emaxHP,
64 G4int zmin, G4int zmax)
65 : G4VCrossSectionDataSet(nameData),
66 fParticle(p),
67 fNeutron(G4Neutron::Neutron()),
68 fManagerHP(G4ParticleHPManager::GetInstance()),
69 emax(emaxHP),
70 emaxT(fManagerHP->GetMaxEnergyDoppler()),
71 elimit(1.0e-11*CLHEP::eV),
72 logElimit(G4Log(elimit)),
73 minZ(zmin),
74 maxZ(zmax),
75 fDataName(nameData),
76 fDataDirectory(nameDir)
77 {
78 // verboseLevel = 1;
79 if (verboseLevel > 1) {
80 G4cout << "G4CrossSectionHP::G4CrossSectionHP: Initialise for "
81 << fDataName << " " << minZ << " < Z < " << maxZ
82 << " EmaxT(MeV)=" << emaxT << G4endl;
83 G4cout << "Data directory: " << fDataDirectory << G4endl;
84 }
85 auto ptr = G4ElementDataRegistry::Instance();
86 auto data = ptr->GetElementDataByName(fDataName);
87 if (nullptr == data) {
88 data = new G4ElementData(maxZ - minZ + 1);
89 data->SetName(fDataName);
90 }
91 fData = data;
92 }
93
94 G4bool G4CrossSectionHP::IsIsoApplicable(const G4DynamicParticle* dp,
95 G4int, G4int,
96 const G4Element*,
97 const G4Material*)
98 {
99 return (dp->GetKineticEnergy() <= emax);
100 }
101
102 G4double G4CrossSectionHP::GetIsoCrossSection(const G4DynamicParticle* dp,
103 G4int Z, G4int A,
104 const G4Isotope*,
105 const G4Element*,
106 const G4Material* mat)
107 {
108 G4double ekin = dp->GetKineticEnergy();
109 G4double loge = dp->GetLogKineticEnergy();
110 if (ekin < elimit) {
111 ekin = elimit;
112 loge = logElimit;
113 }
114 if (mat != fCurrentMat) { PrepareCache(mat); }
115
116 return IsoCrossSection(ekin, loge, Z, A, mat->GetTemperature());
117 }
118
119 G4double
120 G4CrossSectionHP::ComputeIsoCrossSection(G4double e, G4double le,
121 const G4ParticleDefinition*,
122 G4int Z, G4int A,
123 const G4Isotope*,
124 const G4Element*,
125 const G4Material* mat)
126 {
127 G4double ekin = e;
128 G4double loge = le;
129 if (ekin < elimit) {
130 ekin = elimit;
131 loge = logElimit;
132 }
133 if (mat != fCurrentMat) { PrepareCache(mat); }
134
135 return IsoCrossSection(ekin, loge, Z, A, mat->GetTemperature());
136 }
137
138 G4double G4CrossSectionHP::IsoCrossSection(const G4double ekin,
139 const G4double logek,
140 const G4int Z, const G4int A,
141 const G4double T)
142 {
143 // G4cout << "G4CrossSectionHP::IsoCrossSection Z=" << Z << " A=" << A
144 // << " E(MeV)=" << ekin/MeV << " T=" << T << " " << GetName() << G4endl;
145 G4double xs = 0.0;
146 if (ekin > emax || Z > maxZ || Z < minZ || ekin < elimit) { return xs; }
147
148 auto pv0 = fData->GetElementData(Z - minZ);
149 if (nullptr == pv0) {
150 Initialise(Z);
151 pv0 = fData->GetElementData(Z - minZ);
152 }
153
154 // if there is no element data then no iso data
155 if (nullptr == pv0) { return xs; }
156
157 const auto pv = fData->GetComponentDataByID(Z - minZ, A);
158 if (nullptr == pv) { return xs; }
159
160 // no Doppler broading
161 // G4double factT = T/CLHEP::STP_Temperature;
162 if (ekin >= emaxT*T/CLHEP::STP_Temperature || fManagerHP->GetNeglectDoppler()) {
163 xs = pv->LogFreeVectorValue(ekin, logek);
164
165 } else {
166
167 // Doppler broading
168 G4double e0 = CLHEP::k_Boltzmann*T;
169 G4double mass = fParticle->GetPDGMass();
170 G4double massTarget = G4NucleiProperties::GetNuclearMass(A, Z);
171 G4double sig = std::sqrt(2.0*e0/(3.0*massTarget));
172
173 // projectile
174 G4LorentzVector lv(0., 0., std::sqrt(ekin*(ekin + 2*mass)), mass + ekin);
175
176 // limits of integration
177 const G4double lim = 1.01;
178 const G4int nmin = 3;
179 G4int ii = 0;
180 const G4int nn = 20;
181 G4double xs2 = 0.0;
182
183 for (G4int i=0; i<nn; ++i) {
184 G4double vx = G4RandGauss::shoot(0., sig);
185 G4double vy = G4RandGauss::shoot(0., sig);
186 G4double vz = G4RandGauss::shoot(0., sig);
187 fLV.set(massTarget*vx, massTarget*vy, massTarget*vz, massTarget*(1.0 + 0.5*(vx*vx + vy*vy + vz*vz)));
188 fBoost = fLV.boostVector();
189 fLV = lv.boost(-fBoost);
190 if (fLV.pz() <= 0.0) { continue; }
191 ++ii;
192 G4double e = fLV.e() - mass;
193 G4double y = pv->Value(e, index);
194 xs += y;
195 xs2 += y*y;
196 if (ii >= nmin && ii*xs2 <= lim*xs*xs) { break; }
197 }
198 if (ii > 0) { xs /= (G4double)(ii); }
199 }
200 #ifdef G4VERBOSE
201 if (verboseLevel > 1) {
202 G4cout << "G4CrossSectionHP::IsoXS " << fDataName
203 << " Z= " << Z << " A= " << A
204 << " Ekin(MeV)= " << ekin/MeV << " xs(b)= " << xs/barn
205 << " size=" << fZA.size() << G4endl;
206 }
207 #endif
208
209 // save cross section into struct
210 for (std::size_t i=0; i<fZA.size(); ++i) {
211 if (Z == fZA[i].first && A == fZA[i].second) {
212 fIsoXS[i] = xs;
213 break;
214 }
215 }
216 return xs;
217 }
218
219 const G4Isotope* G4CrossSectionHP::SelectIsotope(const G4Element* elm,
220 G4double, G4double)
221 {
222 std::size_t nIso = elm->GetNumberOfIsotopes();
223 const G4Isotope* iso = elm->GetIsotope(0);
224
225 //G4cout << "SelectIsotope NIso= " << nIso << G4endl;
226 if(1 == nIso) { return iso; }
227
228 // more than 1 isotope
229 G4int Z = elm->GetZasInt();
230 if (Z >= minZ && Z <= maxZ && nullptr == fData->GetElementData(Z - minZ)) {
231 Initialise(Z);
232 }
233
234 const G4double* abundVector = elm->GetRelativeAbundanceVector();
235 G4double q = G4UniformRand();
236 G4double sum = 0.0;
237
238 // is there cached isotope wise cross section?
239 std::size_t j;
240 if (Z < minZ || Z > maxZ || !CheckCache(Z) ||
241 0 == fData->GetNumberOfComponents(Z - minZ)) {
242 for (j = 0; j<nIso; ++j) {
243 sum += abundVector[j];
244 if(q <= sum) {
245 iso = elm->GetIsotope((G4int)j);
246 break;
247 }
248 }
249 return iso;
250 }
251 std::size_t nn = fTemp.size();
252 if (nn < nIso) { fTemp.resize(nIso, 0.); }
253
254 // reuse cache
255 for (j=0; j<nIso; ++j) {
256 sum += abundVector[j]*
257 GetCrossSection(Z - minZ, elm->GetIsotope((G4int)j)->GetN());
258 fTemp[j] = sum;
259 }
260 sum *= q;
261 for (j = 0; j<nIso; ++j) {
262 if (fTemp[j] >= sum) {
263 iso = elm->GetIsotope((G4int)j);
264 break;
265 }
266 }
267 return iso;
268 }
269
270 void G4CrossSectionHP::BuildPhysicsTable(const G4ParticleDefinition& p)
271 {
272 if (verboseLevel > 1) {
273 G4cout << "G4CrossSectionHP::BuildPhysicsTable for "
274 << p.GetParticleName() << " and " << fDataName << G4endl;
275 }
276
277 // it is possible re-initialisation for the second run
278 const G4ElementTable* table = G4Element::GetElementTable();
279
280 // Access to elements
281 for ( auto const & elm : *table ) {
282 G4int Z = elm->GetZasInt();
283 if (Z >= minZ && Z <= maxZ && nullptr == fData->GetElementData(Z - minZ)) {
284 Initialise(Z);
285 }
286 }
287
288 // prepare isotope selection
289 std::size_t nmax = 0;
290 std::size_t imax = 0;
291 for ( auto const & mat : *G4Material::GetMaterialTable() ) {
292 std::size_t n = 0;
293 for ( auto const & elm : *mat->GetElementVector() ) {
294 std::size_t niso = elm->GetNumberOfIsotopes();
295 n += niso;
296 if(niso > imax) { imax = niso; }
297 }
298 if (n > nmax) { nmax = n; }
299 }
300 fTemp.resize(imax, 0.0);
301 fZA.clear();
302 fZA.reserve(nmax);
303 fIsoXS.resize(nmax, 0.0);
304 }
305
306 void G4CrossSectionHP::DumpPhysicsTable(const G4ParticleDefinition&)
307 {
308 if (fManagerHP->GetVerboseLevel() == 0 || fPrinted)
309 return;
310
311 //
312 // Dump element based cross section
313 // range 10e-5 eV to 20 MeV
314 // 10 point per decade
315 // in barn
316 //
317 fPrinted = true;
318 G4cout << G4endl;
319 G4cout << "HP Cross Section " << fDataName << " for "
320 << fParticle->GetParticleName() << G4endl;
321 G4cout << "(Pointwise cross-section at 0 Kelvin.)" << G4endl;
322 G4cout << G4endl;
323 G4cout << "Name of Element" << G4endl;
324 G4cout << "Energy[eV] XS[barn]" << G4endl;
325 G4cout << G4endl;
326
327 const G4ElementTable* table = G4Element::GetElementTable();
328 for ( auto const & elm : *table ) {
329 G4int Z = elm->GetZasInt();
330 if (Z >= minZ && Z <= maxZ && nullptr != fData->GetElementData(Z - minZ)) {
331 G4cout << "---------------------------------------------------" << G4endl;
332 G4cout << elm->GetName() << G4endl;
333 std::size_t n = fData->GetNumberOfComponents(Z);
334 for (size_t i=0; i<n; ++i) {
335 G4cout << "## Z=" << Z << " A=" << fData->GetComponentID(Z - minZ, i);
336 G4cout << *(fData->GetComponentDataByIndex(Z - minZ, i)) << G4endl;
337 }
338 }
339 }
340 }
341
342 void G4CrossSectionHP::PrepareCache(const G4Material* mat)
343 {
344 fCurrentMat = mat;
345 fZA.clear();
346 for ( auto const & elm : *(mat->GetElementVector()) ) {
347 G4int Z = elm->GetZasInt();
348 for ( auto const & iso : *(elm->GetIsotopeVector()) ) {
349 G4int A = iso->GetN();
350 fZA.emplace_back(Z, A);
351 }
352 }
353 fIsoXS.resize(fZA.size(), 0.0);
354 }
355
356 void G4CrossSectionHP::Initialise(const G4int Z)
357 {
358 if (fManagerHP->GetVerboseLevel() > 1) {
359 G4cout << " G4CrossSectionHP::Initialise: Z=" << Z
360 << " for " << fDataName
361 << " minZ=" << minZ << " maxZ=" << maxZ << G4endl;
362 }
363 if (Z < minZ || Z > maxZ || nullptr != fData->GetElementData(Z - minZ)) {
364 return;
365 }
366 G4AutoLock l(&theHPXS);
367 if (nullptr != fData->GetElementData(Z - minZ)) {
368 l.unlock();
369 return;
370 }
371
372 // add empty vector to avoid double initialisation
373 fData->InitialiseForElement(Z - minZ, new G4PhysicsVector());
374
375 G4String tnam = "temp";
376 G4bool noComp = true;
377 for (G4int A=amin[Z]; A<=amax[Z]; ++A) {
378 std::ostringstream ost;
379 ost << fDataDirectory;
380 // first check special cases
381 if (6 == Z && 12 == A && fParticle == fNeutron) {
382 ost << Z << "_nat_" << elementName[Z];
383 } else if (18 == Z && 40 != A) {
384 continue;
385 } else if (27 == Z && 62 == A) {
386 ost << Z << "_62m1_" << elementName[Z];
387 } else if (47 == Z && 106 == A) {
388 ost << Z << "_106m1_" << elementName[Z];
389 } else if (48 == Z && 115 == A) {
390 ost << Z << "_115m1_" << elementName[Z];
391 } else if (52 == Z && 127 == A) {
392 ost << Z << "_127m1_" << elementName[Z];
393 } else if (52 == Z && 129 == A) {
394 ost << Z << "_129m1_" << elementName[Z];
395 } else if (52 == Z && 131 == A) {
396 ost << Z << "_131m1_" << elementName[Z];
397 } else if (61 == Z && 145 == A) {
398 ost << Z << "_147_" << elementName[Z];
399 } else if (67 == Z && 166 == A) {
400 ost << Z << "_166m1_" << elementName[Z];
401 } else if (73 == Z && 180 == A) {
402 ost << Z << "_180m1_" << elementName[Z];
403 } else if ((Z == 85 && A == 210) || (Z == 86 && A == 222) || (Z == 87 && A == 223)) {
404 ost << "84_209_" << elementName[84];
405 } else {
406 // the main file name
407 ost << Z << "_" << A << "_" << elementName[Z];
408 }
409 std::istringstream theXSData(tnam, std::ios::in);
410 fManagerHP->GetDataStream(ost.str().c_str(), theXSData);
411 if (theXSData) {
412 G4int i1, i2, n;
413 theXSData >> i1 >> i2 >> n;
414 if (fManagerHP->GetVerboseLevel() > 1) {
415 G4cout << "## G4CrossSectionHP::Initialise for Z=" << Z
416 << " A=" << A << " Npoints=" << n << G4endl;
417 }
418 G4double x, y;
419 G4PhysicsFreeVector* v = new G4PhysicsFreeVector(n);
420 for (G4int i=0; i<n; ++i) {
421 theXSData >> x >> y;
422 x *= CLHEP::eV;
423 y *= CLHEP::barn;
424 //G4cout << " e=" << x << " xs=" << y << G4endl;
425 v->PutValues((std::size_t)i, x, y);
426 }
427 v->EnableLogBinSearch(binSearch);
428 if (noComp) {
429 G4int nmax = amax[Z] - A + 1;
430 fData->InitialiseForComponent(Z - minZ, nmax);
431 noComp = false;
432 }
433 fData->AddComponent(Z - minZ, A, v);
434 }
435 }
436 if (noComp) { fData->InitialiseForComponent(Z - minZ, 0); }
437 l.unlock();
438 }
439
440