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