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