| 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 // -------------------------------------------------------------------
27 //
28 // GEANT4 Class file
29 //
30 //
31 // File name: G4GammaNuclearXS
32 //
33 // Authors V.Ivantchenko, Geant4, 20 October 2020
34 // B.Kutsenko, BINP/NSU, 10 August 2021
35 //
36 // Modifications:
37 //
38
39 #include "G4GammaNuclearXS.hh"
40 #include "G4DynamicParticle.hh"
41 #include "G4ThreeVector.hh"
42 #include "G4ElementTable.hh"
43 #include "G4Material.hh"
44 #include "G4Element.hh"
45 #include "G4ElementData.hh"
46 #include "G4PhysicsVector.hh"
47 #include "G4PhysicsLinearVector.hh"
48 #include "G4PhysicsFreeVector.hh"
49 #include "G4CrossSectionDataSetRegistry.hh"
50 #include "G4PhotoNuclearCrossSection.hh"
51 #include "G4HadronicParameters.hh"
52 #include "Randomize.hh"
53 #include "G4SystemOfUnits.hh"
54 #include "G4Gamma.hh"
55 #include "G4IsotopeList.hh"
56 #include "G4AutoLock.hh"
57
58 #include <fstream>
59 #include <sstream>
60 #include <vector>
61
62 G4ElementData* G4GammaNuclearXS::data = nullptr;
63
64 G4double G4GammaNuclearXS::coeff[3][3];
65 G4double G4GammaNuclearXS::xs150[] = {0.0};
66 const G4int G4GammaNuclearXS::freeVectorException[] = {
67 4, 6, 7, 8, 27, 39, 45, 65, 67, 69, 73};
68 G4String G4GammaNuclearXS::gDataDirectory = "";
69
70 namespace
71 {
72 // Upper limit of the linear transition between IAEA database and CHIPS model
73 const G4double eTransitionBound = 150.*CLHEP::MeV;
74 // A limit energy to correct CHIPS parameterisation for light isotopes
75 const G4double ehigh = 10*CLHEP::GeV;
76 }
77
78 G4GammaNuclearXS::G4GammaNuclearXS()
79 : G4VCrossSectionDataSet(Default_Name()), gamma(G4Gamma::Gamma())
80 {
81 verboseLevel = 0;
82 if (verboseLevel > 0) {
83 G4cout << "G4GammaNuclearXS::G4GammaNuclearXS Initialise for Z < "
84 << MAXZGAMMAXS << G4endl;
85 }
86 ggXsection = dynamic_cast<G4PhotoNuclearCrossSection*>
87 (G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet("PhotoNuclearXS"));
88 if (ggXsection == nullptr) {
89 ggXsection = new G4PhotoNuclearCrossSection();
90 }
91 SetForAllAtomsAndEnergies(true);
92
93 // full data set is uploaded once
94 if (nullptr == data) {
95 data = new G4ElementData(MAXZGAMMAXS);
96 data->SetName("gNuclear");
97 for (G4int Z=1; Z<MAXZGAMMAXS; ++Z) {
98 Initialise(Z);
99 }
100 }
101 }
102
103 void G4GammaNuclearXS::CrossSectionDescription(std::ostream& outFile) const
104 {
105 outFile << "G4GammaNuclearXS calculates the gamma nuclear\n"
106 << "cross-section for GDR energy region on nuclei using "
107 << "data from the high precision\n"
108 << "IAEA photonuclear database (2019). Then liniear connection\n"
109 << "implemented with previous CHIPS photonuclear model." << G4endl;
110 }
111
112 G4bool
113 G4GammaNuclearXS::IsElementApplicable(const G4DynamicParticle*,
114 G4int, const G4Material*)
115 {
116 return true;
117 }
118
119 G4bool G4GammaNuclearXS::IsIsoApplicable(const G4DynamicParticle*,
120 G4int, G4int,
121 const G4Element*, const G4Material*)
122 {
123 return true;
124 }
125
126 G4double
127 G4GammaNuclearXS::GetElementCrossSection(const G4DynamicParticle* aParticle,
128 G4int Z, const G4Material*)
129 {
130 return ElementCrossSection(aParticle->GetKineticEnergy(), Z);
131 }
132
133 G4double
134 G4GammaNuclearXS::ElementCrossSection(const G4double ekin, const G4int ZZ)
135 {
136 // check cache
137 const G4int Z = (ZZ < MAXZGAMMAXS) ? ZZ : MAXZGAMMAXS - 1;
138 if(Z == fZ && ekin == fEkin) { return fXS; }
139 fZ = Z;
140 fEkin = ekin;
141
142 auto pv = data->GetElementData(Z);
143 const G4double limCHIPS1 = 25*CLHEP::MeV;
144 const G4double limCHIPS2 = 16*CLHEP::MeV;
145 if (pv == nullptr || 1 == Z || Z == 40 || Z == 74 ||
146 (Z == 24 && ekin >= limCHIPS1) ||
147 (Z == 39 && ekin >= limCHIPS1) ||
148 (Z == 50 && ekin >= limCHIPS2) ||
149 (Z == 64 && ekin >= limCHIPS2)
150 ) {
151 fXS = ggXsection->ComputeElementXSection(ekin, Z);
152 return fXS;
153 }
154 const G4double emax = pv->GetMaxEnergy();
155
156 // low energy based on data
157 if(ekin <= emax) {
158 fXS = pv->Value(ekin);
159 // high energy CHIPS parameterisation
160 } else if(ekin >= eTransitionBound) {
161 fXS = ggXsection->ComputeElementXSection(ekin, Z);
162 // linear interpolation
163 } else {
164 const G4double rxs = xs150[Z];
165 const G4double lxs = pv->Value(emax);
166 fXS = lxs + (ekin - emax)*(rxs - lxs)/(eTransitionBound - emax);
167 }
168
169 #ifdef G4VERBOSE
170 if(verboseLevel > 1) {
171 G4cout << "Z= " << Z << " Ekin(MeV)= " << ekin/CLHEP::MeV
172 << ", nElmXS(b)= " << fXS/CLHEP::barn
173 << G4endl;
174 }
175 #endif
176 return fXS;
177 }
178
179 G4double G4GammaNuclearXS::LowEnergyCrossSection(G4double ekin, G4int ZZ)
180 {
181 const G4int Z = (ZZ < MAXZGAMMAXS) ? ZZ : MAXZGAMMAXS - 1;
182 auto pv = data->GetElementData(Z);
183 return pv->Value(ekin);
184 }
185
186 G4double G4GammaNuclearXS::GetIsoCrossSection(
187 const G4DynamicParticle* aParticle,
188 G4int Z, G4int A,
189 const G4Isotope*, const G4Element*, const G4Material*)
190 {
191 return IsoCrossSection(aParticle->GetKineticEnergy(), Z, A);
192 }
193
194 G4double
195 G4GammaNuclearXS::IsoCrossSection(const G4double ekin, const G4int ZZ, const G4int A)
196 {
197 const G4int Z = (ZZ < MAXZGAMMAXS) ? ZZ : MAXZGAMMAXS - 1;
198 // cross section per element
199 G4double xs = ElementCrossSection(ekin, Z);
200
201 if (Z > 2) {
202 xs *= A/aeff[Z];
203 } else {
204 G4int AA = A - amin[Z];
205 if(ekin >= ehigh && AA >=0 && AA <=2) {
206 xs *= coeff[Z][AA];
207 } else {
208 xs = ggXsection->ComputeIsoXSection(ekin, Z, A);
209 }
210 }
211
212 #ifdef G4VERBOSE
213 if(verboseLevel > 1) {
214 G4cout << "G4GammaNuclearXS::IsoXS: Z= " << Z << " A= " << A
215 << " Ekin(MeV)= " << ekin/CLHEP::MeV
216 << ", ElmXS(b)= " << xs/CLHEP::barn << G4endl;
217 }
218 #endif
219 return xs;
220 }
221
222 const G4Isotope* G4GammaNuclearXS::SelectIsotope(
223 const G4Element* anElement, G4double kinEnergy, G4double)
224 {
225 std::size_t nIso = anElement->GetNumberOfIsotopes();
226 const G4Isotope* iso = anElement->GetIsotope(0);
227
228 if(1 == nIso) { return iso; }
229
230 const G4double* abundVector = anElement->GetRelativeAbundanceVector();
231 G4double sum = 0.0;
232 G4int Z = anElement->GetZasInt();
233
234 // use isotope cross sections
235 std::size_t nn = temp.size();
236 if(nn < nIso) { temp.resize(nIso, 0.); }
237
238 for (std::size_t j=0; j<nIso; ++j) {
239 //G4cout << j << "-th isotope " << (*isoVector)[j]->GetN()
240 // << " abund= " << abundVector[j] << G4endl;
241 sum += abundVector[j]*
242 IsoCrossSection(kinEnergy, Z, anElement->GetIsotope((G4int)j)->GetN());
243 temp[j] = sum;
244 }
245 sum *= G4UniformRand();
246 for (std::size_t j = 0; j<nIso; ++j) {
247 if(temp[j] >= sum) {
248 iso = anElement->GetIsotope((G4int)j);
249 break;
250 }
251 }
252 return iso;
253 }
254
255 void G4GammaNuclearXS::BuildPhysicsTable(const G4ParticleDefinition& p)
256 {
257 if(verboseLevel > 1) {
258 G4cout << "G4GammaNuclearXS::BuildPhysicsTable for "
259 << p.GetParticleName() << G4endl;
260 }
261 if(p.GetParticleName() != "gamma") {
262 G4ExceptionDescription ed;
263 ed << p.GetParticleName() << " is a wrong particle type -"
264 << " only gamma is allowed";
265 G4Exception("G4GammaNuclearXS::BuildPhysicsTable(..)","had012",
266 FatalException, ed, "");
267 return;
268 }
269
270 // prepare isotope selection
271 const G4ElementTable* table = G4Element::GetElementTable();
272 std::size_t nIso = temp.size();
273 for (auto const & elm : *table ) {
274 std::size_t n = elm->GetNumberOfIsotopes();
275 if (n > nIso) { nIso = n; }
276 }
277 temp.resize(nIso, 0.0);
278 }
279
280 const G4String& G4GammaNuclearXS::FindDirectoryPath()
281 {
282 // build the complete string identifying the file with the data set
283 if(gDataDirectory.empty()) {
284 std::ostringstream ost;
285 ost << G4HadronicParameters::Instance()->GetDirPARTICLEXS() << "/gamma/inel";
286 gDataDirectory = ost.str();
287 }
288 return gDataDirectory;
289 }
290
291 void G4GammaNuclearXS::Initialise(G4int Z)
292 {
293 // upload data from file
294 std::ostringstream ost;
295 ost << FindDirectoryPath() << Z ;
296 G4PhysicsVector* v = RetrieveVector(ost, true, Z);
297
298 data->InitialiseForElement(Z, v);
299 /*
300 G4cout << "G4GammaNuclearXS::Initialise for Z= " << Z
301 << " A= " << Amean << " Amin= " << amin[Z]
302 << " Amax= " << amax[Z] << G4endl;
303 */
304 xs150[Z] = ggXsection->ComputeElementXSection(eTransitionBound, Z);
305
306 // compute corrections for low Z data
307 if(Z <= 2){
308 if(amax[Z] > amin[Z]) {
309 for(G4int A=amin[Z]; A<=amax[Z]; ++A) {
310 G4int AA = A - amin[Z];
311 if(AA >= 0 && AA <= 2) {
312 G4double sig1 = ggXsection->ComputeIsoXSection(ehigh, Z, A);
313 G4double sig2 = ggXsection->ComputeElementXSection(ehigh, Z);
314 if(sig2 > 0.) { coeff[Z][AA] = (sig1/sig2); }
315 else { coeff[Z][AA] = 1.; }
316 }
317 }
318 }
319 }
320 }
321
322 G4PhysicsVector*
323 G4GammaNuclearXS::RetrieveVector(std::ostringstream& ost, G4bool warn, G4int Z)
324 {
325 G4PhysicsVector* v = nullptr;
326
327 std::ifstream filein(ost.str().c_str());
328 if (!filein.is_open()) {
329 if(warn) {
330 G4ExceptionDescription ed;
331 ed << "Data file <" << ost.str().c_str()
332 << "> is not opened!";
333 G4Exception("G4GammaNuclearXS::RetrieveVector(..)","had014",
334 FatalException, ed, "Check G4PARTICLEXSDATA");
335 }
336 } else {
337 if(verboseLevel > 1) {
338 G4cout << "File " << ost.str()
339 << " is opened by G4GammaNuclearXS" << G4endl;
340 }
341 // retrieve data from DB
342 if(std::find(std::begin(freeVectorException), std::end(freeVectorException), Z)
343 == std::end(freeVectorException)) {
344 v = new G4PhysicsLinearVector(false);
345 } else {
346 v = new G4PhysicsFreeVector(false);
347 }
348 if(!v->Retrieve(filein, true)) {
349 G4ExceptionDescription ed;
350 ed << "Data file <" << ost.str().c_str()
351 << "> is not retrieved!";
352 G4Exception("G4GammaNuclearXS::RetrieveVector(..)","had015",
353 FatalException, ed, "Check G4PARTICLEXSDATA");
354 }
355 }
356 return v;
357 }
358
359