| 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 // particle_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans
29 //
30 // 070523 add neglecting doppler broadening on
31 // 070613 fix memory leaking by T. Koi
32 // 071002 enable cross section dump by T. Koi
33 // 080428 change checking point of "neglecting
34 // from GetCrossSection to BuildPhysics
35 // 081024 G4NucleiPropertiesTable:: to G4Nucle
36 //
37 // P. Arce, June-2014 Conversion neutron_hp to
38 //
39 #include "G4ParticleHPInelasticData.hh"
40
41 #include "G4ElementTable.hh"
42 #include "G4HadronicParameters.hh"
43 #include "G4Neutron.hh"
44 #include "G4NucleiProperties.hh"
45 #include "G4ParticleHPData.hh"
46 #include "G4ParticleHPManager.hh"
47 #include "G4Pow.hh"
48
49 G4ParticleHPInelasticData::G4ParticleHPInelast
50 : G4VCrossSectionDataSet("")
51 {
52 const char* dataDirVariable;
53 G4String particleName;
54 if (projectile == G4Neutron::Neutron()) {
55 dataDirVariable = "G4NEUTRONHPDATA";
56 }
57 else if (projectile == G4Proton::Proton()) {
58 dataDirVariable = "G4PROTONHPDATA";
59 particleName = "Proton";
60 }
61 else if (projectile == G4Deuteron::Deuteron(
62 dataDirVariable = "G4DEUTERONHPDATA";
63 particleName = "Deuteron";
64 }
65 else if (projectile == G4Triton::Triton()) {
66 dataDirVariable = "G4TRITONHPDATA";
67 particleName = "Triton";
68 }
69 else if (projectile == G4He3::He3()) {
70 dataDirVariable = "G4HE3HPDATA";
71 particleName = "He3";
72 }
73 else if (projectile == G4Alpha::Alpha()) {
74 dataDirVariable = "G4ALPHAHPDATA";
75 particleName = "Alpha";
76 }
77 else {
78 G4String message(
79 "G4ParticleHPInelasticData may only be c
80 "alpha, while it is called for "
81 + projectile->GetParticleName());
82 throw G4HadronicException(__FILE__, __LINE
83 }
84 G4String dataName = projectile->GetParticleN
85 dataName.at(0) = (char)std::toupper(dataName
86 SetName(dataName);
87
88 if ((G4FindDataDir(dataDirVariable) == nullp
89 {
90 G4String message("Please setenv G4PARTICLE
91 + G4String(dataDirVariabl
92 + projectile->GetParticle
93 throw G4HadronicException(__FILE__, __LINE
94 }
95
96 G4String dirName;
97 if (G4FindDataDir(dataDirVariable) != nullpt
98 dirName = G4FindDataDir(dataDirVariable);
99 }
100 else {
101 G4String baseName = G4FindDataDir("G4PARTI
102 dirName = baseName + "/" + particleName;
103 }
104 #ifdef G4VERBOSE
105 if (G4HadronicParameters::Instance()->GetVer
106 G4cout << "@@@ G4ParticleHPInelasticData i
107 << projectile->GetParticleName() <<
108 << " pointing to " << dirName << G4
109 }
110 #endif
111
112 SetMinKinEnergy(0 * CLHEP::MeV);
113 SetMaxKinEnergy(20 * CLHEP::MeV);
114
115 theCrossSections = nullptr;
116 theProjectile = projectile;
117
118 theHPData = nullptr;
119 instanceOfWorker = false;
120 if (G4Threading::IsMasterThread()) {
121 theHPData = new G4ParticleHPData(theProjec
122 }
123 else {
124 instanceOfWorker = true;
125 }
126 element_cache = nullptr;
127 material_cache = nullptr;
128 ke_cache = 0.0;
129 xs_cache = 0.0;
130 }
131
132 G4ParticleHPInelasticData::~G4ParticleHPInelas
133 {
134 if (theCrossSections != nullptr && !instance
135 theCrossSections->clearAndDestroy();
136 delete theCrossSections;
137 theCrossSections = nullptr;
138 }
139 if (theHPData != nullptr && !instanceOfWorke
140 delete theHPData;
141 theHPData = nullptr;
142 }
143 }
144
145 G4bool G4ParticleHPInelasticData::IsIsoApplica
146
147
148 {
149 G4double eKin = dp->GetKineticEnergy();
150 return eKin <= GetMaxKinEnergy() && eKin >=
151 && dp->GetDefinition() == theProjecti
152 }
153
154 G4double G4ParticleHPInelasticData::GetIsoCros
155
156
157
158 {
159 if (dp->GetKineticEnergy() == ke_cache && el
160 return xs_cache;
161
162 ke_cache = dp->GetKineticEnergy();
163 element_cache = element;
164 material_cache = material;
165 G4double xs = GetCrossSection(dp, element, m
166 xs_cache = xs;
167 return xs;
168 }
169
170 void G4ParticleHPInelasticData::BuildPhysicsTa
171 {
172 if (G4Threading::IsWorkerThread()) {
173 theCrossSections = G4ParticleHPManager::Ge
174 return;
175 }
176 if (theHPData == nullptr)
177 theHPData = G4ParticleHPData::Instance(con
178
179 std::size_t numberOfElements = G4Element::Ge
180 if (theCrossSections == nullptr)
181 theCrossSections = new G4PhysicsTable(numb
182 else
183 theCrossSections->clearAndDestroy();
184
185 // make a PhysicsVector for each element
186
187 auto theElementTable = G4Element::GetElement
188 for (std::size_t i = 0; i < numberOfElements
189 G4PhysicsVector* physVec = theHPData->Make
190 theCrossSections->push_back(physVec);
191 }
192 G4ParticleHPManager::GetInstance()->Register
193 }
194
195 void G4ParticleHPInelasticData::DumpPhysicsTab
196 {
197 #ifdef G4VERBOSE
198 if (G4HadronicParameters::Instance()->GetVer
199
200 // Dump element based cross section
201 // range 10e-5 eV to 20 MeV
202 // 10 point per decade
203 // in barn
204
205 G4cout << G4endl;
206 G4cout << G4endl;
207 G4cout << "Inelastic Cross Section of Neutro
208 G4cout << "(Pointwise cross-section at 0 Kel
209 G4cout << G4endl;
210 G4cout << "Name of Element" << G4endl;
211 G4cout << "Energy[eV] XS[barn]" << G4endl;
212 G4cout << G4endl;
213
214 std::size_t numberOfElements = G4Element::Ge
215 auto theElementTable = G4Element::GetElement
216
217 for (std::size_t i = 0; i < numberOfElements
218 G4cout << (*theElementTable)[i]->GetName()
219
220 G4int ie = 0;
221
222 for (ie = 0; ie < 130; ie++) {
223 G4double eKinetic = 1.0e-5 * G4Pow::GetI
224 G4bool outOfRange = false;
225
226 if (eKinetic < 20 * CLHEP::MeV) {
227 G4cout << eKinetic / CLHEP::eV << " "
228 << (*((*theCrossSections)(i))).
229 << G4endl;
230 }
231 }
232 G4cout << G4endl;
233 }
234 #endif
235 }
236
237 G4double G4ParticleHPInelasticData::GetCrossSe
238
239 {
240 G4double result = 0;
241 G4bool outOfRange;
242 std::size_t index = anE->GetIndex();
243
244 // prepare neutron
245 G4double eKinetic = projectile->GetKineticEn
246
247 if (G4ParticleHPManager::GetInstance()->GetN
248 // NEGLECT_DOPPLER
249 G4double factor = 1.0;
250 if (eKinetic < aT * CLHEP::k_Boltzmann) {
251 // below 0.1 eV neutrons
252 // Have to do some, but now just igonre.
253 // Will take care after performance chec
254 // factor = factor * targetV;
255 }
256 return ((*((*theCrossSections)(index))).Ge
257 }
258
259 G4ReactionProduct theNeutron(projectile->Get
260 theNeutron.SetMomentum(projectile->GetMoment
261 theNeutron.SetKineticEnergy(eKinetic);
262
263 // prepare thermal nucleus
264 G4Nucleus aNuc;
265 G4double eps = 0.0001;
266 G4double theA = anE->GetN();
267 G4double theZ = anE->GetZ();
268 G4double eleMass;
269 eleMass = G4NucleiProperties::GetNuclearMass
270
271
272 G4ReactionProduct boosted;
273 G4double aXsection;
274
275 // MC integration loop
276 G4int counter = 0;
277 G4int failCount = 0;
278 G4double buffer = 0;
279 G4int size = G4int(std::max(10., aT / 60 * C
280 G4ThreeVector neutronVelocity = 1. / theProj
281 G4double neutronVMag = neutronVelocity.mag()
282
283 #ifndef G4PHP_DOPPLER_LOOP_ONCE
284 while (counter == 0
285 || std::abs(buffer - result / std::ma
286 > 0.01 * buffer) // Loop checki
287 {
288 if (counter != 0) buffer = result / counte
289 while (counter < size) // Loop checking,
290 {
291 ++counter;
292 #endif
293 G4ReactionProduct aThermalNuc =
294 aNuc.GetThermalNucleus(eleMass / G4Neu
295 boosted.Lorentz(theNeutron, aThermalNuc)
296 G4double theEkin = boosted.GetKineticEne
297 aXsection = (*((*theCrossSections)(index
298 if (aXsection < 0) {
299 if (failCount < 1000) {
300 ++failCount;
301 #ifndef G4PHP_DOPPLER_LOOP_ONCE
302 --counter;
303 continue;
304 #endif
305 }
306 else {
307 aXsection = 0;
308 }
309 }
310 // velocity correction.
311 G4ThreeVector targetVelocity = 1. / aThe
312 aXsection *= (targetVelocity - neutronVe
313 result += aXsection;
314 }
315 #ifndef G4PHP_DOPPLER_LOOP_ONCE
316 size += size;
317 }
318 result /= counter;
319 #endif
320
321 return result;
322 }
323
324 G4int G4ParticleHPInelasticData::GetVerboseLev
325 {
326 return G4ParticleHPManager::GetInstance()->G
327 }
328
329 void G4ParticleHPInelasticData::SetVerboseLeve
330 {
331 G4ParticleHPManager::GetInstance()->SetVerbo
332 }
333
334 void G4ParticleHPInelasticData::CrossSectionDe
335 {
336 outFile << "Extension of High Precision cros
337 "deuteron, triton, He3 and alpha
338 }
339