| 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 // ------------------------------------------- 26 // -------------------------------------------------------------------
27 // 27 //
28 // GEANT4 Class file 28 // GEANT4 Class file
29 // 29 //
30 // 30 //
31 // File name: G4NeutronElasticXS 31 // File name: G4NeutronElasticXS
32 // 32 //
33 // Author Ivantchenko, Geant4, 3-Aug-09 33 // Author Ivantchenko, Geant4, 3-Aug-09
34 // 34 //
35 // Modifications: 35 // Modifications:
36 // 36 //
37 37
38 #include "G4NeutronElasticXS.hh" 38 #include "G4NeutronElasticXS.hh"
39 #include "G4Neutron.hh" 39 #include "G4Neutron.hh"
40 #include "G4DynamicParticle.hh" 40 #include "G4DynamicParticle.hh"
41 #include "G4ElementTable.hh" << 41 #include "G4ProductionCutsTable.hh"
42 #include "G4Material.hh" 42 #include "G4Material.hh"
43 #include "G4Element.hh" 43 #include "G4Element.hh"
44 #include "G4PhysicsLogVector.hh" 44 #include "G4PhysicsLogVector.hh"
45 #include "G4CrossSectionDataSetRegistry.hh" << 45 #include "G4PhysicsVector.hh"
46 #include "G4ComponentGGHadronNucleusXsc.hh" 46 #include "G4ComponentGGHadronNucleusXsc.hh"
47 #include "G4HadronicParameters.hh" << 47 #include "G4NistManager.hh"
48 #include "Randomize.hh" <<
49 #include "G4SystemOfUnits.hh" 48 #include "G4SystemOfUnits.hh"
50 #include "G4IsotopeList.hh" <<
51 #include "G4AutoLock.hh" <<
52 49
53 #include <fstream> 50 #include <fstream>
54 #include <sstream> 51 #include <sstream>
55 52
>> 53 // factory
>> 54 #include "G4CrossSectionFactory.hh"
>> 55 //
>> 56 G4_DECLARE_XS_FACTORY(G4NeutronElasticXS);
>> 57
>> 58 using namespace std;
>> 59
56 G4PhysicsVector* G4NeutronElasticXS::data[] = 60 G4PhysicsVector* G4NeutronElasticXS::data[] = {nullptr};
57 G4double G4NeutronElasticXS::coeff[] = {0.0}; 61 G4double G4NeutronElasticXS::coeff[] = {0.0};
>> 62 G4double G4NeutronElasticXS::aeff[] = {1.0};
58 G4String G4NeutronElasticXS::gDataDirectory = 63 G4String G4NeutronElasticXS::gDataDirectory = "";
59 G4bool G4NeutronElasticXS::fLock = true; <<
60 64
61 namespace << 65 #ifdef G4MULTITHREADED
62 { << 66 G4Mutex G4NeutronElasticXS::neutronElasticXSMutex = G4MUTEX_INITIALIZER;
63 G4Mutex nElasticXSMutex = G4MUTEX_INITIALIZE << 67 #endif
64 } <<
65 68
66 G4NeutronElasticXS::G4NeutronElasticXS() 69 G4NeutronElasticXS::G4NeutronElasticXS()
67 : G4VCrossSectionDataSet(Default_Name()), 70 : G4VCrossSectionDataSet(Default_Name()),
68 neutron(G4Neutron::Neutron()) << 71 ggXsection(nullptr),
>> 72 neutron(G4Neutron::Neutron()),
>> 73 isMaster(false)
69 { 74 {
70 // verboseLevel = 0; 75 // verboseLevel = 0;
71 if (verboseLevel > 0){ << 76 if(verboseLevel > 0){
72 G4cout << "G4NeutronElasticXS::G4NeutronE 77 G4cout << "G4NeutronElasticXS::G4NeutronElasticXS Initialise for Z < "
73 << MAXZEL << G4endl; 78 << MAXZEL << G4endl;
74 } 79 }
75 ggXsection = << 80 nist = G4NistManager::Instance();
76 G4CrossSectionDataSetRegistry::Instance()- << 81 ggXsection = new G4ComponentGGHadronNucleusXsc();
77 if (ggXsection == nullptr) <<
78 ggXsection = new G4ComponentGGHadronNucleu <<
79 SetForAllAtomsAndEnergies(true); 82 SetForAllAtomsAndEnergies(true);
80 FindDirectoryPath(); << 83 temp.resize(13,0.0);
81 } 84 }
82 85
83 G4NeutronElasticXS::~G4NeutronElasticXS() 86 G4NeutronElasticXS::~G4NeutronElasticXS()
84 { 87 {
85 if (isFirst) { << 88 if(isMaster) {
86 for(G4int i=0; i<MAXZEL; ++i) { 89 for(G4int i=0; i<MAXZEL; ++i) {
87 delete data[i]; 90 delete data[i];
88 data[i] = nullptr; 91 data[i] = nullptr;
89 } 92 }
90 } 93 }
91 } 94 }
92 95
93 void G4NeutronElasticXS::CrossSectionDescripti 96 void G4NeutronElasticXS::CrossSectionDescription(std::ostream& outFile) const
94 { 97 {
95 outFile << "G4NeutronElasticXS calculates th 98 outFile << "G4NeutronElasticXS calculates the neutron elastic scattering\n"
96 << "cross section on nuclei using da 99 << "cross section on nuclei using data from the high precision\n"
97 << "neutron database. These data ar 100 << "neutron database. These data are simplified and smoothed over\n"
98 << "the resonance region in order to 101 << "the resonance region in order to reduce CPU time.\n"
99 << "For high energies Glauber-Gribiv 102 << "For high energies Glauber-Gribiv cross section is used.\n";
100 } 103 }
101 104
102 G4bool 105 G4bool
103 G4NeutronElasticXS::IsElementApplicable(const 106 G4NeutronElasticXS::IsElementApplicable(const G4DynamicParticle*,
104 G4int, const G4Material*) 107 G4int, const G4Material*)
105 { 108 {
106 return true; 109 return true;
107 } 110 }
108 111
109 G4bool G4NeutronElasticXS::IsIsoApplicable(con 112 G4bool G4NeutronElasticXS::IsIsoApplicable(const G4DynamicParticle*,
110 G4i 113 G4int, G4int,
111 con 114 const G4Element*, const G4Material*)
112 { 115 {
113 return false; << 116 return true;
114 } 117 }
115 118
116 G4double 119 G4double
117 G4NeutronElasticXS::GetElementCrossSection(con 120 G4NeutronElasticXS::GetElementCrossSection(const G4DynamicParticle* aParticle,
118 G4int Z, const G4Material*) << 121 G4int ZZ, const G4Material*)
119 { 122 {
120 return ElementCrossSection(aParticle->GetKin << 123 G4double xs = 0.0;
121 aParticle->GetLogKineticEnergy(), Z << 124 G4double ekin = aParticle->GetKineticEnergy();
122 } <<
123 125
124 G4double << 126 G4int Z = (ZZ >= MAXZEL) ? MAXZEL - 1 : ZZ;
125 G4NeutronElasticXS::ComputeCrossSectionPerElem <<
126 const G4ParticleDefinition*, <<
127 const G4Element* elm, <<
128 const G4Material*) <<
129 { <<
130 return ElementCrossSection(ekin, loge, elm-> <<
131 } <<
132 127
133 G4double G4NeutronElasticXS::ElementCrossSecti <<
134 { <<
135 G4int Z = (ZZ >= MAXZEL) ? MAXZEL - 1 : ZZ; <<
136 auto pv = GetPhysicsVector(Z); 128 auto pv = GetPhysicsVector(Z);
>> 129 if(!pv) { return xs; }
>> 130 // G4cout << "G4NeutronElasticXS::GetCrossSection e= " << ekin
>> 131 // << " Z= " << Z << G4endl;
>> 132
>> 133 if(ekin <= pv->Energy(0)) {
>> 134 xs = (*pv)[0];
>> 135 } else if(ekin <= pv->GetMaxEnergy()) {
>> 136 xs = pv->LogVectorValue(ekin, aParticle->GetLogKineticEnergy());
>> 137 } else {
>> 138 xs = coeff[Z]*ggXsection->GetElasticElementCrossSection(neutron,
>> 139 ekin, Z, aeff[Z]);
>> 140 }
137 141
138 G4double xs = (ekin <= pv->GetMaxEnergy()) ? <<
139 : coeff[Z]*ggXsection->GetElasticElementCr <<
140 <<
141 <<
142 #ifdef G4VERBOSE <<
143 if(verboseLevel > 1) { 142 if(verboseLevel > 1) {
144 G4cout << "Z= " << Z << " Ekin(MeV)= " << 143 G4cout << "Z= " << Z << " Ekin(MeV)= " << ekin/CLHEP::MeV
145 << ", nElmXSel(b)= " << xs/CLHEP::barn 144 << ", nElmXSel(b)= " << xs/CLHEP::barn
146 << G4endl; 145 << G4endl;
147 } 146 }
148 #endif <<
149 return xs; 147 return xs;
150 } 148 }
151 149
152 G4double << 150 G4double G4NeutronElasticXS::GetIsoCrossSection(
153 G4NeutronElasticXS::ComputeIsoCrossSection(G4d << 151 const G4DynamicParticle* aParticle,
154 const G4ParticleDefinition* << 152 G4int Z, G4int A,
155 G4int Z, G4int A, << 153 const G4Isotope*, const G4Element*,
156 const G4Isotope*, const G4E << 154 const G4Material*)
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 { 155 {
168 return ElementCrossSection(aParticle->GetKin << 156 return IsoCrossSection(aParticle->GetKineticEnergy(),
169 aParticle->GetLogKineticEnergy(), Z << 157 aParticle->GetLogKineticEnergy(), Z, A);
>> 158 }
>> 159
>> 160 G4double
>> 161 G4NeutronElasticXS::IsoCrossSection(G4double ekin, G4double logekin,
>> 162 G4int ZZ, G4int A)
>> 163 {
>> 164 G4double xs = 0.0;
>> 165 G4int Z = (ZZ >= MAXZEL) ? MAXZEL - 1 : ZZ;
>> 166
>> 167 // tritium and He3
>> 168 if(3 == A) {
>> 169 return ggXsection->GetElasticElementCrossSection(neutron, ekin, Z, A);
>> 170 }
>> 171 /*
>> 172 G4cout << "IsoCrossSection Z= " << Z << " A= " << A
>> 173 << " Amin= " << amin[Z] << " Amax= " << amax[Z]
>> 174 << " E(MeV)= " << ekin << G4endl;
>> 175 */
>> 176 auto pv = GetPhysicsVector(Z);
>> 177 if(!pv) { return xs; }
170 178
>> 179 if(ekin <= pv->Energy(0)) {
>> 180 xs = (*pv)[0];
>> 181 } else if(ekin <= pv->GetMaxEnergy()) {
>> 182 xs = pv->LogVectorValue(ekin, logekin);
>> 183 } else {
>> 184 xs = coeff[Z]*ggXsection->GetElasticElementCrossSection(neutron,
>> 185 ekin, Z, aeff[Z]);
>> 186 }
>> 187 xs *= A/aeff[Z];
>> 188 if(verboseLevel > 1) {
>> 189 G4cout << "G4NeutronElasticXS::IsoXS: Z= " << Z << " A= " << A
>> 190 << " Ekin(MeV)= " << ekin/CLHEP::MeV
>> 191 << ", ElmXS(b)= " << xs/CLHEP::barn << G4endl;
>> 192 }
>> 193 return xs;
171 } 194 }
172 195
173 const G4Isotope* G4NeutronElasticXS::SelectIso 196 const G4Isotope* G4NeutronElasticXS::SelectIsotope(
174 const G4Element* anElement, G4double, G4 << 197 const G4Element* anElement, G4double kinEnergy, G4double logE)
175 { 198 {
176 G4int nIso = (G4int)anElement->GetNumberOfIs << 199 size_t nIso = anElement->GetNumberOfIsotopes();
177 const G4Isotope* iso = anElement->GetIsotope 200 const G4Isotope* iso = anElement->GetIsotope(0);
178 201
179 //G4cout << "SelectIsotope NIso= " << nIso < 202 //G4cout << "SelectIsotope NIso= " << nIso << G4endl;
180 if(1 == nIso) { return iso; } 203 if(1 == nIso) { return iso; }
181 204
>> 205 // more than 1 isotope
>> 206 G4int Z = anElement->GetZasInt();
>> 207 //G4cout << "SelectIsotope Z= " << Z << G4endl;
>> 208
182 const G4double* abundVector = anElement->Get 209 const G4double* abundVector = anElement->GetRelativeAbundanceVector();
183 G4double q = G4UniformRand(); 210 G4double q = G4UniformRand();
184 G4double sum = 0.0; 211 G4double sum = 0.0;
>> 212 size_t j;
185 213
186 // isotope wise cross section not used 214 // isotope wise cross section not used
187 for (G4int j=0; j<nIso; ++j) { << 215 if(anElement->GetNaturalAbundanceFlag()) {
188 sum += abundVector[j]; << 216 for (j=0; j<nIso; ++j) {
189 if(q <= sum) { << 217 sum += abundVector[j];
>> 218 if(q <= sum) {
>> 219 iso = anElement->GetIsotope(j);
>> 220 break;
>> 221 }
>> 222 }
>> 223 return iso;
>> 224 }
>> 225
>> 226 // use isotope cross sections
>> 227 size_t nn = temp.size();
>> 228 if(nn < nIso) { temp.resize(nIso, 0.); }
>> 229
>> 230 for (j=0; j<nIso; ++j) {
>> 231 //G4cout << j << "-th isotope " << (*isoVector)[j]->GetN()
>> 232 // << " abund= " << abundVector[j] << G4endl;
>> 233 sum += abundVector[j]*IsoCrossSection(kinEnergy, logE, Z,
>> 234 anElement->GetIsotope(j)->GetN());
>> 235 temp[j] = sum;
>> 236 }
>> 237 sum *= q;
>> 238 for (j = 0; j<nIso; ++j) {
>> 239 if(temp[j] >= sum) {
190 iso = anElement->GetIsotope(j); 240 iso = anElement->GetIsotope(j);
191 break; 241 break;
192 } 242 }
193 } 243 }
194 return iso; 244 return iso;
195 } 245 }
196 246
197 void 247 void
198 G4NeutronElasticXS::BuildPhysicsTable(const G4 248 G4NeutronElasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
199 { 249 {
200 if(verboseLevel > 0){ 250 if(verboseLevel > 0){
201 G4cout << "G4NeutronElasticXS::BuildPhysic 251 G4cout << "G4NeutronElasticXS::BuildPhysicsTable for "
202 << p.GetParticleName() << G4endl; 252 << p.GetParticleName() << G4endl;
203 } 253 }
204 if(p.GetParticleName() != "neutron") { 254 if(p.GetParticleName() != "neutron") {
205 G4ExceptionDescription ed; 255 G4ExceptionDescription ed;
206 ed << p.GetParticleName() << " is a wrong 256 ed << p.GetParticleName() << " is a wrong particle type -"
207 << " only neutron is allowed"; 257 << " only neutron is allowed";
208 G4Exception("G4NeutronElasticXS::BuildPhys 258 G4Exception("G4NeutronElasticXS::BuildPhysicsTable(..)","had012",
209 FatalException, ed, ""); 259 FatalException, ed, "");
210 return; 260 return;
211 } 261 }
212 if (fLock || isFirst) { << 262 if(0. == coeff[0]) {
213 G4AutoLock l(&nElasticXSMutex); << 263 #ifdef G4MULTITHREADED
214 if (fLock) { << 264 G4MUTEXLOCK(&neutronElasticXSMutex);
215 isFirst = true; << 265 if(0. == coeff[0]) {
216 fLock = false; << 266 #endif
217 FindDirectoryPath(); << 267 coeff[0] = 1.0;
>> 268 isMaster = true;
>> 269 #ifdef G4MULTITHREADED
218 } 270 }
>> 271 G4MUTEXUNLOCK(&neutronElasticXSMutex);
>> 272 #endif
>> 273 }
>> 274
>> 275 // it is possible re-initialisation for the second run
>> 276 if(isMaster) {
219 277
220 // Access to elements 278 // Access to elements
221 const G4ElementTable* table = G4Element::G << 279 auto theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
222 for ( auto & elm : *table ) { << 280 size_t numOfCouples = theCoupleTable->GetTableSize();
223 G4int Z = std::max( 1, std::min( elm->Ge << 281 for(size_t j=0; j<numOfCouples; ++j) {
224 if ( nullptr == data[Z] ) { Initialise(Z << 282 auto mat = theCoupleTable->GetMaterialCutsCouple(j)->GetMaterial();
>> 283 auto elmVec = mat->GetElementVector();
>> 284 size_t numOfElem = mat->GetNumberOfElements();
>> 285 for (size_t ie = 0; ie < numOfElem; ++ie) {
>> 286 G4int Z = std::max(1,std::min(((*elmVec)[ie])->GetZasInt(), MAXZEL-1));
>> 287 if(!data[Z]) { Initialise(Z); }
>> 288 }
225 } 289 }
226 l.unlock(); <<
227 } 290 }
228 } 291 }
229 292
>> 293 G4PhysicsVector* G4NeutronElasticXS::GetPhysicsVector(G4int Z)
>> 294 {
>> 295 if(!data[Z]) { InitialiseOnFly(Z); }
>> 296 return data[Z];
>> 297 }
>> 298
230 const G4String& G4NeutronElasticXS::FindDirect 299 const G4String& G4NeutronElasticXS::FindDirectoryPath()
231 { 300 {
>> 301 // check environment variable
232 // build the complete string identifying the 302 // build the complete string identifying the file with the data set
233 if (gDataDirectory.empty()) { << 303 if(gDataDirectory.empty()) {
234 std::ostringstream ost; << 304 char* path = std::getenv("G4PARTICLEXSDATA");
235 ost << G4HadronicParameters::Instance()->G << 305 if (path) {
236 gDataDirectory = ost.str(); << 306 std::ostringstream ost;
>> 307 ost << path << "/neutron/el";
>> 308 gDataDirectory = ost.str();
>> 309 } else {
>> 310 G4Exception("G4NeutronElasticXS::Initialise(..)","had013",
>> 311 FatalException,
>> 312 "Environment variable G4PARTICLEXSDATA is not defined");
>> 313 }
237 } 314 }
238 return gDataDirectory; 315 return gDataDirectory;
239 } 316 }
240 317
241 void G4NeutronElasticXS::InitialiseOnFly(G4int 318 void G4NeutronElasticXS::InitialiseOnFly(G4int Z)
242 { 319 {
243 G4AutoLock l(&nElasticXSMutex); << 320 #ifdef G4MULTITHREADED
244 Initialise(Z); << 321 G4MUTEXLOCK(&neutronElasticXSMutex);
245 l.unlock(); << 322 if(!data[Z]) {
>> 323 #endif
>> 324 Initialise(Z);
>> 325 #ifdef G4MULTITHREADED
>> 326 }
>> 327 G4MUTEXUNLOCK(&neutronElasticXSMutex);
>> 328 #endif
246 } 329 }
247 330
248 void G4NeutronElasticXS::Initialise(G4int Z) 331 void G4NeutronElasticXS::Initialise(G4int Z)
249 { 332 {
250 if(data[Z] != nullptr) { return; } << 333 if(data[Z]) { return; }
251 334
252 // upload data from file 335 // upload data from file
253 data[Z] = new G4PhysicsLogVector(); 336 data[Z] = new G4PhysicsLogVector();
254 337
255 std::ostringstream ost; 338 std::ostringstream ost;
256 ost << FindDirectoryPath() << Z ; 339 ost << FindDirectoryPath() << Z ;
257 std::ifstream filein(ost.str().c_str()); 340 std::ifstream filein(ost.str().c_str());
258 if (!filein.is_open()) { << 341 if (!(filein)) {
259 G4ExceptionDescription ed; 342 G4ExceptionDescription ed;
260 ed << "Data file <" << ost.str().c_str() 343 ed << "Data file <" << ost.str().c_str()
261 << "> is not opened!"; 344 << "> is not opened!";
262 G4Exception("G4NeutronElasticXS::Initialis 345 G4Exception("G4NeutronElasticXS::Initialise(..)","had014",
263 FatalException, ed, "Check G4P 346 FatalException, ed, "Check G4PARTICLEXSDATA");
264 return; 347 return;
265 } 348 }
266 if(verboseLevel > 1) { 349 if(verboseLevel > 1) {
267 G4cout << "file " << ost.str() 350 G4cout << "file " << ost.str()
268 << " is opened by G4NeutronElasticXS" << 351 << " is opened by G4NeutronElasticXS" << G4endl;
269 } 352 }
270 353
271 // retrieve data from DB 354 // retrieve data from DB
272 if(!data[Z]->Retrieve(filein, true)) { 355 if(!data[Z]->Retrieve(filein, true)) {
273 G4ExceptionDescription ed; 356 G4ExceptionDescription ed;
274 ed << "Data file <" << ost.str().c_str() 357 ed << "Data file <" << ost.str().c_str()
275 << "> is not retrieved!"; 358 << "> is not retrieved!";
276 G4Exception("G4NeutronElasticXS::Initialis 359 G4Exception("G4NeutronElasticXS::Initialise(..)","had015",
277 FatalException, ed, "Check G4PARTICLEXSDAT 360 FatalException, ed, "Check G4PARTICLEXSDATA");
278 return; 361 return;
279 } 362 }
280 // smooth transition 363 // smooth transition
281 G4double sig1 = (*(data[Z]))[data[Z]->GetVe 364 G4double sig1 = (*(data[Z]))[data[Z]->GetVectorLength()-1];
282 G4double ehigh = data[Z]->GetMaxEnergy(); 365 G4double ehigh = data[Z]->GetMaxEnergy();
>> 366 aeff[Z] = nist->GetAtomicMassAmu(Z);
283 G4double sig2 = ggXsection->GetElasticEleme 367 G4double sig2 = ggXsection->GetElasticElementCrossSection(neutron,
284 ehigh, Z, aeff[ 368 ehigh, Z, aeff[Z]);
285 coeff[Z] = (sig2 > 0.) ? sig1/sig2 : 1.0; << 369 if(sig2 > 0.) { coeff[Z] = sig1/sig2; }
286 } 370 }
287 371