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