| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // V. Ivanchenko September 2023
27 //
28 // G4CrossSectionHP is a generic class impleme
29 // cross sections for neutrons, protons and li
30 // It is an alternative to code developed by J
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 G4Par
62 const G4Str
63 const G4Str
64 G4int zmin,
65 : G4VCrossSectionDataSet(nameData),
66 fParticle(p),
67 fNeutron(G4Neutron::Neutron()),
68 fManagerHP(G4ParticleHPManager::GetInstanc
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::G4CrossSectio
81 << fDataName << " " << minZ << " < Z < "
82 << " EmaxT(MeV)=" << emaxT << G4endl;
83 G4cout << "Data directory: " << fDataDirec
84 }
85 auto ptr = G4ElementDataRegistry::Instance()
86 auto data = ptr->GetElementDataByName(fDataN
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
95 G4int, G4int,
96 const
97 const G4Material*)
98 {
99 return (dp->GetKineticEnergy() <= emax);
100 }
101
102 G4double G4CrossSectionHP::GetIsoCrossSection(
103
104
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
117 }
118
119 G4double
120 G4CrossSectionHP::ComputeIsoCrossSection(G4dou
121 const
122 G4int
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
136 }
137
138 G4double G4CrossSectionHP::IsoCrossSection(con
139 con
140 const G4int Z, const G4int A,
141 con
142 {
143 // G4cout << "G4CrossSectionHP::IsoCrossSect
144 // << " E(MeV)=" << ekin/MeV << " T=" << T <
145 G4double xs = 0.0;
146 if (ekin > emax || Z > maxZ || Z < minZ || e
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 d
155 if (nullptr == pv0) { return xs; }
156
157 const auto pv = fData->GetComponentDataByID(
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 |
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::
171 G4double sig = std::sqrt(2.0*e0/(3.0*massT
172
173 // projectile
174 G4LorentzVector lv(0., 0., std::sqrt(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, ma
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) {
197 }
198 if (ii > 0) { xs /= (G4double)(ii); }
199 }
200 #ifdef G4VERBOSE
201 if (verboseLevel > 1) {
202 G4cout << "G4CrossSectionHP::IsoXS " << fD
203 << " Z= " << Z << " A= " << A
204 << " Ekin(MeV)= " << ekin/MeV << " xs(b)
205 << " size=" << fZA.size() << G4end
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].secon
212 fIsoXS[i] = xs;
213 break;
214 }
215 }
216 return xs;
217 }
218
219 const G4Isotope* G4CrossSectionHP::SelectIsoto
220
221 {
222 std::size_t nIso = elm->GetNumberOfIsotopes(
223 const G4Isotope* iso = elm->GetIsotope(0);
224
225 //G4cout << "SelectIsotope NIso= " << nIso <
226 if(1 == nIso) { return iso; }
227
228 // more than 1 isotope
229 G4int Z = elm->GetZasInt();
230 if (Z >= minZ && Z <= maxZ && nullptr == fDa
231 Initialise(Z);
232 }
233
234 const G4double* abundVector = elm->GetRelati
235 G4double q = G4UniformRand();
236 G4double sum = 0.0;
237
238 // is there cached isotope wise cross sectio
239 std::size_t j;
240 if (Z < minZ || Z > maxZ || !CheckCache(Z) |
241 0 == fData->GetNumberOfComponents(Z - mi
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->GetIsotop
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
271 {
272 if (verboseLevel > 1) {
273 G4cout << "G4CrossSectionHP::BuildPhysicsT
274 << p.GetParticleName() << " and " << fDat
275 }
276
277 // it is possible re-initialisation for the
278 const G4ElementTable* table = G4Element::Get
279
280 // Access to elements
281 for ( auto const & elm : *table ) {
282 G4int Z = elm->GetZasInt();
283 if (Z >= minZ && Z <= maxZ && nullptr == f
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::GetMat
292 std::size_t n = 0;
293 for ( auto const & elm : *mat->GetElementV
294 std::size_t niso = elm->GetNumberOfIsoto
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
307 {
308 if (fManagerHP->GetVerboseLevel() == 0 || fP
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 <
320 << fParticle->GetParticleName() << G4endl;
321 G4cout << "(Pointwise cross-section at 0 Kel
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::Get
328 for ( auto const & elm : *table ) {
329 G4int Z = elm->GetZasInt();
330 if (Z >= minZ && Z <= maxZ && nullptr != f
331 G4cout << "-----------------------------
332 G4cout << elm->GetName() << G4endl;
333 std::size_t n = fData->GetNumberOfCompon
334 for (size_t i=0; i<n; ++i) {
335 G4cout << "## Z=" << Z << " A=" << fData->
336 G4cout << *(fData->GetComponentDataByIndex(Z
337 }
338 }
339 }
340 }
341
342 void G4CrossSectionHP::PrepareCache(const G4Ma
343 {
344 fCurrentMat = mat;
345 fZA.clear();
346 for ( auto const & elm : *(mat->GetElementVe
347 G4int Z = elm->GetZasInt();
348 for ( auto const & iso : *(elm->GetIsotope
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
357 {
358 if (fManagerHP->GetVerboseLevel() > 1) {
359 G4cout << " G4CrossSectionHP::Initialise:
360 << " for " << fDataName
361 << " minZ=" << minZ << " maxZ=" << maxZ <
362 }
363 if (Z < minZ || Z > maxZ || nullptr != fData
364 return;
365 }
366 G4AutoLock l(&theHPXS);
367 if (nullptr != fData->GetElementData(Z - min
368 l.unlock();
369 return;
370 }
371
372 // add empty vector to avoid double initiali
373 fData->InitialiseForElement(Z - minZ, new G4
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 == fNeu
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 ==
404 ost << "84_209_" << elementName[84];
405 } else {
406 // the main file name
407 ost << Z << "_" << A << "_" << elementNa
408 }
409 std::istringstream theXSData(tnam, std::io
410 fManagerHP->GetDataStream(ost.str().c_str(
411 if (theXSData) {
412 G4int i1, i2, n;
413 theXSData >> i1 >> i2 >> n;
414 if (fManagerHP->GetVerboseLevel() > 1) {
415 G4cout << "## G4CrossSectionHP::Initialise f
416 << " A=" << A << " Npoints=" << n <<
417 }
418 G4double x, y;
419 G4PhysicsFreeVector* v = new G4PhysicsFr
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 << G
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(
437 l.unlock();
438 }
439
440