| 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 // V. Ivanchenko September 2023 26 // V. Ivanchenko September 2023
27 // 27 //
28 // G4CrossSectionHP is a generic class impleme 28 // G4CrossSectionHP is a generic class implementing
29 // cross sections for neutrons, protons and li 29 // cross sections for neutrons, protons and light ions
30 // It is an alternative to code developed by J 30 // It is an alternative to code developed by J.P. Wellisch & T.Koi
31 // 31 //
32 32
33 #include <fstream> 33 #include <fstream>
34 #include <sstream> 34 #include <sstream>
35 #include <thread> 35 #include <thread>
36 36
37 #include "G4CrossSectionHP.hh" 37 #include "G4CrossSectionHP.hh"
38 #include "G4Material.hh" 38 #include "G4Material.hh"
39 #include "G4Element.hh" 39 #include "G4Element.hh"
40 #include "G4ElementDataRegistry.hh" 40 #include "G4ElementDataRegistry.hh"
41 #include "G4ElementTable.hh" 41 #include "G4ElementTable.hh"
42 #include "G4IsotopeList.hh" 42 #include "G4IsotopeList.hh"
43 #include "G4HadronicParameters.hh" 43 #include "G4HadronicParameters.hh"
44 #include "G4ParticleHPManager.hh" 44 #include "G4ParticleHPManager.hh"
45 #include "G4ParticleDefinition.hh" 45 #include "G4ParticleDefinition.hh"
46 #include "G4PhysicalConstants.hh" 46 #include "G4PhysicalConstants.hh"
47 #include "G4Pow.hh" 47 #include "G4Pow.hh"
48 #include "G4Log.hh" 48 #include "G4Log.hh"
49 #include "Randomize.hh" 49 #include "Randomize.hh"
50 #include "G4RandomDirection.hh" 50 #include "G4RandomDirection.hh"
51 #include "G4SystemOfUnits.hh" 51 #include "G4SystemOfUnits.hh"
52 #include "G4Neutron.hh" 52 #include "G4Neutron.hh"
53 #include "G4NucleiProperties.hh" 53 #include "G4NucleiProperties.hh"
54 #include "G4AutoLock.hh" 54 #include "G4AutoLock.hh"
55 55
56 namespace 56 namespace
57 { 57 {
58 G4Mutex theHPXS = G4MUTEX_INITIALIZER; 58 G4Mutex theHPXS = G4MUTEX_INITIALIZER;
59 } 59 }
60 60
61 G4CrossSectionHP::G4CrossSectionHP(const G4Par 61 G4CrossSectionHP::G4CrossSectionHP(const G4ParticleDefinition* p,
62 const G4Str 62 const G4String& nameData,
63 const G4Str 63 const G4String& nameDir, G4double emaxHP,
64 G4int zmin, 64 G4int zmin, G4int zmax)
65 : G4VCrossSectionDataSet(nameData), 65 : G4VCrossSectionDataSet(nameData),
66 fParticle(p), 66 fParticle(p),
67 fNeutron(G4Neutron::Neutron()), <<
68 fManagerHP(G4ParticleHPManager::GetInstanc 67 fManagerHP(G4ParticleHPManager::GetInstance()),
69 emax(emaxHP), 68 emax(emaxHP),
70 emaxT(fManagerHP->GetMaxEnergyDoppler()), 69 emaxT(fManagerHP->GetMaxEnergyDoppler()),
71 elimit(1.0e-11*CLHEP::eV), 70 elimit(1.0e-11*CLHEP::eV),
72 logElimit(G4Log(elimit)), 71 logElimit(G4Log(elimit)),
73 minZ(zmin), 72 minZ(zmin),
74 maxZ(zmax), 73 maxZ(zmax),
75 fDataName(nameData), 74 fDataName(nameData),
76 fDataDirectory(nameDir) 75 fDataDirectory(nameDir)
77 { 76 {
78 // verboseLevel = 1; 77 // verboseLevel = 1;
79 if (verboseLevel > 1) { 78 if (verboseLevel > 1) {
80 G4cout << "G4CrossSectionHP::G4CrossSectio 79 G4cout << "G4CrossSectionHP::G4CrossSectionHP: Initialise for "
81 << fDataName << " " << minZ << " < Z < " 80 << fDataName << " " << minZ << " < Z < " << maxZ
82 << " EmaxT(MeV)=" << emaxT << G4endl; 81 << " EmaxT(MeV)=" << emaxT << G4endl;
83 G4cout << "Data directory: " << fDataDirec 82 G4cout << "Data directory: " << fDataDirectory << G4endl;
84 } 83 }
85 auto ptr = G4ElementDataRegistry::Instance() 84 auto ptr = G4ElementDataRegistry::Instance();
86 auto data = ptr->GetElementDataByName(fDataN 85 auto data = ptr->GetElementDataByName(fDataName);
87 if (nullptr == data) { 86 if (nullptr == data) {
88 data = new G4ElementData(maxZ - minZ + 1); 87 data = new G4ElementData(maxZ - minZ + 1);
89 data->SetName(fDataName); 88 data->SetName(fDataName);
90 } 89 }
91 fData = data; 90 fData = data;
92 } 91 }
93 92
94 G4bool G4CrossSectionHP::IsIsoApplicable(const 93 G4bool G4CrossSectionHP::IsIsoApplicable(const G4DynamicParticle* dp,
95 G4int, G4int, 94 G4int, G4int,
96 const 95 const G4Element*,
97 const G4Material*) 96 const G4Material*)
98 { 97 {
99 return (dp->GetKineticEnergy() <= emax); 98 return (dp->GetKineticEnergy() <= emax);
100 } 99 }
101 100
102 G4double G4CrossSectionHP::GetIsoCrossSection( 101 G4double G4CrossSectionHP::GetIsoCrossSection(const G4DynamicParticle* dp,
103 102 G4int Z, G4int A,
104 103 const G4Isotope*,
105 const G4Element*, 104 const G4Element*,
106 const G4Material* mat) 105 const G4Material* mat)
107 { 106 {
108 G4double ekin = dp->GetKineticEnergy(); 107 G4double ekin = dp->GetKineticEnergy();
109 G4double loge = dp->GetLogKineticEnergy(); 108 G4double loge = dp->GetLogKineticEnergy();
110 if (ekin < elimit) { 109 if (ekin < elimit) {
111 ekin = elimit; 110 ekin = elimit;
112 loge = logElimit; 111 loge = logElimit;
113 } 112 }
114 if (mat != fCurrentMat) { PrepareCache(mat); 113 if (mat != fCurrentMat) { PrepareCache(mat); }
115 114
116 return IsoCrossSection(ekin, loge, Z, A, mat 115 return IsoCrossSection(ekin, loge, Z, A, mat->GetTemperature());
117 } 116 }
118 117
119 G4double 118 G4double
120 G4CrossSectionHP::ComputeIsoCrossSection(G4dou 119 G4CrossSectionHP::ComputeIsoCrossSection(G4double e, G4double le,
121 const 120 const G4ParticleDefinition*,
122 G4int 121 G4int Z, G4int A,
123 const G4Isotope*, 122 const G4Isotope*,
124 const G4Element*, 123 const G4Element*,
125 const G4Material* mat) 124 const G4Material* mat)
126 { 125 {
127 G4double ekin = e; 126 G4double ekin = e;
128 G4double loge = le; 127 G4double loge = le;
129 if (ekin < elimit) { 128 if (ekin < elimit) {
130 ekin = elimit; 129 ekin = elimit;
131 loge = logElimit; 130 loge = logElimit;
132 } 131 }
133 if (mat != fCurrentMat) { PrepareCache(mat); 132 if (mat != fCurrentMat) { PrepareCache(mat); }
134 133
135 return IsoCrossSection(ekin, loge, Z, A, mat 134 return IsoCrossSection(ekin, loge, Z, A, mat->GetTemperature());
136 } 135 }
137 136
138 G4double G4CrossSectionHP::IsoCrossSection(con 137 G4double G4CrossSectionHP::IsoCrossSection(const G4double ekin,
139 con 138 const G4double logek,
140 const G4int Z, const G4int A, 139 const G4int Z, const G4int A,
141 con 140 const G4double T)
142 { 141 {
143 // G4cout << "G4CrossSectionHP::IsoCrossSect << 142 //G4cout << "G4CrossSectionHP::IsoCrossSection Z=" << Z << " A=" << A
144 // << " E(MeV)=" << ekin/MeV << " T=" << T < 143 // << " E(MeV)=" << ekin/MeV << " T=" << T << " " << GetName() << G4endl;
145 G4double xs = 0.0; 144 G4double xs = 0.0;
146 if (ekin > emax || Z > maxZ || Z < minZ || e 145 if (ekin > emax || Z > maxZ || Z < minZ || ekin < elimit) { return xs; }
147 146
148 auto pv0 = fData->GetElementData(Z - minZ); << 147 const G4PhysicsVector* pv0 = fData->GetElementData(Z - minZ);
149 if (nullptr == pv0) { 148 if (nullptr == pv0) {
150 Initialise(Z); << 149 InitialiseOnFly(Z);
151 pv0 = fData->GetElementData(Z - minZ); 150 pv0 = fData->GetElementData(Z - minZ);
152 } 151 }
153 <<
154 // if there is no element data then no iso d <<
155 if (nullptr == pv0) { return xs; } 152 if (nullptr == pv0) { return xs; }
156 << 153 const G4PhysicsVector* pv = fData->GetComponentDataByID(Z - minZ, A);
157 const auto pv = fData->GetComponentDataByID( <<
158 if (nullptr == pv) { return xs; } 154 if (nullptr == pv) { return xs; }
159 155
160 // no Doppler broading 156 // no Doppler broading
161 // G4double factT = T/CLHEP::STP_Temperature << 157 G4double factT = T/CLHEP::STP_Temperature;
162 if (ekin >= emaxT*T/CLHEP::STP_Temperature | << 158 if (ekin >= emaxT*factT || fManagerHP->GetNeglectDoppler()) {
163 xs = pv->LogFreeVectorValue(ekin, logek); 159 xs = pv->LogFreeVectorValue(ekin, logek);
164 160
165 } else { 161 } else {
166 162
167 // Doppler broading 163 // Doppler broading
168 G4double e0 = CLHEP::k_Boltzmann*T; << 164 G4double lambda = 1.0/(CLHEP::k_Boltzmann*T);
169 G4double mass = fParticle->GetPDGMass(); 165 G4double mass = fParticle->GetPDGMass();
170 G4double massTarget = G4NucleiProperties:: 166 G4double massTarget = G4NucleiProperties::GetNuclearMass(A, Z);
171 G4double sig = std::sqrt(2.0*e0/(3.0*massT << 167 G4LorentzVector lv(0., 0., 0., mass + ekin);
172 <<
173 // projectile <<
174 G4LorentzVector lv(0., 0., std::sqrt(ekin* <<
175 168
176 // limits of integration 169 // limits of integration
177 const G4double lim = 1.01; 170 const G4double lim = 1.01;
178 const G4int nmin = 3; 171 const G4int nmin = 3;
179 G4int ii = 0; << 172 G4int i;
180 const G4int nn = 20; 173 const G4int nn = 20;
181 G4double xs2 = 0.0; 174 G4double xs2 = 0.0;
182 175
183 for (G4int i=0; i<nn; ++i) { << 176 for (i=1; i<nn; ++i) {
184 G4double vx = G4RandGauss::shoot(0., sig << 177 G4double erand = G4RandGamma::shoot(2.0, lambda);
185 G4double vy = G4RandGauss::shoot(0., sig << 178 auto mom = G4RandomDirection()*std::sqrt(2*massTarget*erand);
186 G4double vz = G4RandGauss::shoot(0., sig << 179 fLV.set(mom.x(), mom.y(), mom.z(), mass + erand);
187 fLV.set(massTarget*vx, massTarget*vy, ma <<
188 fBoost = fLV.boostVector(); 180 fBoost = fLV.boostVector();
189 fLV = lv.boost(-fBoost); << 181 G4double e = lv.boost(fBoost).e() - mass;
190 if (fLV.pz() <= 0.0) { continue; } <<
191 ++ii; <<
192 G4double e = fLV.e() - mass; <<
193 G4double y = pv->Value(e, index); 182 G4double y = pv->Value(e, index);
194 xs += y; 183 xs += y;
195 xs2 += y*y; 184 xs2 += y*y;
196 if (ii >= nmin && ii*xs2 <= lim*xs*xs) { << 185 if (i >= nmin && i*xs2 <= lim*xs*xs) { break; }
197 } 186 }
198 if (ii > 0) { xs /= (G4double)(ii); } << 187 xs /= (G4double)std::min(i, nn-1);
199 } 188 }
200 #ifdef G4VERBOSE 189 #ifdef G4VERBOSE
201 if (verboseLevel > 1) { 190 if (verboseLevel > 1) {
202 G4cout << "G4CrossSectionHP::IsoXS " << fD 191 G4cout << "G4CrossSectionHP::IsoXS " << fDataName
203 << " Z= " << Z << " A= " << A 192 << " Z= " << Z << " A= " << A
204 << " Ekin(MeV)= " << ekin/MeV << " xs(b) 193 << " Ekin(MeV)= " << ekin/MeV << " xs(b)= " << xs/barn
205 << " size=" << fZA.size() << G4end 194 << " size=" << fZA.size() << G4endl;
206 } 195 }
207 #endif 196 #endif
208 197
209 // save cross section into struct 198 // save cross section into struct
210 for (std::size_t i=0; i<fZA.size(); ++i) { 199 for (std::size_t i=0; i<fZA.size(); ++i) {
211 if (Z == fZA[i].first && A == fZA[i].secon 200 if (Z == fZA[i].first && A == fZA[i].second) {
212 fIsoXS[i] = xs; 201 fIsoXS[i] = xs;
213 break; 202 break;
214 } 203 }
215 } 204 }
216 return xs; 205 return xs;
217 } 206 }
218 207
219 const G4Isotope* G4CrossSectionHP::SelectIsoto 208 const G4Isotope* G4CrossSectionHP::SelectIsotope(const G4Element* elm,
220 209 G4double, G4double)
221 { 210 {
222 std::size_t nIso = elm->GetNumberOfIsotopes( 211 std::size_t nIso = elm->GetNumberOfIsotopes();
223 const G4Isotope* iso = elm->GetIsotope(0); 212 const G4Isotope* iso = elm->GetIsotope(0);
224 213
225 //G4cout << "SelectIsotope NIso= " << nIso < 214 //G4cout << "SelectIsotope NIso= " << nIso << G4endl;
226 if(1 == nIso) { return iso; } 215 if(1 == nIso) { return iso; }
227 216
228 // more than 1 isotope 217 // more than 1 isotope
229 G4int Z = elm->GetZasInt(); 218 G4int Z = elm->GetZasInt();
230 if (Z >= minZ && Z <= maxZ && nullptr == fDa 219 if (Z >= minZ && Z <= maxZ && nullptr == fData->GetElementData(Z - minZ)) {
231 Initialise(Z); << 220 InitialiseOnFly(Z);
232 } 221 }
233 222
234 const G4double* abundVector = elm->GetRelati 223 const G4double* abundVector = elm->GetRelativeAbundanceVector();
235 G4double q = G4UniformRand(); 224 G4double q = G4UniformRand();
236 G4double sum = 0.0; 225 G4double sum = 0.0;
237 226
238 // is there cached isotope wise cross sectio 227 // is there cached isotope wise cross section?
239 std::size_t j; 228 std::size_t j;
240 if (Z < minZ || Z > maxZ || !CheckCache(Z) | 229 if (Z < minZ || Z > maxZ || !CheckCache(Z) ||
241 0 == fData->GetNumberOfComponents(Z - mi 230 0 == fData->GetNumberOfComponents(Z - minZ)) {
242 for (j = 0; j<nIso; ++j) { 231 for (j = 0; j<nIso; ++j) {
243 sum += abundVector[j]; 232 sum += abundVector[j];
244 if(q <= sum) { 233 if(q <= sum) {
245 iso = elm->GetIsotope((G4int)j); 234 iso = elm->GetIsotope((G4int)j);
246 break; 235 break;
247 } 236 }
248 } 237 }
249 return iso; 238 return iso;
250 } 239 }
251 std::size_t nn = fTemp.size(); 240 std::size_t nn = fTemp.size();
252 if (nn < nIso) { fTemp.resize(nIso, 0.); } 241 if (nn < nIso) { fTemp.resize(nIso, 0.); }
253 242
254 // reuse cache 243 // reuse cache
255 for (j=0; j<nIso; ++j) { 244 for (j=0; j<nIso; ++j) {
256 sum += abundVector[j]* 245 sum += abundVector[j]*
257 GetCrossSection(Z - minZ, elm->GetIsotop 246 GetCrossSection(Z - minZ, elm->GetIsotope((G4int)j)->GetN());
258 fTemp[j] = sum; 247 fTemp[j] = sum;
259 } 248 }
260 sum *= q; 249 sum *= q;
261 for (j = 0; j<nIso; ++j) { 250 for (j = 0; j<nIso; ++j) {
262 if (fTemp[j] >= sum) { 251 if (fTemp[j] >= sum) {
263 iso = elm->GetIsotope((G4int)j); 252 iso = elm->GetIsotope((G4int)j);
264 break; 253 break;
265 } 254 }
266 } 255 }
267 return iso; 256 return iso;
268 } 257 }
269 258
270 void G4CrossSectionHP::BuildPhysicsTable(const 259 void G4CrossSectionHP::BuildPhysicsTable(const G4ParticleDefinition& p)
271 { 260 {
272 if (verboseLevel > 1) { << 261 if (verboseLevel > 1){
273 G4cout << "G4CrossSectionHP::BuildPhysicsT 262 G4cout << "G4CrossSectionHP::BuildPhysicsTable for "
274 << p.GetParticleName() << " and " << fDat 263 << p.GetParticleName() << " and " << fDataName << G4endl;
275 } 264 }
276 265
277 // it is possible re-initialisation for the 266 // it is possible re-initialisation for the second run
278 const G4ElementTable* table = G4Element::Get 267 const G4ElementTable* table = G4Element::GetElementTable();
279 268
280 // Access to elements 269 // Access to elements
281 for ( auto const & elm : *table ) { 270 for ( auto const & elm : *table ) {
282 G4int Z = elm->GetZasInt(); 271 G4int Z = elm->GetZasInt();
283 if (Z >= minZ && Z <= maxZ && nullptr == f << 272 if (Z >= minZ && Z <= maxZ &&
284 Initialise(Z); << 273 nullptr == fData->GetElementData(Z - minZ)) {
>> 274 InitialiseOnFly(Z);
285 } 275 }
286 } 276 }
287 277
288 // prepare isotope selection 278 // prepare isotope selection
289 std::size_t nmax = 0; 279 std::size_t nmax = 0;
290 std::size_t imax = 0; 280 std::size_t imax = 0;
291 for ( auto const & mat : *G4Material::GetMat 281 for ( auto const & mat : *G4Material::GetMaterialTable() ) {
292 std::size_t n = 0; 282 std::size_t n = 0;
293 for ( auto const & elm : *mat->GetElementV 283 for ( auto const & elm : *mat->GetElementVector() ) {
294 std::size_t niso = elm->GetNumberOfIsoto 284 std::size_t niso = elm->GetNumberOfIsotopes();
295 n += niso; 285 n += niso;
296 if(niso > imax) { imax = niso; } 286 if(niso > imax) { imax = niso; }
297 } 287 }
298 if (n > nmax) { nmax = n; } 288 if (n > nmax) { nmax = n; }
299 } 289 }
300 fTemp.resize(imax, 0.0); 290 fTemp.resize(imax, 0.0);
301 fZA.clear(); 291 fZA.clear();
302 fZA.reserve(nmax); 292 fZA.reserve(nmax);
303 fIsoXS.resize(nmax, 0.0); 293 fIsoXS.resize(nmax, 0.0);
304 } 294 }
305 295
306 void G4CrossSectionHP::DumpPhysicsTable(const 296 void G4CrossSectionHP::DumpPhysicsTable(const G4ParticleDefinition&)
307 { 297 {
308 if (fManagerHP->GetVerboseLevel() == 0 || fP 298 if (fManagerHP->GetVerboseLevel() == 0 || fPrinted)
309 return; 299 return;
310 300
311 // 301 //
312 // Dump element based cross section 302 // Dump element based cross section
313 // range 10e-5 eV to 20 MeV 303 // range 10e-5 eV to 20 MeV
314 // 10 point per decade 304 // 10 point per decade
315 // in barn 305 // in barn
316 // 306 //
317 fPrinted = true; 307 fPrinted = true;
318 G4cout << G4endl; 308 G4cout << G4endl;
319 G4cout << "HP Cross Section " << fDataName < 309 G4cout << "HP Cross Section " << fDataName << " for "
320 << fParticle->GetParticleName() << G4endl; 310 << fParticle->GetParticleName() << G4endl;
321 G4cout << "(Pointwise cross-section at 0 Kel 311 G4cout << "(Pointwise cross-section at 0 Kelvin.)" << G4endl;
322 G4cout << G4endl; 312 G4cout << G4endl;
323 G4cout << "Name of Element" << G4endl; 313 G4cout << "Name of Element" << G4endl;
324 G4cout << "Energy[eV] XS[barn]" << G4endl; 314 G4cout << "Energy[eV] XS[barn]" << G4endl;
325 G4cout << G4endl; 315 G4cout << G4endl;
326 316
327 const G4ElementTable* table = G4Element::Get 317 const G4ElementTable* table = G4Element::GetElementTable();
328 for ( auto const & elm : *table ) { 318 for ( auto const & elm : *table ) {
329 G4int Z = elm->GetZasInt(); 319 G4int Z = elm->GetZasInt();
330 if (Z >= minZ && Z <= maxZ && nullptr != f << 320 if (Z >= minZ && Z <= maxZ &&
>> 321 nullptr != fData->GetElementData(Z - minZ)) {
331 G4cout << "----------------------------- 322 G4cout << "---------------------------------------------------" << G4endl;
332 G4cout << elm->GetName() << G4endl; 323 G4cout << elm->GetName() << G4endl;
333 std::size_t n = fData->GetNumberOfCompon 324 std::size_t n = fData->GetNumberOfComponents(Z);
334 for (size_t i=0; i<n; ++i) { 325 for (size_t i=0; i<n; ++i) {
335 G4cout << "## Z=" << Z << " A=" << fData-> 326 G4cout << "## Z=" << Z << " A=" << fData->GetComponentID(Z - minZ, i);
336 G4cout << *(fData->GetComponentDataByIndex(Z 327 G4cout << *(fData->GetComponentDataByIndex(Z - minZ, i)) << G4endl;
337 } 328 }
338 } 329 }
339 } 330 }
340 } 331 }
341 332
342 void G4CrossSectionHP::PrepareCache(const G4Ma 333 void G4CrossSectionHP::PrepareCache(const G4Material* mat)
343 { 334 {
344 fCurrentMat = mat; 335 fCurrentMat = mat;
345 fZA.clear(); 336 fZA.clear();
346 for ( auto const & elm : *(mat->GetElementVe 337 for ( auto const & elm : *(mat->GetElementVector()) ) {
347 G4int Z = elm->GetZasInt(); 338 G4int Z = elm->GetZasInt();
348 for ( auto const & iso : *(elm->GetIsotope 339 for ( auto const & iso : *(elm->GetIsotopeVector()) ) {
349 G4int A = iso->GetN(); 340 G4int A = iso->GetN();
350 fZA.emplace_back(Z, A); 341 fZA.emplace_back(Z, A);
351 } 342 }
352 } 343 }
353 fIsoXS.resize(fZA.size(), 0.0); 344 fIsoXS.resize(fZA.size(), 0.0);
354 } 345 }
355 346
>> 347 void G4CrossSectionHP::InitialiseOnFly(const G4int Z)
>> 348 {
>> 349 G4AutoLock l(&theHPXS);
>> 350 Initialise(Z);
>> 351 l.unlock();
>> 352 }
>> 353
356 void G4CrossSectionHP::Initialise(const G4int 354 void G4CrossSectionHP::Initialise(const G4int Z)
357 { 355 {
358 if (fManagerHP->GetVerboseLevel() > 1) { 356 if (fManagerHP->GetVerboseLevel() > 1) {
359 G4cout << " G4CrossSectionHP::Initialise: 357 G4cout << " G4CrossSectionHP::Initialise: Z=" << Z
360 << " for " << fDataName 358 << " for " << fDataName
361 << " minZ=" << minZ << " maxZ=" << maxZ < 359 << " minZ=" << minZ << " maxZ=" << maxZ << G4endl;
362 } 360 }
363 if (Z < minZ || Z > maxZ || nullptr != fData 361 if (Z < minZ || Z > maxZ || nullptr != fData->GetElementData(Z - minZ)) {
364 return; 362 return;
365 } 363 }
366 G4AutoLock l(&theHPXS); <<
367 if (nullptr != fData->GetElementData(Z - min <<
368 l.unlock(); <<
369 return; <<
370 } <<
371 364
372 // add empty vector to avoid double initiali 365 // add empty vector to avoid double initialisation
373 fData->InitialiseForElement(Z - minZ, new G4 366 fData->InitialiseForElement(Z - minZ, new G4PhysicsVector());
374 367
375 G4String tnam = "temp"; 368 G4String tnam = "temp";
376 G4bool noComp = true; 369 G4bool noComp = true;
377 for (G4int A=amin[Z]; A<=amax[Z]; ++A) { 370 for (G4int A=amin[Z]; A<=amax[Z]; ++A) {
378 std::ostringstream ost; 371 std::ostringstream ost;
379 ost << fDataDirectory; << 372 ost << fDataDirectory << Z << "_";
380 // first check special cases << 373 if (6 == Z && 12 == A) {
381 if (6 == Z && 12 == A && fParticle == fNeu << 374 ost << "nat_";
382 ost << Z << "_nat_" << elementName[Z]; <<
383 } else if (18 == Z && 40 != A) { <<
384 continue; <<
385 } else if (27 == Z && 62 == A) { 375 } else if (27 == Z && 62 == A) {
386 ost << Z << "_62m1_" << elementName[Z]; << 376 ost << "62m1_";
387 } else if (47 == Z && 106 == A) { 377 } else if (47 == Z && 106 == A) {
388 ost << Z << "_106m1_" << elementName[Z]; << 378 ost << "106m1_";
389 } else if (48 == Z && 115 == A) { 379 } else if (48 == Z && 115 == A) {
390 ost << Z << "_115m1_" << elementName[Z]; << 380 ost << "115m1_";
391 } else if (52 == Z && 127 == A) { 381 } else if (52 == Z && 127 == A) {
392 ost << Z << "_127m1_" << elementName[Z]; << 382 ost << "127m1_";
393 } else if (52 == Z && 129 == A) { 383 } else if (52 == Z && 129 == A) {
394 ost << Z << "_129m1_" << elementName[Z]; << 384 ost << "129m1_";
395 } else if (52 == Z && 131 == A) { 385 } else if (52 == Z && 131 == A) {
396 ost << Z << "_131m1_" << elementName[Z]; << 386 ost << "131m1_";
397 } else if (61 == Z && 145 == A) { <<
398 ost << Z << "_147_" << elementName[Z]; <<
399 } else if (67 == Z && 166 == A) { 387 } else if (67 == Z && 166 == A) {
400 ost << Z << "_166m1_" << elementName[Z]; << 388 ost << "166m1_";
401 } else if (73 == Z && 180 == A) { 389 } else if (73 == Z && 180 == A) {
402 ost << Z << "_180m1_" << elementName[Z]; << 390 ost << "180m1_";
403 } else if ((Z == 85 && A == 210) || (Z == <<
404 ost << "84_209_" << elementName[84]; <<
405 } else { 391 } else {
406 // the main file name << 392 ost << A << "_";
407 ost << Z << "_" << A << "_" << elementNa <<
408 } 393 }
>> 394 ost << elementName[Z];
>> 395 std::ifstream filein(ost.str().c_str());
>> 396 //G4cout << "File: " << ost.str() << " " << G4endl;
409 std::istringstream theXSData(tnam, std::io 397 std::istringstream theXSData(tnam, std::ios::in);
410 fManagerHP->GetDataStream(ost.str().c_str( 398 fManagerHP->GetDataStream(ost.str().c_str(), theXSData);
411 if (theXSData) { 399 if (theXSData) {
412 G4int i1, i2, n; 400 G4int i1, i2, n;
413 theXSData >> i1 >> i2 >> n; 401 theXSData >> i1 >> i2 >> n;
414 if (fManagerHP->GetVerboseLevel() > 1) { 402 if (fManagerHP->GetVerboseLevel() > 1) {
415 G4cout << "## G4CrossSectionHP::Initialise f 403 G4cout << "## G4CrossSectionHP::Initialise for Z=" << Z
416 << " A=" << A << " Npoints=" << n << 404 << " A=" << A << " Npoints=" << n << G4endl;
417 } 405 }
418 G4double x, y; 406 G4double x, y;
419 G4PhysicsFreeVector* v = new G4PhysicsFr 407 G4PhysicsFreeVector* v = new G4PhysicsFreeVector(n);
420 for (G4int i=0; i<n; ++i) { 408 for (G4int i=0; i<n; ++i) {
421 theXSData >> x >> y; 409 theXSData >> x >> y;
422 x *= CLHEP::eV; 410 x *= CLHEP::eV;
423 y *= CLHEP::barn; << 411 y *= CLHEP::barn;
424 //G4cout << " e=" << x << " xs=" << y << G 412 //G4cout << " e=" << x << " xs=" << y << G4endl;
425 v->PutValues((std::size_t)i, x, y); 413 v->PutValues((std::size_t)i, x, y);
426 } 414 }
427 v->EnableLogBinSearch(binSearch); << 415 v->EnableLogBinSearch(2);
428 if (noComp) { 416 if (noComp) {
429 G4int nmax = amax[Z] - A + 1; 417 G4int nmax = amax[Z] - A + 1;
430 fData->InitialiseForComponent(Z - minZ, nmax 418 fData->InitialiseForComponent(Z - minZ, nmax);
431 noComp = false; 419 noComp = false;
432 } 420 }
433 fData->AddComponent(Z - minZ, A, v); 421 fData->AddComponent(Z - minZ, A, v);
434 } 422 }
435 } 423 }
436 if (noComp) { fData->InitialiseForComponent( 424 if (noComp) { fData->InitialiseForComponent(Z - minZ, 0); }
437 l.unlock(); <<
438 } 425 }
439 426
440 427