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