| 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: 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 "G4ElementTable.hh"
42 #include "G4Material.hh"
43 #include "G4Element.hh"
44 #include "G4PhysicsLogVector.hh"
45 #include "G4CrossSectionDataSetRegistry.hh"
46 #include "G4ComponentGGHadronNucleusXsc.hh"
47 #include "G4ComponentGGNuclNuclXsc.hh"
48 #include "G4HadronicParameters.hh"
49 #include "G4Proton.hh"
50 #include "Randomize.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4IsotopeList.hh"
53 #include "G4HadronicParameters.hh"
54 #include "G4AutoLock.hh"
55
56 #include <fstream>
57 #include <sstream>
58
59 G4ElementData* G4ParticleInelasticXS::data[] = {nullptr, nullptr, nullptr, nullptr, nullptr};
60 G4double G4ParticleInelasticXS::coeff[MAXZINELP][5] = {{1.0}, {1.0}, {1.0}, {1.0}, {1.0}};
61 G4String G4ParticleInelasticXS::gDataDirectory = {""};
62
63 namespace
64 {
65 G4Mutex pInelasticXSMutex = G4MUTEX_INITIALIZER;
66 G4String pname[5] = {"proton", "deuteron", "triton", "He3", "alpha"};
67 }
68
69 G4ParticleInelasticXS::G4ParticleInelasticXS(const G4ParticleDefinition* part)
70 : G4VCrossSectionDataSet("G4ParticleInelasticXS"),
71 particle(part),
72 elimit(20*CLHEP::MeV)
73 {
74 if (nullptr == part) {
75 G4Exception("G4ParticleInelasticXS::G4ParticleInelasticXS(..)","had015",
76 FatalException, "NO particle definition in constructor");
77 } else {
78 verboseLevel = 0;
79 const G4String& particleName = particle->GetParticleName();
80 if(verboseLevel > 1) {
81 G4cout << "G4ParticleInelasticXS::G4ParticleInelasticXS for "
82 << particleName << " on atoms with Z < " << MAXZINELP << G4endl;
83 }
84 auto xsr = G4CrossSectionDataSetRegistry::Instance();
85 if (particleName == "proton") {
86 highEnergyXsection = xsr->GetComponentCrossSection("Glauber-Gribov");
87 if(highEnergyXsection == nullptr) {
88 highEnergyXsection = new G4ComponentGGHadronNucleusXsc();
89 }
90 } else {
91 highEnergyXsection =
92 xsr->GetComponentCrossSection("Glauber-Gribov Nucl-nucl");
93 if(highEnergyXsection == nullptr) {
94 highEnergyXsection = new G4ComponentGGNuclNuclXsc();
95 }
96 for (index=1; index<5; ++index) {
97 if (particleName == pname[index]) { break; }
98 }
99 index = std::min(index, 4);
100 if (1 < index) { SetMaxKinEnergy(25.6*CLHEP::PeV); }
101 }
102 }
103 SetForAllAtomsAndEnergies(true);
104 if (gDataDirectory.empty()) {
105 gDataDirectory = G4HadronicParameters::Instance()->GetDirPARTICLEXS();
106 }
107 G4String ss = pname[index] + "ParticleXS";
108 SetName(ss);
109 if (data[index] == nullptr) {
110 data[index] = new G4ElementData(MAXZINELP);
111 data[index]->SetName(pname[index] + "PartInel");
112 }
113 }
114
115 void G4ParticleInelasticXS::CrossSectionDescription(std::ostream& outFile) const
116 {
117 outFile << "G4ParticleInelasticXS calculates n, p, d, t, he3, he4 inelastic\n"
118 << "cross section on nuclei using data from the high precision\n"
119 << "neutron database. These data are simplified and smoothed over\n"
120 << "the resonance region in order to reduce CPU time.\n"
121 << "For high energy Glauber-Gribov cross section model is used.\n";
122 }
123
124 G4bool
125 G4ParticleInelasticXS::IsElementApplicable(const G4DynamicParticle*,
126 G4int, const G4Material*)
127 {
128 return true;
129 }
130
131 G4bool
132 G4ParticleInelasticXS::IsIsoApplicable(const G4DynamicParticle*,
133 G4int, G4int,
134 const G4Element*, const G4Material*)
135 {
136 return true;
137 }
138
139 G4double
140 G4ParticleInelasticXS::GetElementCrossSection(const G4DynamicParticle* aParticle,
141 G4int Z, const G4Material*)
142 {
143 return ElementCrossSection(aParticle->GetKineticEnergy(),
144 aParticle->GetLogKineticEnergy(), Z);
145 }
146
147 G4double
148 G4ParticleInelasticXS::ComputeCrossSectionPerElement(G4double ekin, G4double loge,
149 const G4ParticleDefinition*,
150 const G4Element* elm,
151 const G4Material*)
152 {
153 return ElementCrossSection(ekin, loge, elm->GetZasInt());
154 }
155
156 G4double G4ParticleInelasticXS::ElementCrossSection(G4double ekin, G4double loge, G4int ZZ)
157 {
158 G4int Z = (ZZ >= MAXZINELP) ? MAXZINELP - 1 : ZZ;
159
160 // element data is always valid pointer by construction of XS
161 auto pv = GetPhysicsVector(Z);
162
163 // set to null x-section below lowest energy in the table
164 G4double xs = 0.0;
165 if (ekin > pv->Energy(0)) {
166 xs = (ekin <= pv->GetMaxEnergy()) ? pv->LogVectorValue(ekin, loge)
167 : coeff[Z][index]*highEnergyXsection->GetInelasticElementCrossSection(particle,
168 ekin, Z, aeff[Z]);
169 }
170
171 #ifdef G4VERBOSE
172 if(verboseLevel > 1) {
173 G4cout << "ElmXS: Z= " << Z << " Ekin(MeV)= " << ekin/CLHEP::MeV
174 << " xs(bn)= " << xs/CLHEP::barn << " element data for "
175 << particle->GetParticleName()
176 << " idx= " << index << G4endl;
177 }
178 #endif
179 return xs;
180 }
181
182 G4double
183 G4ParticleInelasticXS::ComputeIsoCrossSection(G4double ekin, G4double loge,
184 const G4ParticleDefinition*,
185 G4int Z, G4int A, const G4Isotope*,
186 const G4Element*, const G4Material*)
187 {
188 return IsoCrossSection(ekin, loge, Z, A);
189 }
190
191 G4double
192 G4ParticleInelasticXS::GetIsoCrossSection(const G4DynamicParticle* aParticle,
193 G4int Z, G4int A, const G4Isotope*,
194 const G4Element*, const G4Material*)
195 {
196 return IsoCrossSection(aParticle->GetKineticEnergy(),
197 aParticle->GetLogKineticEnergy(), Z, A);
198 }
199
200 G4double
201 G4ParticleInelasticXS::IsoCrossSection(G4double ekin, G4double logE,
202 G4int ZZ, G4int A)
203 {
204 G4double xs = 0.0;
205 G4int Z = (ZZ >= MAXZINELP) ? MAXZINELP - 1 : ZZ;
206
207 // needed here to gurantee upload data for Z
208 auto pv = GetPhysicsVector(Z);
209
210 // compute isotope cross section if applicable
211 if (ekin <= elimit && data[index]->GetNumberOfComponents(Z) > 0) {
212 auto pviso = data[index]->GetComponentDataByID(Z, A);
213 if (pviso != nullptr && ekin > pviso->Energy(0)) {
214 xs = pviso->LogVectorValue(ekin, logE);
215 #ifdef G4VERBOSE
216 if(verboseLevel > 1) {
217 G4cout << "G4ParticleInelasticXS::IsoXS: for "
218 << particle->GetParticleName() << " Ekin(MeV)= "
219 << ekin/CLHEP::MeV << " xs(b)= " << xs/CLHEP::barn
220 << " Z= " << Z << " A= " << A
221 << " idx= " << index << G4endl;
222 }
223 #endif
224 return xs;
225 }
226 }
227 // use element x-section
228 if (ekin > pv->Energy(0)) {
229 xs = (ekin <= pv->GetMaxEnergy()) ? pv->LogVectorValue(ekin, logE) :
230 coeff[Z][index] *
231 highEnergyXsection->GetInelasticElementCrossSection(particle, ekin, Z, aeff[Z])
232 * A/aeff[Z];
233 }
234 #ifdef G4VERBOSE
235 if(verboseLevel > 1) {
236 G4cout << "IsoXS for " << particle->GetParticleName()
237 << " Target Z= " << Z << " A= " << A
238 << " Ekin(MeV)= " << ekin/CLHEP::MeV
239 << " xs(bn)= " << xs/CLHEP::barn
240 << " idx= " << index << G4endl;
241 }
242 #endif
243 return xs;
244 }
245
246 const G4Isotope* G4ParticleInelasticXS::SelectIsotope(
247 const G4Element* anElement, G4double kinEnergy, G4double logE)
248 {
249 G4int nIso = (G4int)anElement->GetNumberOfIsotopes();
250 const G4Isotope* iso = anElement->GetIsotope(0);
251
252 if (1 == nIso) { return iso; }
253
254 // more than 1 isotope
255 G4int Z = anElement->GetZasInt();
256
257 // initialisation for given Z
258 GetPhysicsVector(Z);
259
260 const G4double* abundVector = anElement->GetRelativeAbundanceVector();
261 G4double q = G4UniformRand();
262 G4double sum = 0.0;
263 G4int j;
264
265 // isotope wise cross section not available
266 if (Z >= MAXZINELP || 0 == data[index]->GetNumberOfComponents(Z)) {
267 for (j=0; j<nIso; ++j) {
268 sum += abundVector[j];
269 if(q <= sum) {
270 iso = anElement->GetIsotope(j);
271 break;
272 }
273 }
274 return iso;
275 }
276
277 G4int nn = (G4int)temp.size();
278 if (nn < nIso) { temp.resize(nIso, 0.); }
279
280 for (j=0; j<nIso; ++j) {
281 sum += abundVector[j]*IsoCrossSection(kinEnergy, logE, Z,
282 anElement->GetIsotope(j)->GetN());
283 temp[j] = sum;
284 }
285 sum *= q;
286 for (j=0; j<nIso; ++j) {
287 if (temp[j] >= sum) {
288 iso = anElement->GetIsotope(j);
289 break;
290 }
291 }
292 return iso;
293 }
294
295 void
296 G4ParticleInelasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
297 {
298 if (verboseLevel > 0){
299 G4cout << "G4ParticleInelasticXS::BuildPhysicsTable for "
300 << p.GetParticleName() << G4endl;
301 }
302 if (&p != particle) {
303 G4ExceptionDescription ed;
304 ed << p.GetParticleName() << " is a wrong particle type -"
305 << particle->GetParticleName() << " is expected";
306 G4Exception("G4ParticleInelasticXS::BuildPhysicsTable(..)","had012",
307 FatalException, ed, "");
308 return;
309 }
310
311 // it is possible re-initialisation for the new run
312 const G4ElementTable* table = G4Element::GetElementTable();
313
314 // prepare isotope selection
315 std::size_t nIso = temp.size();
316
317 // Access to elements
318 for ( auto const & elm : *table ) {
319 std::size_t n = elm->GetNumberOfIsotopes();
320 if (n > nIso) { nIso = n; }
321 G4int Z = std::min( elm->GetZasInt(), MAXZINELP-1);
322 if ( nullptr == (data[index])->GetElementData(Z) ) {
323 Initialise(Z);
324 }
325 }
326 temp.resize(nIso, 0.0);
327 }
328
329 void G4ParticleInelasticXS::Initialise(G4int Z)
330 {
331 if ( nullptr != (data[index])->GetElementData(Z) ) { return; }
332
333 G4AutoLock l(&pInelasticXSMutex);
334 if ( nullptr == (data[index])->GetElementData(Z) ) {
335 // upload element data
336 std::ostringstream ost;
337 ost << gDataDirectory << "/" << pname[index] << "/inel" << Z;
338 G4PhysicsVector* v = RetrieveVector(ost, true);
339 data[index]->InitialiseForElement(Z, v);
340
341 // upload isotope data
342 G4bool noComp = true;
343 if (amin[Z] < amax[Z]) {
344
345 for (G4int A=amin[Z]; A<=amax[Z]; ++A) {
346 std::ostringstream ost1;
347 ost1 << gDataDirectory << "/" << pname[index] << "/inel" << Z << "_" << A;
348 G4PhysicsVector* v1 = RetrieveVector(ost1, false);
349 if (nullptr != v1) {
350 if (noComp) {
351 G4int nmax = amax[Z] - A + 1;
352 data[index]->InitialiseForComponent(Z, nmax);
353 noComp = false;
354 }
355 data[index]->AddComponent(Z, A, v1);
356 }
357 }
358 }
359 // no components case
360 if (noComp) { data[index]->InitialiseForComponent(Z, 0); }
361
362 // smooth transition
363 G4double sig1 = (*v)[v->GetVectorLength()-1];
364 G4double ehigh = v->GetMaxEnergy();
365 G4double sig2 = highEnergyXsection->GetInelasticElementCrossSection(
366 particle, ehigh, Z, aeff[Z]);
367 coeff[Z][index] = (sig2 > 0.) ? sig1/sig2 : 1.0;
368 }
369 l.unlock();
370 }
371
372 G4PhysicsVector*
373 G4ParticleInelasticXS::RetrieveVector(std::ostringstream& ost, G4bool warn)
374 {
375 G4PhysicsLogVector* v = nullptr;
376 std::ifstream filein(ost.str().c_str());
377 if (!filein.is_open()) {
378 if (warn) {
379 G4ExceptionDescription ed;
380 ed << "Data file <" << ost.str().c_str()
381 << "> is not opened! index=" << index
382 << " dir: <" << gDataDirectory << ">. ";
383 G4Exception("G4ParticleInelasticXS::RetrieveVector(..)","had014",
384 FatalException, ed, "Check G4PARTICLEXSDATA");
385 }
386 } else {
387 if(verboseLevel > 1) {
388 G4cout << "File " << ost.str()
389 << " is opened by G4ParticleInelasticXS" << G4endl;
390 }
391 // retrieve data from DB
392 v = new G4PhysicsLogVector();
393 if(!v->Retrieve(filein, true)) {
394 G4ExceptionDescription ed;
395 ed << "Data file <" << ost.str().c_str()
396 << "> is not retrieved!";
397 G4Exception("G4ParticleInelasticXS::RetrieveVector(..)","had015",
398 FatalException, ed, "Check G4PARTICLEXSDATA");
399 }
400 }
401 return v;
402 }
403