| 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: 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[] =
60 G4double G4ParticleInelasticXS::coeff[MAXZINEL
61 G4String G4ParticleInelasticXS::gDataDirectory
62
63 namespace
64 {
65 G4Mutex pInelasticXSMutex = G4MUTEX_INITIALI
66 G4String pname[5] = {"proton", "deuteron", "
67 }
68
69 G4ParticleInelasticXS::G4ParticleInelasticXS(c
70 : G4VCrossSectionDataSet("G4ParticleInelasti
71 particle(part),
72 elimit(20*CLHEP::MeV)
73 {
74 if (nullptr == part) {
75 G4Exception("G4ParticleInelasticXS::G4Part
76 FatalException, "NO particle definition in
77 } else {
78 verboseLevel = 0;
79 const G4String& particleName = particle->G
80 if(verboseLevel > 1) {
81 G4cout << "G4ParticleInelasticXS::G4Part
82 << particleName << " on atoms with Z <
83 }
84 auto xsr = G4CrossSectionDataSetRegistry::
85 if (particleName == "proton") {
86 highEnergyXsection = xsr->GetComponentCr
87 if(highEnergyXsection == nullptr) {
88 highEnergyXsection = new G4ComponentGGHadron
89 }
90 } else {
91 highEnergyXsection =
92 xsr->GetComponentCrossSection("Glauber
93 if(highEnergyXsection == nullptr) {
94 highEnergyXsection = new G4ComponentGGNuclNu
95 }
96 for (index=1; index<5; ++index) {
97 if (particleName == pname[index]) { br
98 }
99 index = std::min(index, 4);
100 if (1 < index) { SetMaxKinEnergy(25.6*CL
101 }
102 }
103 SetForAllAtomsAndEnergies(true);
104 if (gDataDirectory.empty()) {
105 gDataDirectory = G4HadronicParameters::Ins
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] + "PartI
112 }
113 }
114
115 void G4ParticleInelasticXS::CrossSectionDescri
116 {
117 outFile << "G4ParticleInelasticXS calculates
118 << "cross section on nuclei using da
119 << "neutron database. These data ar
120 << "the resonance region in order to
121 << "For high energy Glauber-Gribov c
122 }
123
124 G4bool
125 G4ParticleInelasticXS::IsElementApplicable(con
126 G4int, const G4Material*)
127 {
128 return true;
129 }
130
131 G4bool
132 G4ParticleInelasticXS::IsIsoApplicable(const G
133 G4int, G4int,
134 const G4Element*, const G4Mater
135 {
136 return true;
137 }
138
139 G4double
140 G4ParticleInelasticXS::GetElementCrossSection(
141
142 {
143 return ElementCrossSection(aParticle->GetKin
144 aParticle->GetLog
145 }
146
147 G4double
148 G4ParticleInelasticXS::ComputeCrossSectionPerE
149
150
151
152 {
153 return ElementCrossSection(ekin, loge, elm->
154 }
155
156 G4double G4ParticleInelasticXS::ElementCrossSe
157 {
158 G4int Z = (ZZ >= MAXZINELP) ? MAXZINELP - 1
159
160 // element data is always valid pointer by c
161 auto pv = GetPhysicsVector(Z);
162
163 // set to null x-section below lowest energy
164 G4double xs = 0.0;
165 if (ekin > pv->Energy(0)) {
166 xs = (ekin <= pv->GetMaxEnergy()) ? pv->Lo
167 : coeff[Z][index]*highEnergyXsection->GetI
168 ekin, Z, aeff[Z]);
169 }
170
171 #ifdef G4VERBOSE
172 if(verboseLevel > 1) {
173 G4cout << "ElmXS: Z= " << Z << " Ekin(MeV
174 << " xs(bn)= " << xs/CLHEP::barn << " el
175 << particle->GetParticleName()
176 << " idx= " << index << G4endl;
177 }
178 #endif
179 return xs;
180 }
181
182 G4double
183 G4ParticleInelasticXS::ComputeIsoCrossSection(
184
185
186
187 {
188 return IsoCrossSection(ekin, loge, Z, A);
189 }
190
191 G4double
192 G4ParticleInelasticXS::GetIsoCrossSection(cons
193 G4in
194 cons
195 {
196 return IsoCrossSection(aParticle->GetKinetic
197 aParticle->GetLogKine
198 }
199
200 G4double
201 G4ParticleInelasticXS::IsoCrossSection(G4doubl
202 G4int Z
203 {
204 G4double xs = 0.0;
205 G4int Z = (ZZ >= MAXZINELP) ? MAXZINELP - 1
206
207 // needed here to gurantee upload data for Z
208 auto pv = GetPhysicsVector(Z);
209
210 // compute isotope cross section if applicab
211 if (ekin <= elimit && data[index]->GetNumber
212 auto pviso = data[index]->GetComponentData
213 if (pviso != nullptr && ekin > pviso->Ener
214 xs = pviso->LogVectorValue(ekin, logE);
215 #ifdef G4VERBOSE
216 if(verboseLevel > 1) {
217 G4cout << "G4ParticleInelasticXS::IsoXS: for
218 << particle->GetParticleName()
219 << ekin/CLHEP::MeV << " xs(b)=
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->Lo
230 coeff[Z][index] *
231 highEnergyXsection->GetInelasticElementC
232 * A/aeff[Z];
233 }
234 #ifdef G4VERBOSE
235 if(verboseLevel > 1) {
236 G4cout << "IsoXS for " << particle->GetPa
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::Select
247 const G4Element* anElement, G4double kinE
248 {
249 G4int nIso = (G4int)anElement->GetNumberOfIs
250 const G4Isotope* iso = anElement->GetIsotope
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->Get
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]->GetN
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(kinE
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
297 {
298 if (verboseLevel > 0){
299 G4cout << "G4ParticleInelasticXS::BuildPhy
300 << p.GetParticleName() << G4endl;
301 }
302 if (&p != particle) {
303 G4ExceptionDescription ed;
304 ed << p.GetParticleName() << " is a wrong
305 << particle->GetParticleName() << " is
306 G4Exception("G4ParticleInelasticXS::BuildP
307 FatalException, ed, "");
308 return;
309 }
310
311 // it is possible re-initialisation for the
312 const G4ElementTable* table = G4Element::Get
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(), MAXZ
322 if ( nullptr == (data[index])->GetElementD
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])->GetElementDat
332
333 G4AutoLock l(&pInelasticXSMutex);
334 if ( nullptr == (data[index])->GetElementDat
335 // upload element data
336 std::ostringstream ost;
337 ost << gDataDirectory << "/" << pname[inde
338 G4PhysicsVector* v = RetrieveVector(ost, t
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
348 G4PhysicsVector* v1 = RetrieveVector(ost1, f
349 if (nullptr != v1) {
350 if (noComp) {
351 G4int nmax = amax[Z] - A + 1;
352 data[index]->InitialiseForComponent(Z, n
353 noComp = false;
354 }
355 data[index]->AddComponent(Z, A, v1);
356 }
357 }
358 }
359 // no components case
360 if (noComp) { data[index]->InitialiseForCo
361
362 // smooth transition
363 G4double sig1 = (*v)[v->GetVectorLength()
364 G4double ehigh = v->GetMaxEnergy();
365 G4double sig2 = highEnergyXsection->GetIne
366 particle, ehigh, Z, aeff[Z]);
367 coeff[Z][index] = (sig2 > 0.) ? sig1/sig2
368 }
369 l.unlock();
370 }
371
372 G4PhysicsVector*
373 G4ParticleInelasticXS::RetrieveVector(std::ost
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::Retr
384 FatalException, ed, "Check G4PARTICLEXSD
385 }
386 } else {
387 if(verboseLevel > 1) {
388 G4cout << "File " << ost.str()
389 << " is opened by G4ParticleInelasticXS
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::Retr
398 FatalException, ed, "Check G4PARTICLEXSD
399 }
400 }
401 return v;
402 }
403