| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4NeutronCaptureXS.cc,v 1.8 2011-01-09 02:37:48 dennis Exp $
>> 27 // GEANT4 tag $Name: geant4-09-04-patch-01 $
>> 28 //
26 // ------------------------------------------- 29 // -------------------------------------------------------------------
27 // 30 //
28 // GEANT4 Class file 31 // GEANT4 Class file
29 // 32 //
30 // 33 //
31 // File name: G4NeutronCaptureXS 34 // File name: G4NeutronCaptureXS
32 // 35 //
33 // Author Ivantchenko, Geant4, 3-Aug-09 36 // Author Ivantchenko, Geant4, 3-Aug-09
34 // 37 //
35 // Modifications: 38 // Modifications:
36 // 39 //
37 40
38 #include <fstream> << 41 #include "G4HadronicException.hh"
39 #include <sstream> <<
40 #include <thread> <<
41 <<
42 #include "G4SystemOfUnits.hh" <<
43 #include "G4NeutronCaptureXS.hh" 42 #include "G4NeutronCaptureXS.hh"
44 #include "G4Material.hh" <<
45 #include "G4Element.hh" 43 #include "G4Element.hh"
46 #include "G4PhysicsLogVector.hh" <<
47 #include "G4DynamicParticle.hh" <<
48 #include "G4ElementTable.hh" 44 #include "G4ElementTable.hh"
49 #include "G4IsotopeList.hh" << 45 #include "G4PhysicsLogVector.hh"
50 #include "G4HadronicParameters.hh" << 46 #include "G4PhysicsVector.hh"
51 #include "Randomize.hh" <<
52 #include "G4Log.hh" <<
53 #include "G4AutoLock.hh" <<
54 <<
55 G4ElementData* G4NeutronCaptureXS::data = null <<
56 G4String G4NeutronCaptureXS::gDataDirectory = <<
57 47
58 static std::once_flag applyOnce; << 48 #include <iostream>
>> 49 #include <fstream>
>> 50 #include <sstream>
59 51
60 namespace << 52 using namespace std;
61 { <<
62 G4Mutex neutronCaptureXSMutex = G4MUTEX_INIT <<
63 const G4int MAXZCAPTURE = 92; <<
64 } <<
65 53
66 G4NeutronCaptureXS::G4NeutronCaptureXS() 54 G4NeutronCaptureXS::G4NeutronCaptureXS()
67 : G4VCrossSectionDataSet(Default_Name()), << 55 : G4VCrossSectionDataSet("G4NeutronCaptureXS"),
68 emax(20*CLHEP::MeV), elimit(1.0e-5*CLHEP::e << 56 emax(20*MeV),maxZ(92)
69 { 57 {
70 verboseLevel = 0; << 58 // verboseLevel = 0;
71 if (verboseLevel > 0) { << 59 if(verboseLevel > 0){
72 G4cout << "G4NeutronCaptureXS::G4NeutronC 60 G4cout << "G4NeutronCaptureXS::G4NeutronCaptureXS: Initialise for Z < "
73 << MAXZCAPTURE << G4endl; << 61 << maxZ + 1 << G4endl;
74 } <<
75 logElimit = G4Log(elimit); <<
76 if (nullptr == data) { <<
77 data = new G4ElementData(MAXZCAPTURE+1); <<
78 data->SetName("nCapture"); <<
79 FindDirectoryPath(); <<
80 } 62 }
>> 63 data.resize(maxZ+1, 0);
>> 64 isInitialized = false;
81 } 65 }
82 66
83 void G4NeutronCaptureXS::CrossSectionDescripti << 67 G4NeutronCaptureXS::~G4NeutronCaptureXS()
84 { 68 {
85 outFile << "G4NeutronCaptureXS calculates th << 69 for(G4int i=0; i<=maxZ; ++i) {
86 << "on nuclei using data from the hi << 70 delete data[i];
87 << "These data are simplified and sm << 71 }
88 << "in order to reduce CPU time. G4N <<
89 << "above 20 MeV for all targets. Fo <<
90 << "Uranium is used.\n"; <<
91 } 72 }
92 << 73
93 G4bool 74 G4bool
94 G4NeutronCaptureXS::IsElementApplicable(const << 75 G4NeutronCaptureXS::IsApplicable(const G4DynamicParticle*,
95 G4int, const G4Material*) << 76 const G4Element*)
96 { 77 {
97 return true; 78 return true;
98 } 79 }
99 80
100 G4bool 81 G4bool
101 G4NeutronCaptureXS::IsIsoApplicable(const G4Dy 82 G4NeutronCaptureXS::IsIsoApplicable(const G4DynamicParticle*,
102 G4int, G4int, << 83 G4int /*ZZ*/, G4int /*AA*/)
103 const G4Element*, const G4Material <<
104 { 84 {
105 return true; << 85 return false;
106 } 86 }
107 87
108 G4double 88 G4double
109 G4NeutronCaptureXS::GetElementCrossSection(con << 89 G4NeutronCaptureXS::GetCrossSection(const G4DynamicParticle* aParticle,
110 G4int Z, const G4Material*) << 90 const G4Element* elm,
>> 91 G4double)
111 { 92 {
112 G4double xs = 0.0; 93 G4double xs = 0.0;
113 G4double ekin = aParticle->GetKineticEnergy( 94 G4double ekin = aParticle->GetKineticEnergy();
114 if (ekin < emax) { << 95 if(ekin > emax) { return xs; }
115 xs = ElementCrossSection(ekin, aParticle-> << 96 const G4double elimit = 1.0e-10*eV;
116 } << 97 if(ekin < elimit) { ekin = elimit; }
117 return xs; << 98
118 } << 99 G4int Z = G4int(elm->GetZ());
>> 100 G4PhysicsVector* pv = data[Z];
>> 101
>> 102 // element was not initialised
>> 103 if(!pv) {
>> 104 Initialise(Z);
>> 105 pv = data[Z];
>> 106 if(!pv) { return xs; }
>> 107 }
>> 108
>> 109 G4int n = pv->GetVectorLength() - 1;
>> 110 G4double e1 = pv->Energy(0);
>> 111 G4double e2 = pv->Energy(n);
>> 112 if(ekin < e1) { xs = (*pv)[0]*std::sqrt(e1/ekin); }
>> 113 else if(ekin <= e2) { xs = pv->Value(ekin); }
119 114
120 G4double << 115 if(verboseLevel > 0){
121 G4NeutronCaptureXS::ComputeCrossSectionPerElem << 116 G4cout << "ekin= " << ekin << ", xs= " << xs << G4endl;
122 const G4ParticleDefi <<
123 const G4Element* elm <<
124 const G4Material*) <<
125 { <<
126 G4double xs = 0.0; <<
127 if (ekin < emax) { <<
128 xs = ElementCrossSection(ekin, loge, elm-> <<
129 } 117 }
130 return xs; 118 return xs;
131 } 119 }
132 120
133 G4double <<
134 G4NeutronCaptureXS::ElementCrossSection(G4doub <<
135 { <<
136 G4int Z = std::min(ZZ, MAXZCAPTURE); <<
137 G4double ekin = eKin; <<
138 G4double logEkin = logE; <<
139 if (ekin < elimit) { <<
140 ekin = elimit; <<
141 logEkin = logElimit; <<
142 } <<
143 <<
144 auto pv = GetPhysicsVector(Z); <<
145 const G4double e0 = pv->Energy(0); <<
146 G4double xs = (ekin >= e0) ? pv->LogVectorVa <<
147 : (*pv)[0]*std::sqrt(e0/ekin); <<
148 <<
149 #ifdef G4VERBOSE <<
150 if (verboseLevel > 1){ <<
151 G4cout << "Ekin= " << ekin/CLHEP::MeV <<
152 << " ElmXScap(b)= " << xs/CLHEP::b <<
153 } <<
154 #endif <<
155 return xs; <<
156 } <<
157 <<
158 G4double <<
159 G4NeutronCaptureXS::ComputeIsoCrossSection(G4d <<
160 const G4ParticleDefinition* <<
161 G4int Z, G4int A, <<
162 const G4Isotope*, const G4E <<
163 const G4Material*) <<
164 { <<
165 return IsoCrossSection(ekin, loge, Z, A); <<
166 } <<
167 <<
168 G4double <<
169 G4NeutronCaptureXS::GetIsoCrossSection(const G <<
170 G4int Z, G4int A, <<
171 const G4Isotope*, const G4Eleme <<
172 const G4Material*) <<
173 { <<
174 return IsoCrossSection(aParticle->GetKinetic <<
175 aParticle->GetLogKine <<
176 Z, A); <<
177 } <<
178 <<
179 G4double G4NeutronCaptureXS::IsoCrossSection(G <<
180 G <<
181 { <<
182 G4double xs = 0.0; <<
183 if (eKin > emax) { return xs; } <<
184 <<
185 G4int Z = std::min(ZZ, MAXZCAPTURE); <<
186 G4double ekin = eKin; <<
187 G4double logEkin = logE; <<
188 if (ekin < elimit) { <<
189 ekin = elimit; <<
190 logEkin = logElimit; <<
191 } <<
192 <<
193 auto pv = GetPhysicsVector(Z); <<
194 if (pv == nullptr) { return xs; } <<
195 <<
196 // use isotope x-section if possible <<
197 if (data->GetNumberOfComponents(Z) > 0) { <<
198 G4PhysicsVector* pviso = data->GetComponen <<
199 if(pviso != nullptr) { <<
200 const G4double e0 = pviso->Energy(0); <<
201 xs = (ekin >= e0) ? pviso->LogVectorValu <<
202 : (*pviso)[0]*std::sqrt(e0/ekin); <<
203 #ifdef G4VERBOSE <<
204 if(verboseLevel > 0) { <<
205 G4cout << "G4NeutronCaptureXS::IsoXS: Ekin(M <<
206 << " xs(b)= " << xs/barn <<
207 << " Z= " << Z << " A= " << A << G4 <<
208 } <<
209 #endif <<
210 return xs; <<
211 } <<
212 } <<
213 // isotope data are not available or applica <<
214 const G4double e0 = pv->Energy(0); <<
215 xs = (ekin >= e0) ? pv->LogVectorValue(ekin, <<
216 : (*pv)[0]*std::sqrt(e0/ekin); <<
217 #ifdef G4VERBOSE <<
218 if (verboseLevel > 0) { <<
219 G4cout << "G4NeutronCaptureXS::IsoXS: Ekin <<
220 << " xs(b)= " << xs/barn <<
221 << " Z= " << Z << " A= " << A << " no i <<
222 } <<
223 #endif <<
224 return xs; <<
225 } <<
226 <<
227 const G4Isotope* <<
228 G4NeutronCaptureXS::SelectIsotope(const G4Elem <<
229 G4double kinEnergy, G4double logE) <<
230 { <<
231 G4int nIso = (G4int)anElement->GetNumberOfIs <<
232 const G4Isotope* iso = anElement->GetIsotope <<
233 <<
234 //G4cout << "SelectIsotope NIso= " << nIso < <<
235 if(1 == nIso) { return iso; } <<
236 <<
237 // more than 1 isotope <<
238 G4int Z = anElement->GetZasInt(); <<
239 if (nullptr == data->GetElementData(Z)) { In <<
240 <<
241 const G4double* abundVector = anElement->Get <<
242 G4double q = G4UniformRand(); <<
243 G4double sum = 0.0; <<
244 <<
245 // is there isotope wise cross section? <<
246 G4int j; <<
247 if (Z > MAXZCAPTURE || 0 == data->GetNumberO <<
248 for (j = 0; j<nIso; ++j) { <<
249 sum += abundVector[j]; <<
250 if(q <= sum) { <<
251 iso = anElement->GetIsotope(j); <<
252 break; <<
253 } <<
254 } <<
255 return iso; <<
256 } <<
257 G4int nn = (G4int)temp.size(); <<
258 if (nn < nIso) { temp.resize(nIso, 0.); } <<
259 <<
260 for (j=0; j<nIso; ++j) { <<
261 sum += abundVector[j]*IsoCrossSection(kinE <<
262 anElement->GetIsotope(j)->GetN()); <<
263 temp[j] = sum; <<
264 } <<
265 sum *= q; <<
266 for (j = 0; j<nIso; ++j) { <<
267 if (temp[j] >= sum) { <<
268 iso = anElement->GetIsotope(j); <<
269 break; <<
270 } <<
271 } <<
272 return iso; <<
273 } <<
274 <<
275 void 121 void
276 G4NeutronCaptureXS::BuildPhysicsTable(const G4 122 G4NeutronCaptureXS::BuildPhysicsTable(const G4ParticleDefinition& p)
277 { 123 {
278 if (verboseLevel > 0){ << 124 if(verboseLevel > 0){
279 G4cout << "G4NeutronCaptureXS::BuildPhysic 125 G4cout << "G4NeutronCaptureXS::BuildPhysicsTable for "
280 << p.GetParticleName() << G4endl; 126 << p.GetParticleName() << G4endl;
281 } 127 }
282 if (p.GetParticleName() != "neutron") { << 128 if(isInitialized || p.GetParticleName() != "neutron") { return; }
283 G4ExceptionDescription ed; << 129 isInitialized = true;
284 ed << p.GetParticleName() << " is a wrong <<
285 << " only neutron is allowed"; <<
286 G4Exception("G4NeutronCaptureXS::BuildPhys <<
287 FatalException, ed, ""); <<
288 return; <<
289 } <<
290 <<
291 // it is possible re-initialisation for the <<
292 const G4ElementTable* table = G4Element::Get <<
293 130
294 // initialise static tables only once << 131 // check environment variable
295 std::call_once(applyOnce, [this]() { isIniti << 132 // Build the complete string identifying the file with the data set
296 << 133 char* path = getenv("G4NEUTRONXSDATA");
297 if (isInitializer) { << 134 if (!path){
298 G4AutoLock l(&neutronCaptureXSMutex); << 135 G4cout << "G4NEUTRONXSDATA environment variable not set" << G4endl;
299 // Access to elements << 136 }
300 for ( auto const & elm : *table ) { << 137
301 G4int Z = std::max( 1, std::min( elm->Ge << 138 // Access to elements
302 if ( nullptr == data->GetElementData(Z) << 139 const G4ElementTable* theElmTable = G4Element::GetElementTable();
>> 140 size_t numOfElm = G4Element::GetNumberOfElements();
>> 141 if(numOfElm > 0) {
>> 142 for(size_t i=0; i<numOfElm; ++i) {
>> 143 G4int Z = G4int(((*theElmTable)[i])->GetZ());
>> 144 if(Z < 1) { Z = 1; }
>> 145 else if(Z > maxZ) { Z = maxZ; }
>> 146 //G4cout << "Z= " << Z << G4endl;
>> 147 // Initialisation
>> 148 if(!data[Z]) { Initialise(Z, path); }
303 } 149 }
304 l.unlock(); <<
305 } <<
306 <<
307 // prepare isotope selection <<
308 std::size_t nIso = temp.size(); <<
309 for ( auto const & elm : *table ) { <<
310 std::size_t n = elm->GetNumberOfIsotopes() <<
311 if (n > nIso) { nIso = n; } <<
312 } 150 }
313 temp.resize(nIso, 0.0); <<
314 } 151 }
315 152
316 const G4String& G4NeutronCaptureXS::FindDirect << 153 void
317 { << 154 G4NeutronCaptureXS::DumpPhysicsTable(const G4ParticleDefinition&)
318 // build the complete string identifying the << 155 {}
319 if(gDataDirectory.empty()) { <<
320 std::ostringstream ost; <<
321 ost << G4HadronicParameters::Instance()->G <<
322 gDataDirectory = ost.str(); <<
323 } <<
324 return gDataDirectory; <<
325 } <<
326 156
327 void G4NeutronCaptureXS::InitialiseOnFly(G4int << 157 void
>> 158 G4NeutronCaptureXS::Initialise(G4int Z, const char* p)
328 { 159 {
329 G4AutoLock l(&neutronCaptureXSMutex); << 160 if(data[Z]) { return; }
330 Initialise(Z); << 161 const char* path = p;
331 l.unlock(); << 162 if(!p) {
332 } << 163 // check environment variable
>> 164 // Build the complete string identifying the file with the data set
>> 165 path = getenv("G4NEUTRONXSDATA");
>> 166 if (!path) {
>> 167 if(verboseLevel > 1) {
>> 168 G4cout << "G4NEUTRONXSDATA environment variable not set" << G4endl;
>> 169 }
>> 170 return;
>> 171 }
>> 172 }
333 173
334 void G4NeutronCaptureXS::Initialise(G4int Z) << 174 // upload data from file
335 { << 175 data[Z] = new G4PhysicsLogVector();
336 if (nullptr != data->GetElementData(Z)) { re <<
337 176
338 // upload element data <<
339 std::ostringstream ost; 177 std::ostringstream ost;
340 ost << FindDirectoryPath() << Z ; << 178 ost << path << "/cap" << Z ;
341 G4PhysicsVector* v = RetrieveVector(ost, tru <<
342 data->InitialiseForElement(Z, v); <<
343 <<
344 // upload isotope data <<
345 G4bool noComp = true; <<
346 if (amin[Z] < amax[Z]) { <<
347 for(G4int A=amin[Z]; A<=amax[Z]; ++A) { <<
348 std::ostringstream ost1; <<
349 ost1 << gDataDirectory << Z << "_" << A; <<
350 G4PhysicsVector* v1 = RetrieveVector(ost <<
351 if (nullptr != v1) { <<
352 if (noComp) { <<
353 G4int nmax = amax[Z] - A + 1; <<
354 data->InitialiseForComponent(Z, nmax); <<
355 noComp = false; <<
356 } <<
357 data->AddComponent(Z, A, v1); <<
358 } <<
359 } <<
360 } <<
361 // no components case <<
362 if (noComp) { data->InitialiseForComponent(Z <<
363 } <<
364 <<
365 G4PhysicsVector* <<
366 G4NeutronCaptureXS::RetrieveVector(std::ostrin <<
367 { <<
368 G4PhysicsLogVector* v = nullptr; <<
369 std::ifstream filein(ost.str().c_str()); 179 std::ifstream filein(ost.str().c_str());
370 if (!filein.is_open()) { << 180 if (!(filein)) {
371 if (warn) { << 181 G4cout << ost.str() << " is not opened by G4NeutronCaptureXS" << G4endl;
372 G4ExceptionDescription ed; << 182 throw G4HadronicException(__FILE__, __LINE__,
373 ed << "Data file <" << ost.str().c_str() << 183 "G4NeutronCaptureXS: no data sets registered");
374 << "> is not opened!"; << 184 return;
375 G4Exception("G4NeutronCaptureXS::Retriev << 185 }else{
376 FatalException, ed, "Check G4PARTICLEXSD << 186 if(verboseLevel > 1) {
377 } <<
378 } else { <<
379 if (verboseLevel > 1) { <<
380 G4cout << "File " << ost.str() 187 G4cout << "File " << ost.str()
381 << " is opened by G4NeutronCaptureXS" < 188 << " is opened by G4NeutronCaptureXS" << G4endl;
382 } 189 }
383 // retrieve data from DB 190 // retrieve data from DB
384 v = new G4PhysicsLogVector(); << 191 data[Z]->Retrieve(filein, true);
385 if (!v->Retrieve(filein, true)) { << 192 }
386 G4ExceptionDescription ed; <<
387 ed << "Data file <" << ost.str().c_str() <<
388 << "> is not retrieved!"; <<
389 G4Exception("G4NeutronCaptureXS::Retriev <<
390 FatalException, ed, "Check G4PARTICLEXSD <<
391 } <<
392 } <<
393 return v; <<
394 } 193 }
>> 194
395 195