| 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: G4NeutronElasticXS
32 //
33 // Author Ivantchenko, Geant4, 3-Aug-09
34 //
35 // Modifications:
36 //
37
38 #include "G4NeutronElasticXS.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 "G4HadronicParameters.hh"
48 #include "Randomize.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "G4IsotopeList.hh"
51 #include "G4AutoLock.hh"
52
53 #include <fstream>
54 #include <sstream>
55
56 G4PhysicsVector* G4NeutronElasticXS::data[] =
57 G4double G4NeutronElasticXS::coeff[] = {0.0};
58 G4String G4NeutronElasticXS::gDataDirectory =
59 G4bool G4NeutronElasticXS::fLock = true;
60
61 namespace
62 {
63 G4Mutex nElasticXSMutex = G4MUTEX_INITIALIZE
64 }
65
66 G4NeutronElasticXS::G4NeutronElasticXS()
67 : G4VCrossSectionDataSet(Default_Name()),
68 neutron(G4Neutron::Neutron())
69 {
70 // verboseLevel = 0;
71 if (verboseLevel > 0){
72 G4cout << "G4NeutronElasticXS::G4NeutronE
73 << MAXZEL << G4endl;
74 }
75 ggXsection =
76 G4CrossSectionDataSetRegistry::Instance()-
77 if (ggXsection == nullptr)
78 ggXsection = new G4ComponentGGHadronNucleu
79 SetForAllAtomsAndEnergies(true);
80 FindDirectoryPath();
81 }
82
83 G4NeutronElasticXS::~G4NeutronElasticXS()
84 {
85 if (isFirst) {
86 for(G4int i=0; i<MAXZEL; ++i) {
87 delete data[i];
88 data[i] = nullptr;
89 }
90 }
91 }
92
93 void G4NeutronElasticXS::CrossSectionDescripti
94 {
95 outFile << "G4NeutronElasticXS calculates th
96 << "cross section on nuclei using da
97 << "neutron database. These data ar
98 << "the resonance region in order to
99 << "For high energies Glauber-Gribiv
100 }
101
102 G4bool
103 G4NeutronElasticXS::IsElementApplicable(const
104 G4int, const G4Material*)
105 {
106 return true;
107 }
108
109 G4bool G4NeutronElasticXS::IsIsoApplicable(con
110 G4i
111 con
112 {
113 return false;
114 }
115
116 G4double
117 G4NeutronElasticXS::GetElementCrossSection(con
118 G4int Z, const G4Material*)
119 {
120 return ElementCrossSection(aParticle->GetKin
121 aParticle->GetLogKineticEnergy(), Z
122 }
123
124 G4double
125 G4NeutronElasticXS::ComputeCrossSectionPerElem
126 const G4ParticleDefinition*,
127 const G4Element* elm,
128 const G4Material*)
129 {
130 return ElementCrossSection(ekin, loge, elm->
131 }
132
133 G4double G4NeutronElasticXS::ElementCrossSecti
134 {
135 G4int Z = (ZZ >= MAXZEL) ? MAXZEL - 1 : ZZ;
136 auto pv = GetPhysicsVector(Z);
137
138 G4double xs = (ekin <= pv->GetMaxEnergy()) ?
139 : coeff[Z]*ggXsection->GetElasticElementCr
140
141
142 #ifdef G4VERBOSE
143 if(verboseLevel > 1) {
144 G4cout << "Z= " << Z << " Ekin(MeV)= " <<
145 << ", nElmXSel(b)= " << xs/CLHEP::barn
146 << G4endl;
147 }
148 #endif
149 return xs;
150 }
151
152 G4double
153 G4NeutronElasticXS::ComputeIsoCrossSection(G4d
154 const G4ParticleDefinition*
155 G4int Z, G4int A,
156 const G4Isotope*, const G4E
157 const G4Material*)
158 {
159 return ElementCrossSection(ekin, loge, Z)*A/
160 }
161
162 G4double
163 G4NeutronElasticXS::GetIsoCrossSection(const G
164 G4int Z, G4int A,
165 const G4Isotope*, const G4Eleme
166 const G4Material*)
167 {
168 return ElementCrossSection(aParticle->GetKin
169 aParticle->GetLogKineticEnergy(), Z
170
171 }
172
173 const G4Isotope* G4NeutronElasticXS::SelectIso
174 const G4Element* anElement, G4double, G4
175 {
176 G4int nIso = (G4int)anElement->GetNumberOfIs
177 const G4Isotope* iso = anElement->GetIsotope
178
179 //G4cout << "SelectIsotope NIso= " << nIso <
180 if(1 == nIso) { return iso; }
181
182 const G4double* abundVector = anElement->Get
183 G4double q = G4UniformRand();
184 G4double sum = 0.0;
185
186 // isotope wise cross section not used
187 for (G4int j=0; j<nIso; ++j) {
188 sum += abundVector[j];
189 if(q <= sum) {
190 iso = anElement->GetIsotope(j);
191 break;
192 }
193 }
194 return iso;
195 }
196
197 void
198 G4NeutronElasticXS::BuildPhysicsTable(const G4
199 {
200 if(verboseLevel > 0){
201 G4cout << "G4NeutronElasticXS::BuildPhysic
202 << p.GetParticleName() << G4endl;
203 }
204 if(p.GetParticleName() != "neutron") {
205 G4ExceptionDescription ed;
206 ed << p.GetParticleName() << " is a wrong
207 << " only neutron is allowed";
208 G4Exception("G4NeutronElasticXS::BuildPhys
209 FatalException, ed, "");
210 return;
211 }
212 if (fLock || isFirst) {
213 G4AutoLock l(&nElasticXSMutex);
214 if (fLock) {
215 isFirst = true;
216 fLock = false;
217 FindDirectoryPath();
218 }
219
220 // Access to elements
221 const G4ElementTable* table = G4Element::G
222 for ( auto & elm : *table ) {
223 G4int Z = std::max( 1, std::min( elm->Ge
224 if ( nullptr == data[Z] ) { Initialise(Z
225 }
226 l.unlock();
227 }
228 }
229
230 const G4String& G4NeutronElasticXS::FindDirect
231 {
232 // build the complete string identifying the
233 if (gDataDirectory.empty()) {
234 std::ostringstream ost;
235 ost << G4HadronicParameters::Instance()->G
236 gDataDirectory = ost.str();
237 }
238 return gDataDirectory;
239 }
240
241 void G4NeutronElasticXS::InitialiseOnFly(G4int
242 {
243 G4AutoLock l(&nElasticXSMutex);
244 Initialise(Z);
245 l.unlock();
246 }
247
248 void G4NeutronElasticXS::Initialise(G4int Z)
249 {
250 if(data[Z] != nullptr) { return; }
251
252 // upload data from file
253 data[Z] = new G4PhysicsLogVector();
254
255 std::ostringstream ost;
256 ost << FindDirectoryPath() << Z ;
257 std::ifstream filein(ost.str().c_str());
258 if (!filein.is_open()) {
259 G4ExceptionDescription ed;
260 ed << "Data file <" << ost.str().c_str()
261 << "> is not opened!";
262 G4Exception("G4NeutronElasticXS::Initialis
263 FatalException, ed, "Check G4P
264 return;
265 }
266 if(verboseLevel > 1) {
267 G4cout << "file " << ost.str()
268 << " is opened by G4NeutronElasticXS" <<
269 }
270
271 // retrieve data from DB
272 if(!data[Z]->Retrieve(filein, true)) {
273 G4ExceptionDescription ed;
274 ed << "Data file <" << ost.str().c_str()
275 << "> is not retrieved!";
276 G4Exception("G4NeutronElasticXS::Initialis
277 FatalException, ed, "Check G4PARTICLEXSDAT
278 return;
279 }
280 // smooth transition
281 G4double sig1 = (*(data[Z]))[data[Z]->GetVe
282 G4double ehigh = data[Z]->GetMaxEnergy();
283 G4double sig2 = ggXsection->GetElasticEleme
284 ehigh, Z, aeff[
285 coeff[Z] = (sig2 > 0.) ? sig1/sig2 : 1.0;
286 }
287