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