| 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 // -------------------------------------------
27 //
28 // GEANT4 Class file
29 //
30 //
31 // File name: G4GammaNuclearXS
32 //
33 // Authors V.Ivantchenko, Geant4, 20 October
34 // B.Kutsenko, BINP/NSU, 10 August 20
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 = nullpt
63
64 G4double G4GammaNuclearXS::coeff[3][3];
65 G4double G4GammaNuclearXS::xs150[] = {0.0};
66 const G4int G4GammaNuclearXS::freeVectorExcept
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 betw
73 const G4double eTransitionBound = 150.*CLHEP
74 // A limit energy to correct CHIPS parameter
75 const G4double ehigh = 10*CLHEP::GeV;
76 }
77
78 G4GammaNuclearXS::G4GammaNuclearXS()
79 : G4VCrossSectionDataSet(Default_Name()), ga
80 {
81 verboseLevel = 0;
82 if (verboseLevel > 0) {
83 G4cout << "G4GammaNuclearXS::G4GammaNuclea
84 << MAXZGAMMAXS << G4endl;
85 }
86 ggXsection = dynamic_cast<G4PhotoNuclearCros
87 (G4CrossSectionDataSetRegistry::Instance()
88 if (ggXsection == nullptr) {
89 ggXsection = new G4PhotoNuclearCrossSectio
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
104 {
105 outFile << "G4GammaNuclearXS calculates the
106 << "cross-section for GDR energy reg
107 << "data from the high precision\n"
108 << "IAEA photonuclear database (2019
109 << "implemented with previous CHIPS photon
110 }
111
112 G4bool
113 G4GammaNuclearXS::IsElementApplicable(const G4
114 G4int, c
115 {
116 return true;
117 }
118
119 G4bool G4GammaNuclearXS::IsIsoApplicable(const
120 G4int
121 const
122 {
123 return true;
124 }
125
126 G4double
127 G4GammaNuclearXS::GetElementCrossSection(const
128 G4int
129 {
130 return ElementCrossSection(aParticle->GetKin
131 }
132
133 G4double
134 G4GammaNuclearXS::ElementCrossSection(const G4
135 {
136 // check cache
137 const G4int Z = (ZZ < MAXZGAMMAXS) ? ZZ : MA
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
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(e
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(e
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)/(eTr
167 }
168
169 #ifdef G4VERBOSE
170 if(verboseLevel > 1) {
171 G4cout << "Z= " << Z << " Ekin(MeV)= " <<
172 << ", nElmXS(b)= " << fXS/CLHEP::barn
173 << G4endl;
174 }
175 #endif
176 return fXS;
177 }
178
179 G4double G4GammaNuclearXS::LowEnergyCrossSecti
180 {
181 const G4int Z = (ZZ < MAXZGAMMAXS) ? ZZ : MA
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 G
190 {
191 return IsoCrossSection(aParticle->GetKinetic
192 }
193
194 G4double
195 G4GammaNuclearXS::IsoCrossSection(const G4doub
196 {
197 const G4int Z = (ZZ < MAXZGAMMAXS) ? ZZ : MA
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
209 }
210 }
211
212 #ifdef G4VERBOSE
213 if(verboseLevel > 1) {
214 G4cout << "G4GammaNuclearXS::IsoXS: Z= "
215 << " Ekin(MeV)= " << ekin/CLHEP::MeV
216 << ", ElmXS(b)= " << xs/CLHEP::barn << G
217 }
218 #endif
219 return xs;
220 }
221
222 const G4Isotope* G4GammaNuclearXS::SelectIsoto
223 const G4Element* anElement, G4double ki
224 {
225 std::size_t nIso = anElement->GetNumberOfIso
226 const G4Isotope* iso = anElement->GetIsotope
227
228 if(1 == nIso) { return iso; }
229
230 const G4double* abundVector = anElement->Get
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 " << (*isoVe
240 // << " abund= " << abundVector[j]
241 sum += abundVector[j]*
242 IsoCrossSection(kinEnergy, Z, anElement-
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
256 {
257 if(verboseLevel > 1) {
258 G4cout << "G4GammaNuclearXS::BuildPhysicsT
259 << p.GetParticleName() << G4endl;
260 }
261 if(p.GetParticleName() != "gamma") {
262 G4ExceptionDescription ed;
263 ed << p.GetParticleName() << " is a wrong
264 << " only gamma is allowed";
265 G4Exception("G4GammaNuclearXS::BuildPhysic
266 FatalException, ed, "");
267 return;
268 }
269
270 // prepare isotope selection
271 const G4ElementTable* table = G4Element::Get
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::FindDirector
281 {
282 // build the complete string identifying the
283 if(gDataDirectory.empty()) {
284 std::ostringstream ost;
285 ost << G4HadronicParameters::Instance()->G
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, tru
297
298 data->InitialiseForElement(Z, v);
299 /*
300 G4cout << "G4GammaNuclearXS::Initialise for
301 << " A= " << Amean << " Amin= " << amin[Z]
302 << " Amax= " << amax[Z] << G4endl;
303 */
304 xs150[Z] = ggXsection->ComputeElementXSectio
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->ComputeIsoXSec
313 G4double sig2 = ggXsection->ComputeElement
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::ostrings
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::RetrieveV
334 FatalException, ed, "Check G4PARTICLEXSD
335 }
336 } else {
337 if(verboseLevel > 1) {
338 G4cout << "File " << ost.str()
339 << " is opened by G4GammaNuclearXS" <<
340 }
341 // retrieve data from DB
342 if(std::find(std::begin(freeVectorExceptio
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::RetrieveV
353 FatalException, ed, "Check G4PARTICLEXSD
354 }
355 }
356 return v;
357 }
358
359