| 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: G4NeutronInelasticXS.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: G4NeutronInelasticXS 33 // File name: G4NeutronInelasticXS
32 // 34 //
33 // Author Ivantchenko, Geant4, 3-Aug-09 35 // Author Ivantchenko, Geant4, 3-Aug-09
34 // 36 //
>> 37 // Modifications:
>> 38 //
35 39
36 #include "G4NeutronInelasticXS.hh" 40 #include "G4NeutronInelasticXS.hh"
37 #include "G4Neutron.hh" 41 #include "G4Neutron.hh"
38 #include "G4DynamicParticle.hh" 42 #include "G4DynamicParticle.hh"
39 #include "G4ElementTable.hh" <<
40 #include "G4Material.hh" <<
41 #include "G4Element.hh" 43 #include "G4Element.hh"
>> 44 #include "G4ElementTable.hh"
42 #include "G4PhysicsLogVector.hh" 45 #include "G4PhysicsLogVector.hh"
43 #include "G4CrossSectionDataSetRegistry.hh" << 46 #include "G4PhysicsVector.hh"
44 #include "G4ComponentGGHadronNucleusXsc.hh" 47 #include "G4ComponentGGHadronNucleusXsc.hh"
45 #include "G4HadronicParameters.hh" << 48 #include "G4HadronNucleonXsc.hh"
>> 49 #include "G4NistManager.hh"
>> 50 #include "G4Proton.hh"
46 #include "Randomize.hh" 51 #include "Randomize.hh"
47 #include "G4Neutron.hh" <<
48 #include "G4SystemOfUnits.hh" <<
49 #include "G4IsotopeList.hh" <<
50 #include "G4NuclearRadii.hh" <<
51 #include "G4AutoLock.hh" <<
52 52
>> 53 #include <iostream>
53 #include <fstream> 54 #include <fstream>
54 #include <sstream> 55 #include <sstream>
55 #include <thread> <<
56 56
57 G4double G4NeutronInelasticXS::coeff[] = {1.0} << 57 // factory
58 G4double G4NeutronInelasticXS::lowcoeff[] = {1 << 58 #include "G4CrossSectionFactory.hh"
59 G4ElementData* G4NeutronInelasticXS::data = nu << 59 //
60 G4String G4NeutronInelasticXS::gDataDirectory << 60 G4_DECLARE_XS_FACTORY(G4NeutronInelasticXS);
61 61
62 static std::once_flag applyOnce; << 62 using namespace std;
63 63
64 namespace << 64 const G4int G4NeutronInelasticXS::amin[] = {
65 { << 65 0,
66 G4Mutex nInelasticXSMutex = G4MUTEX_INITIALI << 66 0, 0, 6, 0,10,12,14,16, 0, 0, //1-10
67 } << 67 0, 0, 0,28, 0, 0, 0,36, 0,40, //11-20
>> 68 0, 0, 0, 0, 0,54, 0,58,63,64, //21-30
>> 69 0,70, 0, 0, 0, 0, 0, 0, 0,90, //31-40
>> 70 0, 0, 0, 0, 0, 0,107,106, 0,112, //41-50
>> 71 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //51-60
>> 72 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //61-70
>> 73 0, 0, 0,180, 0, 0, 0, 0, 0, 0, //71-80
>> 74 0,204, 0, 0, 0, 0, 0, 0, 0, 0, //81-90
>> 75 0,235};
>> 76 const G4int G4NeutronInelasticXS::amax[] = {
>> 77 0,
>> 78 0, 0, 7, 0,11,13,15,18, 0, 0, //1-10
>> 79 0, 0, 0,30, 0, 0, 0,40, 0,48, //11-20
>> 80 0, 0, 0, 0, 0,58, 0,64,65,70, //21-30
>> 81 0,76, 0, 0, 0, 0, 0, 0, 0,96, //31-40
>> 82 0, 0, 0, 0, 0, 0,109,116, 0,124, //41-50
>> 83 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //51-60
>> 84 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //61-70
>> 85 0, 0, 0,186, 0, 0, 0, 0, 0, 0, //71-80
>> 86 0,208, 0, 0, 0, 0, 0, 0, 0, 0, //81-90
>> 87 0,238};
>> 88
>> 89 G4double G4NeutronInelasticXS::coeff[] = {1.0};
>> 90
>> 91 G4ElementData* G4NeutronInelasticXS::data = 0;
68 92
69 G4NeutronInelasticXS::G4NeutronInelasticXS() 93 G4NeutronInelasticXS::G4NeutronInelasticXS()
70 : G4VCrossSectionDataSet(Default_Name()), 94 : G4VCrossSectionDataSet(Default_Name()),
71 neutron(G4Neutron::Neutron()), << 95 proton(G4Proton::Proton())
72 elimit(20*CLHEP::MeV), <<
73 lowElimit(1.0e-7*CLHEP::eV) <<
74 { 96 {
75 verboseLevel = 0; << 97 // verboseLevel = 0;
76 if (verboseLevel > 0) { << 98 if(verboseLevel > 0){
77 G4cout << "G4NeutronInelasticXS::G4Neutron 99 G4cout << "G4NeutronInelasticXS::G4NeutronInelasticXS Initialise for Z < "
78 << MAXZINEL << G4endl; << 100 << MAXZINEL << G4endl;
79 } 101 }
80 loglowElimit = G4Log(lowElimit); << 102 ggXsection = new G4ComponentGGHadronNucleusXsc();
81 if (nullptr == data) { << 103 fNucleon = new G4HadronNucleonXsc();
82 data = new G4ElementData(MAXZINEL); << 104 isMaster = false;
83 data->SetName("nInelastic"); << 105 }
84 FindDirectoryPath(); <<
85 } <<
86 ggXsection = <<
87 G4CrossSectionDataSetRegistry::Instance()- <<
88 if(ggXsection == nullptr) <<
89 ggXsection = new G4ComponentGGHadronNucleu <<
90 106
91 SetForAllAtomsAndEnergies(true); << 107 G4NeutronInelasticXS::~G4NeutronInelasticXS()
>> 108 {
>> 109 //G4cout << "G4NeutronInelasticXS::~G4NeutronInelasticXS() "
>> 110 // << " isMaster= " << isMaster << " data: " << data << G4endl;
>> 111 delete fNucleon;
>> 112 delete ggXsection;
>> 113 if(isMaster) { delete data; data = 0; }
92 } 114 }
93 115
94 void G4NeutronInelasticXS::CrossSectionDescrip 116 void G4NeutronInelasticXS::CrossSectionDescription(std::ostream& outFile) const
95 { 117 {
96 outFile << "G4NeutronInelasticXS calculates 118 outFile << "G4NeutronInelasticXS calculates the neutron inelastic scattering\n"
97 << "cross section on nuclei using da 119 << "cross section on nuclei using data from the high precision\n"
98 << "neutron database. These data ar 120 << "neutron database. These data are simplified and smoothed over\n"
99 << "the resonance region in order to 121 << "the resonance region in order to reduce CPU time.\n"
100 << "For high energy Glauber-Gribov c << 122 << "G4NeutronInelasticXS is valid for energies up to 20 MeV, for\n"
>> 123 << "nuclei through U.\n";
101 } 124 }
102 125
103 G4bool 126 G4bool
104 G4NeutronInelasticXS::IsElementApplicable(cons 127 G4NeutronInelasticXS::IsElementApplicable(const G4DynamicParticle*,
105 G4in << 128 G4int, const G4Material*)
106 { 129 {
107 return true; 130 return true;
108 } 131 }
109 132
110 G4bool 133 G4bool
111 G4NeutronInelasticXS::IsIsoApplicable(const G4 134 G4NeutronInelasticXS::IsIsoApplicable(const G4DynamicParticle*,
112 G4int, G << 135 G4int /*ZZ*/, G4int /*AA*/,
113 const G4 << 136 const G4Element*, const G4Material*)
114 { 137 {
115 return true; 138 return true;
116 } 139 }
117 140
118 G4double << 141 G4double G4NeutronInelasticXS::GetElementCrossSection(
119 G4NeutronInelasticXS::GetElementCrossSection(c << 142 const G4DynamicParticle* aParticle,
120 G << 143 G4int Z, const G4Material*)
121 { 144 {
122 return ElementCrossSection(aParticle->GetKin << 145 G4double xs = 0.0;
123 aParticle->GetLog << 146 G4double ekin = aParticle->GetKineticEnergy();
124 } <<
125 147
126 G4double << 148 if(Z >= MAXZINEL) { Z = MAXZINEL - 1; }
127 G4NeutronInelasticXS::ComputeCrossSectionPerEl << 149 else if(Z < 1) { Z = 1; }
128 <<
129 <<
130 <<
131 { <<
132 return ElementCrossSection(ekin, loge, elm-> <<
133 } <<
134 150
135 G4double << 151 G4int Amean = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
136 G4NeutronInelasticXS::ElementCrossSection(G4do <<
137 { <<
138 G4int Z = std::min(ZZ, MAXZINEL-1); <<
139 G4double ekin = eKin; <<
140 G4double loge = logE; <<
141 G4double xs; <<
142 152
143 // very low energy limit << 153 G4PhysicsVector* pv = data->GetElementData(Z);
144 if (ekin < lowElimit) { << 154 // G4cout << "G4NeutronInelasticXS::GetCrossSection e= " << ekin
145 ekin = lowElimit; << 155 // << " Z= " << Z << G4endl;
146 loge = loglowElimit; << 156
>> 157 // element was not initialised
>> 158 if(!pv) {
>> 159 Initialise(Z);
>> 160 pv = data->GetElementData(Z);
>> 161 if(!pv) { return xs; }
147 } 162 }
148 // pv should exist <<
149 auto pv = GetPhysicsVector(Z); <<
150 163
151 const G4double e0 = pv->Energy(0); << 164 G4double e1 = pv->Energy(0);
152 if (ekin <= e0) { << 165 if(ekin <= e1) { return xs; }
153 xs = (*pv)[0]; << 166
154 if (xs > 0.0) { xs *= std::sqrt(e0/ekin); << 167 G4double e2 = pv->GetMaxEnergy();
155 } else if (ekin <= pv->GetMaxEnergy()) { << 168
156 xs = pv->LogVectorValue(ekin, loge); << 169 if(ekin <= e2) {
157 } else { << 170 xs = pv->Value(ekin);
158 xs = coeff[Z]*ggXsection->GetInelasticElem << 171 } else if(1 == Z) {
159 << 172 fNucleon->GetHadronNucleonXscPDG(aParticle, proton);
>> 173 xs = coeff[1]*fNucleon->GetInelasticHadronNucleonXsc();
>> 174 } else {
>> 175 ggXsection->GetIsoCrossSection(aParticle, Z, Amean);
>> 176 xs = coeff[Z]*ggXsection->GetInelasticGlauberGribovXsc();
160 } 177 }
161 178
162 #ifdef G4VERBOSE << 179 if(verboseLevel > 0) {
163 if(verboseLevel > 1) { << 180 G4cout << "ekin= " << ekin << ", XSinel= " << xs << G4endl;
164 G4cout << "G4NeutronInelasticXS::ElementC <<
165 << " Ekin(MeV)= " << ekin/CLHEP::M <<
166 << ", ElmXSinel(b)= " << xs/CLHEP: <<
167 << G4endl; <<
168 } 181 }
169 #endif <<
170 return xs; 182 return xs;
171 } 183 }
172 184
173 G4double << 185 G4double G4NeutronInelasticXS::GetIsoCrossSection(
174 G4NeutronInelasticXS::ComputeIsoCrossSection(G << 186 const G4DynamicParticle* aParticle,
175 c << 187 G4int Z, G4int A,
176 G << 188 const G4Isotope*, const G4Element*,
177 c << 189 const G4Material*)
178 c << 190 {
179 { << 191 if(Z >= MAXZINEL) { Z = MAXZINEL - 1; }
180 return IsoCrossSection(ekin, loge, Z, A); << 192 else if(Z < 1) { Z = 1; }
181 } <<
182 <<
183 G4double <<
184 G4NeutronInelasticXS::GetIsoCrossSection(const <<
185 G4int <<
186 const <<
187 const <<
188 { <<
189 return IsoCrossSection(aParticle->GetKinetic <<
190 aParticle->GetLogKine <<
191 } <<
192 <<
193 G4double <<
194 G4NeutronInelasticXS::IsoCrossSection(G4double <<
195 G4int ZZ <<
196 { <<
197 G4double xs; <<
198 G4int Z = std::min(ZZ, MAXZINEL-1); <<
199 G4double ekin = eKin; <<
200 G4double loge = logE; <<
201 <<
202 /* <<
203 G4cout << "G4NeutronInelasticXS::IsoCrossSec <<
204 << Z << " A= " << A << G4endl; <<
205 G4cout << " Amin= " << amin[Z] << " Amax= " <<
206 << " E(MeV)= " << ekin << " Ncomp=" <<
207 << data->GetNumberOfComponents(Z) << <<
208 */ <<
209 GetPhysicsVector(Z); <<
210 193
211 // use isotope cross section if applicable << 194 return IsoCrossSection(aParticle->GetKineticEnergy(), Z, A);
212 if (ekin <= elimit && data->GetNumberOfCompo << 195 }
213 auto pviso = data->GetComponentDataByID(Z, <<
214 if (nullptr != pviso) { <<
215 const G4double e0 = pviso->Energy(0); <<
216 if (ekin > e0) { <<
217 xs = pviso->LogVectorValue(ekin, loge) <<
218 } else { <<
219 xs = (*pviso)[0]; <<
220 if (xs > 0.0) { xs *= std::sqrt(e0/eki <<
221 } <<
222 #ifdef G4VERBOSE <<
223 if(verboseLevel > 1) { <<
224 G4cout << "G4NeutronInelasticXS::IsoXS <<
225 << ekin/CLHEP::MeV <<
226 << " xs(b)= " << xs/CLHEP::bar <<
227 << " Z= " << Z << " A= " << A <<
228 } <<
229 #endif <<
230 return xs; <<
231 } <<
232 } <<
233 196
234 // use element x-section << 197 G4double
235 xs = ElementCrossSection(ekin, loge, Z)*A/ae << 198 G4NeutronInelasticXS::IsoCrossSection(G4double ekin, G4int Z, G4int A)
>> 199 {
>> 200 G4double xs = 0.0;
>> 201 /*
>> 202 G4cout << "IsoCrossSection Z= " << Z << " A= " << A
>> 203 << " Amin= " << amin[Z] << " Amax= " << amax[Z]
>> 204 << " E(MeV)= " << ekin << G4endl;
>> 205 */
>> 206 // element was not initialised
>> 207 G4PhysicsVector* pv = data->GetElementData(Z);
>> 208 if(!pv) {
>> 209 Initialise(Z);
>> 210 pv = data->GetElementData(Z);
>> 211 }
236 212
237 #ifdef G4VERBOSE << 213 // isotope cross section exist
238 if(verboseLevel > 1) { << 214 if(pv && amin[Z] > 0 && A >= amin[Z] && A <= amax[Z]) {
239 G4cout << "G4NeutronInelasticXS::IsoXS: Z << 215 pv = data->GetComponentDataByIndex(Z, A - amin[Z]);
240 << " Ekin(MeV)= " << ekin/CLHEP::M << 216 if(pv && ekin > pv->Energy(0)) { xs = pv->Value(ekin); }
241 << ", ElmXS(b)= " << xs/CLHEP::bar << 217 }
>> 218 if(verboseLevel > 0){
>> 219 G4cout << "ekin= " << ekin << ", xs= " << xs << G4endl;
242 } 220 }
243 #endif <<
244 return xs; 221 return xs;
245 } 222 }
246 223
247 const G4Isotope* G4NeutronInelasticXS::SelectI << 224 G4Isotope* G4NeutronInelasticXS::SelectIsotope(
248 const G4Element* anElement, G4double kin << 225 const G4Element* anElement, G4double kinEnergy)
249 { 226 {
250 std::size_t nIso = anElement->GetNumberOfIso << 227 size_t nIso = anElement->GetNumberOfIsotopes();
251 const G4Isotope* iso = anElement->GetIsotope << 228 G4IsotopeVector* isoVector = anElement->GetIsotopeVector();
252 if(1 == nIso) { return iso; } << 229 G4Isotope* iso = (*isoVector)[0];
253 230
254 // more than 1 isotope << 231 //G4cout << "SelectIsotope NIso= " << nIso << G4endl;
255 G4int Z = anElement->GetZasInt(); <<
256 if (nullptr == data->GetElementData(Z)) { In <<
257 232
258 const G4double* abundVector = anElement->Get << 233 // more than 1 isotope
259 G4double q = G4UniformRand(); << 234 if(1 < nIso) {
260 G4double sum = 0.0; << 235 G4int Z = G4lrint(anElement->GetZ());
261 std::size_t j; << 236 //G4cout << "SelectIsotope Z= " << Z << G4endl;
262 << 237
263 // isotope wise cross section not available << 238 G4double* abundVector = anElement->GetRelativeAbundanceVector();
264 if (Z >= MAXZINEL || 0 == data->GetNumberOfC << 239 G4double q = G4UniformRand();
265 for (j=0; j<nIso; ++j) { << 240 G4double sum = 0.0;
266 sum += abundVector[j]; << 241
267 if(q <= sum) { << 242 // is there isotope wise cross section?
268 iso = anElement->GetIsotope((G4int)j); << 243 size_t j;
269 break; << 244 if(0 == amin[Z] || Z >= MAXZINEL) {
>> 245 for (j = 0; j<nIso; ++j) {
>> 246 sum += abundVector[j];
>> 247 if(q <= sum) {
>> 248 iso = (*isoVector)[j];
>> 249 break;
>> 250 }
270 } 251 }
271 } << 252 } else {
272 return iso; <<
273 } <<
274 253
275 // use isotope cross sections << 254 // element may be not initialised in unit test
276 auto nn = temp.size(); << 255 if(!data->GetElementData(Z)) { Initialise(Z); }
277 if(nn < nIso) { temp.resize(nIso, 0.); } << 256 size_t nn = temp.size();
278 << 257 if(nn < nIso) { temp.resize(nIso, 0.); }
279 for (j=0; j<nIso; ++j) { << 258
280 // G4cout << j << "-th isotope " << anElem << 259 for (j=0; j<nIso; ++j) {
281 // << " abund= " << abundVector[j] << 260 //G4cout << j << "-th isotope " << (*isoVector)[j]->GetN()
282 sum += abundVector[j]*IsoCrossSection(kinE << 261 // << " abund= " << abundVector[j] << G4endl;
283 anEl << 262 sum += abundVector[j]*IsoCrossSection(kinEnergy, Z,
284 temp[j] = sum; << 263 (*isoVector)[j]->GetN());
285 } << 264 temp[j] = sum;
286 sum *= q; << 265 }
287 for (j = 0; j<nIso; ++j) { << 266 sum *= q;
288 if (temp[j] >= sum) { << 267 for (j = 0; j<nIso; ++j) {
289 iso = anElement->GetIsotope((G4int)j); << 268 if(temp[j] >= sum) {
290 break; << 269 iso = (*isoVector)[j];
>> 270 break;
>> 271 }
>> 272 }
291 } 273 }
292 } 274 }
293 return iso; 275 return iso;
294 } 276 }
295 277
296 void 278 void
297 G4NeutronInelasticXS::BuildPhysicsTable(const 279 G4NeutronInelasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
298 { 280 {
299 if (verboseLevel > 0) { << 281 if(verboseLevel > 0){
300 G4cout << "G4NeutronInelasticXS::BuildPhys 282 G4cout << "G4NeutronInelasticXS::BuildPhysicsTable for "
301 << p.GetParticleName() << G4endl; << 283 << p.GetParticleName() << G4endl;
302 } 284 }
303 if (p.GetParticleName() != "neutron") { << 285 if(p.GetParticleName() != "neutron") {
304 G4ExceptionDescription ed; 286 G4ExceptionDescription ed;
305 ed << p.GetParticleName() << " is a wrong 287 ed << p.GetParticleName() << " is a wrong particle type -"
306 << " only neutron is allowed"; 288 << " only neutron is allowed";
307 G4Exception("G4NeutronInelasticXS::BuildPh 289 G4Exception("G4NeutronInelasticXS::BuildPhysicsTable(..)","had012",
308 FatalException, ed, ""); << 290 FatalException, ed, "");
309 return; 291 return;
310 } 292 }
311 // it is possible re-initialisation for the <<
312 const G4ElementTable* table = G4Element::Get <<
313 293
314 // initialise static tables only once << 294 if(!data) {
315 std::call_once(applyOnce, [this]() { isIniti << 295 isMaster = true;
>> 296 data = new G4ElementData();
>> 297 data->SetName("NeutronInelastic");
>> 298 temp.resize(13,0.0);
>> 299 }
316 300
317 if (isInitializer) { << 301 // it is possible re-initialisation for the new run
318 G4AutoLock l(&nInelasticXSMutex); << 302 if(isMaster) {
319 303
320 // Upload data for elements used in geomet << 304 // check environment variable
321 for ( auto const & elm : *table ) { << 305 // Build the complete string identifying the file with the data set
322 G4int Z = std::max( 1, std::min( elm->Ge << 306 char* path = getenv("G4NEUTRONXSDATA");
323 if ( nullptr == data->GetElementData(Z) << 307
>> 308 G4DynamicParticle* dynParticle =
>> 309 new G4DynamicParticle(G4Neutron::Neutron(),G4ThreeVector(1,0,0),1);
>> 310
>> 311 // Access to elements
>> 312 const G4ElementTable* theElmTable = G4Element::GetElementTable();
>> 313 size_t numOfElm = G4Element::GetNumberOfElements();
>> 314 if(numOfElm > 0) {
>> 315 for(size_t i=0; i<numOfElm; ++i) {
>> 316 G4int Z = G4lrint(((*theElmTable)[i])->GetZ());
>> 317 if(Z < 1) { Z = 1; }
>> 318 else if(Z >= MAXZINEL) { Z = MAXZINEL-1; }
>> 319 //G4cout << "Z= " << Z << G4endl;
>> 320 // Initialisation
>> 321 if(!(data->GetElementData(Z))) {
>> 322 Initialise(Z, dynParticle, path);
>> 323 }
>> 324 }
324 } 325 }
325 l.unlock(); << 326 delete dynParticle;
326 } 327 }
327 // prepare isotope selection <<
328 std::size_t nIso = temp.size(); <<
329 for ( auto const & elm : *table ) { <<
330 std::size_t n = elm->GetNumberOfIsotopes() <<
331 if (n > nIso) { nIso = n; } <<
332 } <<
333 temp.resize(nIso, 0.0); <<
334 } 328 }
335 329
336 const G4String& G4NeutronInelasticXS::FindDire << 330 void
>> 331 G4NeutronInelasticXS::Initialise(G4int Z, G4DynamicParticle* dp,
>> 332 const char* p)
337 { 333 {
338 // build the complete string identifying the << 334 if(data->GetElementData(Z) || Z < 1 || Z >= MAXZINEL) { return; }
339 if (gDataDirectory.empty()) { << 335 const char* path = p;
340 std::ostringstream ost; << 336 if(!p) {
341 ost << G4HadronicParameters::Instance()->G << 337 // check environment variable
342 gDataDirectory = ost.str(); << 338 // Build the complete string identifying the file with the data set
>> 339 path = getenv("G4NEUTRONXSDATA");
>> 340 if (!path) {
>> 341 G4Exception("G4NeutronInelasticXS::Initialise(..)","had013",
>> 342 FatalException,
>> 343 "Environment variable G4NEUTRONXSDATA is not defined");
>> 344 return;
>> 345 }
>> 346 }
>> 347 G4DynamicParticle* dynParticle = dp;
>> 348 if(!dp) {
>> 349 dynParticle =
>> 350 new G4DynamicParticle(G4Neutron::Neutron(),G4ThreeVector(1,0,0),1);
343 } 351 }
344 return gDataDirectory; <<
345 } <<
346 <<
347 void G4NeutronInelasticXS::InitialiseOnFly(G4i <<
348 { <<
349 G4AutoLock l(&nInelasticXSMutex); <<
350 Initialise(Z); <<
351 l.unlock(); <<
352 } <<
353 352
354 void G4NeutronInelasticXS::Initialise(G4int Z) << 353 G4int Amean = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
355 { <<
356 if (nullptr != data->GetElementData(Z)) { re <<
357 354
358 // upload element data 355 // upload element data
359 std::ostringstream ost; 356 std::ostringstream ost;
360 ost << FindDirectoryPath() << Z; << 357 ost << path << "/inelast" << Z ;
361 G4PhysicsVector* v = RetrieveVector(ost, tru 358 G4PhysicsVector* v = RetrieveVector(ost, true);
362 data->InitialiseForElement(Z, v); 359 data->InitialiseForElement(Z, v);
363 if (verboseLevel > 1) { << 360 /*
364 G4cout << "G4NeutronInelasticXS::Initiali << 361 G4cout << "G4NeutronInelasticXS::Initialise for Z= " << Z
365 << " A= " << aeff[Z] << " Amin= " << 362 << " A= " << Amean << " Amin= " << amin[Z]
366 << " Amax= " << amax[Z] << G4endl << 363 << " Amax= " << amax[Z] << G4endl;
367 } << 364 */
368 // upload isotope data 365 // upload isotope data
369 G4bool noComp = true; << 366 if(amin[Z] > 0) {
370 if (amin[Z] < amax[Z]) { << 367 size_t nmax = (size_t)(amax[Z]-amin[Z]+1);
>> 368 data->InitialiseForComponent(Z, nmax);
371 369
372 for (G4int A=amin[Z]; A<=amax[Z]; ++A) { << 370 for(G4int A=amin[Z]; A<=amax[Z]; ++A) {
373 std::ostringstream ost1; 371 std::ostringstream ost1;
374 ost1 << gDataDirectory << Z << "_" << A; << 372 ost1 << path << "/inelast" << Z << "_" << A;
375 G4PhysicsVector* v1 = RetrieveVector(ost 373 G4PhysicsVector* v1 = RetrieveVector(ost1, false);
376 if (nullptr != v1) { << 374 data->AddComponent(Z, A, v1);
377 if (noComp) { <<
378 G4int nmax = amax[Z] - A + 1; <<
379 data->InitialiseForComponent(Z, nmax <<
380 noComp = false; <<
381 } <<
382 data->AddComponent(Z, A, v1); <<
383 } <<
384 } 375 }
385 } 376 }
386 // no components case <<
387 if (noComp) { data->InitialiseForComponent(Z <<
388 377
389 // smooth transition 378 // smooth transition
390 G4double sig1 = (*v)[v->GetVectorLength()-1] << 379 G4double emax = v->GetMaxEnergy();
391 G4double ehigh= v->GetMaxEnergy(); << 380 G4double sig1 = (*v)[v->GetVectorLength() - 1];
392 G4double sig2 = ggXsection->GetInelasticElem << 381 dynParticle->SetKineticEnergy(emax);
393 ehigh, Z, aeff[Z << 382 G4double sig2 = 0.0;
394 coeff[Z] = (sig2 > 0.) ? sig1/sig2 : 1.0; << 383 if(1 == Z) {
>> 384 fNucleon->GetHadronNucleonXscPDG(dynParticle, proton);
>> 385 sig2 = fNucleon->GetInelasticHadronNucleonXsc();
>> 386 } else {
>> 387 ggXsection->GetIsoCrossSection(dynParticle, Z, Amean);
>> 388 sig2 = ggXsection->GetInelasticGlauberGribovXsc();
>> 389 }
>> 390 if(sig2 > 0.) { coeff[Z] = sig1/sig2; }
>> 391 if(!dp) { delete dynParticle; }
395 } 392 }
396 393
397 G4PhysicsVector* 394 G4PhysicsVector*
398 G4NeutronInelasticXS::RetrieveVector(std::ostr 395 G4NeutronInelasticXS::RetrieveVector(std::ostringstream& ost, G4bool warn)
399 { 396 {
400 G4PhysicsLogVector* v = nullptr; << 397 G4PhysicsLogVector* v = 0;
401 std::ifstream filein(ost.str().c_str()); 398 std::ifstream filein(ost.str().c_str());
402 if (!filein.is_open()) { << 399 if (!(filein)) {
403 if(warn) { 400 if(warn) {
404 G4ExceptionDescription ed; 401 G4ExceptionDescription ed;
405 ed << "Data file <" << ost.str().c_str() 402 ed << "Data file <" << ost.str().c_str()
406 << "> is not opened!"; << 403 << "> is not opened!";
407 G4Exception("G4NeutronInelasticXS::Retri 404 G4Exception("G4NeutronInelasticXS::RetrieveVector(..)","had014",
408 FatalException, ed, "Check G << 405 FatalException, ed, "Check G4NEUTRONXSDATA");
409 } 406 }
410 } else { 407 } else {
411 if(verboseLevel > 1) { 408 if(verboseLevel > 1) {
412 G4cout << "File " << ost.str() 409 G4cout << "File " << ost.str()
413 << " is opened by G4NeutronInelas << 410 << " is opened by G4NeutronInelasticXS" << G4endl;
414 } 411 }
415 // retrieve data from DB 412 // retrieve data from DB
416 v = new G4PhysicsLogVector(); 413 v = new G4PhysicsLogVector();
417 if(!v->Retrieve(filein, true)) { 414 if(!v->Retrieve(filein, true)) {
418 G4ExceptionDescription ed; 415 G4ExceptionDescription ed;
419 ed << "Data file <" << ost.str().c_str() 416 ed << "Data file <" << ost.str().c_str()
420 << "> is not retrieved!"; << 417 << "> is not retrieved!";
421 G4Exception("G4NeutronInelasticXS::Retri 418 G4Exception("G4NeutronInelasticXS::RetrieveVector(..)","had015",
422 FatalException, ed, "Check G << 419 FatalException, ed, "Check G4NEUTRONXSDATA");
423 } 420 }
424 } 421 }
425 return v; 422 return v;
426 } 423 }
427 424