| 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 // ------------------------------------------- 26 // -------------------------------------------------------------------
27 // 27 //
28 // GEANT4 Class file 28 // GEANT4 Class file
29 // 29 //
30 // 30 //
31 // File name: G4ParticleInelasticXS 31 // File name: G4ParticleInelasticXS
32 // 32 //
33 // Author Ivantchenko, Geant4, 3-Aug-09 33 // Author Ivantchenko, Geant4, 3-Aug-09
34 // 34 //
35 // Modifications: 35 // Modifications:
36 // 36 //
37 37
38 #include "G4ParticleInelasticXS.hh" 38 #include "G4ParticleInelasticXS.hh"
39 #include "G4Neutron.hh" 39 #include "G4Neutron.hh"
40 #include "G4DynamicParticle.hh" 40 #include "G4DynamicParticle.hh"
41 #include "G4ElementTable.hh" << 41 #include "G4ProductionCutsTable.hh"
42 #include "G4Material.hh" 42 #include "G4Material.hh"
43 #include "G4Element.hh" 43 #include "G4Element.hh"
44 #include "G4PhysicsLogVector.hh" 44 #include "G4PhysicsLogVector.hh"
45 #include "G4CrossSectionDataSetRegistry.hh" 45 #include "G4CrossSectionDataSetRegistry.hh"
46 #include "G4ComponentGGHadronNucleusXsc.hh" 46 #include "G4ComponentGGHadronNucleusXsc.hh"
47 #include "G4ComponentGGNuclNuclXsc.hh" 47 #include "G4ComponentGGNuclNuclXsc.hh"
48 #include "G4HadronicParameters.hh" <<
49 #include "G4Proton.hh" 48 #include "G4Proton.hh"
50 #include "Randomize.hh" 49 #include "Randomize.hh"
51 #include "G4SystemOfUnits.hh" 50 #include "G4SystemOfUnits.hh"
52 #include "G4IsotopeList.hh" 51 #include "G4IsotopeList.hh"
53 #include "G4HadronicParameters.hh" 52 #include "G4HadronicParameters.hh"
54 #include "G4AutoLock.hh" <<
55 53
56 #include <fstream> 54 #include <fstream>
57 #include <sstream> 55 #include <sstream>
58 56
59 G4ElementData* G4ParticleInelasticXS::data[] = << 57 G4ElementData* G4ParticleInelasticXS::data[] = {nullptr};
60 G4double G4ParticleInelasticXS::coeff[MAXZINEL 58 G4double G4ParticleInelasticXS::coeff[MAXZINELP][5] = {{1.0}, {1.0}, {1.0}, {1.0}, {1.0}};
61 G4String G4ParticleInelasticXS::gDataDirectory << 59 G4String G4ParticleInelasticXS::gDataDirectory[] = {""};
62 60
63 namespace << 61 #ifdef G4MULTITHREADED
64 { << 62 G4Mutex G4ParticleInelasticXS::particleInelasticXSMutex = G4MUTEX_INITIALIZER;
65 G4Mutex pInelasticXSMutex = G4MUTEX_INITIALI << 63 #endif
66 G4String pname[5] = {"proton", "deuteron", " <<
67 } <<
68 64
69 G4ParticleInelasticXS::G4ParticleInelasticXS(c 65 G4ParticleInelasticXS::G4ParticleInelasticXS(const G4ParticleDefinition* part)
70 : G4VCrossSectionDataSet("G4ParticleInelasti 66 : G4VCrossSectionDataSet("G4ParticleInelasticXS"),
>> 67 highEnergyXsection(nullptr),
71 particle(part), 68 particle(part),
72 elimit(20*CLHEP::MeV) << 69 index(0),
>> 70 isMaster(false)
73 { 71 {
74 if (nullptr == part) { << 72 if(!part) {
75 G4Exception("G4ParticleInelasticXS::G4Part 73 G4Exception("G4ParticleInelasticXS::G4ParticleInelasticXS(..)","had015",
76 FatalException, "NO particle definition in 74 FatalException, "NO particle definition in constructor");
77 } else { 75 } else {
78 verboseLevel = 0; 76 verboseLevel = 0;
79 const G4String& particleName = particle->G 77 const G4String& particleName = particle->GetParticleName();
80 if(verboseLevel > 1) { 78 if(verboseLevel > 1) {
81 G4cout << "G4ParticleInelasticXS::G4Part 79 G4cout << "G4ParticleInelasticXS::G4ParticleInelasticXS for "
82 << particleName << " on atoms with Z < 80 << particleName << " on atoms with Z < " << MAXZINELP << G4endl;
83 } 81 }
84 auto xsr = G4CrossSectionDataSetRegistry:: << 82 if(particleName == "proton") {
85 if (particleName == "proton") { << 83 highEnergyXsection = G4CrossSectionDataSetRegistry::Instance()->GetComponentCrossSection("Glauber-Gribov");
86 highEnergyXsection = xsr->GetComponentCr <<
87 if(highEnergyXsection == nullptr) { 84 if(highEnergyXsection == nullptr) {
88 highEnergyXsection = new G4ComponentGGHadron 85 highEnergyXsection = new G4ComponentGGHadronNucleusXsc();
89 } 86 }
90 } else { 87 } else {
91 highEnergyXsection = << 88 highEnergyXsection = G4CrossSectionDataSetRegistry::Instance()->GetComponentCrossSection("Glauber-Gribov Nucl-nucl");
92 xsr->GetComponentCrossSection("Glauber <<
93 if(highEnergyXsection == nullptr) { 89 if(highEnergyXsection == nullptr) {
94 highEnergyXsection = new G4ComponentGGNuclNu 90 highEnergyXsection = new G4ComponentGGNuclNuclXsc();
95 } 91 }
96 for (index=1; index<5; ++index) { << 92 if(particleName == "deuteron") index = 1;
97 if (particleName == pname[index]) { br << 93 else if(particleName == "triton") index = 2;
98 } << 94 else if(particleName == "He3") index = 3;
99 index = std::min(index, 4); << 95 else if(particleName == "alpha") index = 4;
100 if (1 < index) { SetMaxKinEnergy(25.6*CL << 96 else {
>> 97 G4ExceptionDescription ed;
>> 98 ed << particleName << " is a wrong particle type";
>> 99 G4Exception("G4ParticleInelasticXS::BuildPhysicsTable(..)","had012",
>> 100 FatalException, ed, "");
>> 101 }
101 } 102 }
102 } 103 }
103 SetForAllAtomsAndEnergies(true); 104 SetForAllAtomsAndEnergies(true);
104 if (gDataDirectory.empty()) { << 105 temp.resize(13,0.0);
105 gDataDirectory = G4HadronicParameters::Ins << 106 }
106 } << 107
107 G4String ss = pname[index] + "ParticleXS"; << 108 G4ParticleInelasticXS::~G4ParticleInelasticXS()
108 SetName(ss); << 109 {
109 if (data[index] == nullptr) { << 110 if(isMaster) {
110 data[index] = new G4ElementData(MAXZINELP) << 111 delete data[index];
111 data[index]->SetName(pname[index] + "PartI << 112 data[index] = nullptr;
112 } 113 }
113 } 114 }
114 115
115 void G4ParticleInelasticXS::CrossSectionDescri 116 void G4ParticleInelasticXS::CrossSectionDescription(std::ostream& outFile) const
116 { 117 {
117 outFile << "G4ParticleInelasticXS calculates 118 outFile << "G4ParticleInelasticXS calculates n, p, d, t, he3, he4 inelastic\n"
118 << "cross section on nuclei using da 119 << "cross section on nuclei using data from the high precision\n"
119 << "neutron database. These data ar 120 << "neutron database. These data are simplified and smoothed over\n"
120 << "the resonance region in order to 121 << "the resonance region in order to reduce CPU time.\n"
121 << "For high energy Glauber-Gribov c << 122 << "For high energy Glauber-Gribov cross section model is used\n";
122 } 123 }
123 124
124 G4bool 125 G4bool
125 G4ParticleInelasticXS::IsElementApplicable(con 126 G4ParticleInelasticXS::IsElementApplicable(const G4DynamicParticle*,
126 G4int, const G4Material*) 127 G4int, const G4Material*)
127 { 128 {
128 return true; 129 return true;
129 } 130 }
130 131
131 G4bool 132 G4bool
132 G4ParticleInelasticXS::IsIsoApplicable(const G 133 G4ParticleInelasticXS::IsIsoApplicable(const G4DynamicParticle*,
133 G4int, G4int, 134 G4int, G4int,
134 const G4Element*, const G4Mater 135 const G4Element*, const G4Material*)
135 { 136 {
136 return true; 137 return true;
137 } 138 }
138 139
139 G4double << 140 G4double G4ParticleInelasticXS::GetElementCrossSection(
140 G4ParticleInelasticXS::GetElementCrossSection( << 141 const G4DynamicParticle* aParticle,
141 << 142 G4int ZZ, const G4Material*)
142 { 143 {
143 return ElementCrossSection(aParticle->GetKin << 144 G4double xs = 0.0;
144 aParticle->GetLog << 145 G4double ekin = aParticle->GetKineticEnergy();
145 } <<
146 <<
147 G4double <<
148 G4ParticleInelasticXS::ComputeCrossSectionPerE <<
149 <<
150 <<
151 <<
152 { <<
153 return ElementCrossSection(ekin, loge, elm-> <<
154 } <<
155 146
156 G4double G4ParticleInelasticXS::ElementCrossSe <<
157 { <<
158 G4int Z = (ZZ >= MAXZINELP) ? MAXZINELP - 1 147 G4int Z = (ZZ >= MAXZINELP) ? MAXZINELP - 1 : ZZ;
159 148
160 // element data is always valid pointer by c <<
161 auto pv = GetPhysicsVector(Z); 149 auto pv = GetPhysicsVector(Z);
>> 150 if(nullptr == pv) { return xs; }
>> 151 // G4cout << "G4ParticleInelasticXS::GetCrossSection e= " << ekin
>> 152 // << " Z= " << Z << G4endl;
162 153
163 // set to null x-section below lowest energy << 154 if(ekin <= pv->GetMaxEnergy()) {
164 G4double xs = 0.0; << 155 xs = pv->LogVectorValue(ekin, aParticle->GetLogKineticEnergy());
165 if (ekin > pv->Energy(0)) { << 156 } else {
166 xs = (ekin <= pv->GetMaxEnergy()) ? pv->Lo << 157 xs = coeff[Z][index]*highEnergyXsection->GetInelasticElementCrossSection(particle,
167 : coeff[Z][index]*highEnergyXsection->GetI <<
168 ekin, Z, aeff[Z]); 158 ekin, Z, aeff[Z]);
169 } 159 }
170 160
171 #ifdef G4VERBOSE 161 #ifdef G4VERBOSE
172 if(verboseLevel > 1) { 162 if(verboseLevel > 1) {
173 G4cout << "ElmXS: Z= " << Z << " Ekin(MeV 163 G4cout << "ElmXS: Z= " << Z << " Ekin(MeV)= " << ekin/CLHEP::MeV
174 << " xs(bn)= " << xs/CLHEP::barn << " el 164 << " xs(bn)= " << xs/CLHEP::barn << " element data for "
175 << particle->GetParticleName() 165 << particle->GetParticleName()
176 << " idx= " << index << G4endl; 166 << " idx= " << index << G4endl;
177 } 167 }
178 #endif 168 #endif
179 return xs; 169 return xs;
180 } 170 }
181 171
182 G4double << 172 G4double G4ParticleInelasticXS::GetIsoCrossSection(
183 G4ParticleInelasticXS::ComputeIsoCrossSection( << 173 const G4DynamicParticle* aParticle,
184 << 174 G4int Z, G4int A,
185 << 175 const G4Isotope*, const G4Element*, const G4Material*)
186 <<
187 { <<
188 return IsoCrossSection(ekin, loge, Z, A); <<
189 } <<
190 <<
191 G4double <<
192 G4ParticleInelasticXS::GetIsoCrossSection(cons <<
193 G4in <<
194 cons <<
195 { 176 {
196 return IsoCrossSection(aParticle->GetKinetic 177 return IsoCrossSection(aParticle->GetKineticEnergy(),
197 aParticle->GetLogKine << 178 aParticle->GetLogKineticEnergy(),Z, A);
198 } 179 }
199 180
200 G4double 181 G4double
201 G4ParticleInelasticXS::IsoCrossSection(G4doubl 182 G4ParticleInelasticXS::IsoCrossSection(G4double ekin, G4double logE,
202 G4int Z 183 G4int ZZ, G4int A)
203 { 184 {
204 G4double xs = 0.0; 185 G4double xs = 0.0;
205 G4int Z = (ZZ >= MAXZINELP) ? MAXZINELP - 1 186 G4int Z = (ZZ >= MAXZINELP) ? MAXZINELP - 1 : ZZ;
206 << 187 /*
207 // needed here to gurantee upload data for Z << 188 G4cout << "G4ParticleInelasticXS: IsoCrossSection Z= "
>> 189 << Z << " A= " << A
>> 190 << " Amin= " << amin[Z] << " Amax= " << amax[Z]
>> 191 << " E(MeV)= " << ekin << G4endl;
>> 192 */
208 auto pv = GetPhysicsVector(Z); 193 auto pv = GetPhysicsVector(Z);
>> 194 if(pv == nullptr) { return xs; }
209 195
210 // compute isotope cross section if applicab 196 // compute isotope cross section if applicable
211 if (ekin <= elimit && data[index]->GetNumber << 197 const G4double emax = pv->GetMaxEnergy();
212 auto pviso = data[index]->GetComponentData << 198 if(ekin <= emax && amin[Z]>0 && A >= amin[Z] && A <= amax[Z]) {
213 if (pviso != nullptr && ekin > pviso->Ener << 199 auto pviso = data[index]->GetComponentDataByIndex(Z, A - amin[Z]);
>> 200 if(pviso != nullptr) {
214 xs = pviso->LogVectorValue(ekin, logE); 201 xs = pviso->LogVectorValue(ekin, logE);
215 #ifdef G4VERBOSE 202 #ifdef G4VERBOSE
216 if(verboseLevel > 1) { 203 if(verboseLevel > 1) {
217 G4cout << "G4ParticleInelasticXS::IsoXS: for 204 G4cout << "G4ParticleInelasticXS::IsoXS: for "
218 << particle->GetParticleName() 205 << particle->GetParticleName() << " Ekin(MeV)= "
219 << ekin/CLHEP::MeV << " xs(b)= 206 << ekin/CLHEP::MeV << " xs(b)= " << xs/CLHEP::barn
220 << " Z= " << Z << " A= " << A 207 << " Z= " << Z << " A= " << A
221 << " idx= " << index << G4endl; 208 << " idx= " << index << G4endl;
222 } 209 }
223 #endif 210 #endif
224 return xs; 211 return xs;
225 } 212 }
226 } 213 }
227 // use element x-section 214 // use element x-section
228 if (ekin > pv->Energy(0)) { << 215 if(ekin <= emax) {
229 xs = (ekin <= pv->GetMaxEnergy()) ? pv->Lo << 216 xs = pv->LogVectorValue(ekin, logE);
230 coeff[Z][index] * << 217 } else {
231 highEnergyXsection->GetInelasticElementC << 218 xs = coeff[Z][index] *
232 * A/aeff[Z]; << 219 highEnergyXsection->GetInelasticElementCrossSection(particle,
>> 220 ekin, Z, aeff[Z]);
233 } 221 }
>> 222 xs *= A/aeff[Z];
234 #ifdef G4VERBOSE 223 #ifdef G4VERBOSE
235 if(verboseLevel > 1) { 224 if(verboseLevel > 1) {
236 G4cout << "IsoXS for " << particle->GetPa 225 G4cout << "IsoXS for " << particle->GetParticleName()
237 << " Target Z= " << Z << " A= " << A 226 << " Target Z= " << Z << " A= " << A
238 << " Ekin(MeV)= " << ekin/CLHEP::MeV 227 << " Ekin(MeV)= " << ekin/CLHEP::MeV
239 << " xs(bn)= " << xs/CLHEP::barn 228 << " xs(bn)= " << xs/CLHEP::barn
240 << " idx= " << index << G4endl; 229 << " idx= " << index << G4endl;
241 } 230 }
242 #endif 231 #endif
243 return xs; 232 return xs;
244 } 233 }
245 234
246 const G4Isotope* G4ParticleInelasticXS::Select 235 const G4Isotope* G4ParticleInelasticXS::SelectIsotope(
247 const G4Element* anElement, G4double kinE 236 const G4Element* anElement, G4double kinEnergy, G4double logE)
248 { 237 {
249 G4int nIso = (G4int)anElement->GetNumberOfIs << 238 size_t nIso = anElement->GetNumberOfIsotopes();
250 const G4Isotope* iso = anElement->GetIsotope 239 const G4Isotope* iso = anElement->GetIsotope(0);
251 240
252 if (1 == nIso) { return iso; } << 241 //G4cout << "SelectIsotope NIso= " << nIso << G4endl;
>> 242 if(1 == nIso) { return iso; }
253 243
254 // more than 1 isotope 244 // more than 1 isotope
255 G4int Z = anElement->GetZasInt(); 245 G4int Z = anElement->GetZasInt();
256 << 246 //G4cout << "SelectIsotope Z= " << Z << G4endl;
257 // initialisation for given Z <<
258 GetPhysicsVector(Z); <<
259 247
260 const G4double* abundVector = anElement->Get 248 const G4double* abundVector = anElement->GetRelativeAbundanceVector();
261 G4double q = G4UniformRand(); 249 G4double q = G4UniformRand();
262 G4double sum = 0.0; 250 G4double sum = 0.0;
263 G4int j; << 251 size_t j;
264 252
265 // isotope wise cross section not available 253 // isotope wise cross section not available
266 if (Z >= MAXZINELP || 0 == data[index]->GetN << 254 if(0 == amin[Z] || Z >= MAXZINELP) {
267 for (j=0; j<nIso; ++j) { 255 for (j=0; j<nIso; ++j) {
268 sum += abundVector[j]; 256 sum += abundVector[j];
269 if(q <= sum) { 257 if(q <= sum) {
270 iso = anElement->GetIsotope(j); 258 iso = anElement->GetIsotope(j);
271 break; 259 break;
272 } 260 }
273 } 261 }
274 return iso; 262 return iso;
275 } 263 }
276 264
277 G4int nn = (G4int)temp.size(); << 265 size_t nn = temp.size();
278 if (nn < nIso) { temp.resize(nIso, 0.); } << 266 if(nn < nIso) { temp.resize(nIso, 0.); }
279 267
280 for (j=0; j<nIso; ++j) { 268 for (j=0; j<nIso; ++j) {
>> 269 //G4cout << j << "-th isotope " << (*isoVector)[j]->GetN()
>> 270 // << " abund= " << abundVector[j] << G4endl;
281 sum += abundVector[j]*IsoCrossSection(kinE 271 sum += abundVector[j]*IsoCrossSection(kinEnergy, logE, Z,
282 anElement->GetIsotope(j)->GetN()); 272 anElement->GetIsotope(j)->GetN());
283 temp[j] = sum; 273 temp[j] = sum;
284 } 274 }
285 sum *= q; 275 sum *= q;
286 for (j=0; j<nIso; ++j) { 276 for (j=0; j<nIso; ++j) {
287 if (temp[j] >= sum) { << 277 if(temp[j] >= sum) {
288 iso = anElement->GetIsotope(j); 278 iso = anElement->GetIsotope(j);
289 break; 279 break;
290 } 280 }
291 } 281 }
292 return iso; 282 return iso;
293 } 283 }
294 284
295 void 285 void
296 G4ParticleInelasticXS::BuildPhysicsTable(const 286 G4ParticleInelasticXS::BuildPhysicsTable(const G4ParticleDefinition& p)
297 { 287 {
298 if (verboseLevel > 0){ << 288 if(verboseLevel > 0){
299 G4cout << "G4ParticleInelasticXS::BuildPhy 289 G4cout << "G4ParticleInelasticXS::BuildPhysicsTable for "
300 << p.GetParticleName() << G4endl; 290 << p.GetParticleName() << G4endl;
301 } 291 }
302 if (&p != particle) { << 292 if(&p != particle) {
303 G4ExceptionDescription ed; 293 G4ExceptionDescription ed;
304 ed << p.GetParticleName() << " is a wrong 294 ed << p.GetParticleName() << " is a wrong particle type -"
305 << particle->GetParticleName() << " is 295 << particle->GetParticleName() << " is expected";
306 G4Exception("G4ParticleInelasticXS::BuildP 296 G4Exception("G4ParticleInelasticXS::BuildPhysicsTable(..)","had012",
307 FatalException, ed, ""); 297 FatalException, ed, "");
308 return; 298 return;
309 } 299 }
310 300
>> 301 G4int fact = (p.GetParticleName() == "proton") ? 1 : 256;
>> 302 SetMaxKinEnergy(G4HadronicParameters::Instance()->GetMaxEnergy() * fact);
>> 303
>> 304 if(data[index] == nullptr) {
>> 305 #ifdef G4MULTITHREADED
>> 306 G4MUTEXLOCK(&particleInelasticXSMutex);
>> 307 if(data[index] == nullptr) {
>> 308 #endif
>> 309 isMaster = true;
>> 310 data[index] = new G4ElementData();
>> 311 data[index]->SetName(particle->GetParticleName() + "Inelastic");
>> 312 FindDirectoryPath();
>> 313 #ifdef G4MULTITHREADED
>> 314 }
>> 315 G4MUTEXUNLOCK(&particleInelasticXSMutex);
>> 316 #endif
>> 317 }
>> 318
311 // it is possible re-initialisation for the 319 // it is possible re-initialisation for the new run
312 const G4ElementTable* table = G4Element::Get << 320 if(isMaster) {
313 321
314 // prepare isotope selection << 322 // Access to elements
315 std::size_t nIso = temp.size(); << 323 auto theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
>> 324 size_t numOfCouples = theCoupleTable->GetTableSize();
>> 325 for(size_t j=0; j<numOfCouples; ++j) {
>> 326 auto mat = theCoupleTable->GetMaterialCutsCouple(j)->GetMaterial();
>> 327 auto elmVec = mat->GetElementVector();
>> 328 size_t numOfElem = mat->GetNumberOfElements();
>> 329 for (size_t ie = 0; ie < numOfElem; ++ie) {
>> 330 G4int Z = std::max(1,std::min(((*elmVec)[ie])->GetZasInt(), MAXZINELP-1));
>> 331 if(nullptr == data[index]->GetElementData(Z)) { Initialise(Z); }
>> 332 }
>> 333 }
>> 334 }
>> 335 }
316 336
317 // Access to elements << 337 const G4String& G4ParticleInelasticXS::FindDirectoryPath()
318 for ( auto const & elm : *table ) { << 338 {
319 std::size_t n = elm->GetNumberOfIsotopes() << 339 // check environment variable
320 if (n > nIso) { nIso = n; } << 340 // build the complete string identifying the file with the data set
321 G4int Z = std::min( elm->GetZasInt(), MAXZ << 341 if(gDataDirectory[index].empty()) {
322 if ( nullptr == (data[index])->GetElementD << 342 char* path = std::getenv("G4PARTICLEXSDATA");
323 Initialise(Z); << 343 if (path) {
>> 344 std::ostringstream ost;
>> 345 ost << path << "/" << particle->GetParticleName() << "/inel";
>> 346 gDataDirectory[index] = ost.str();
>> 347 } else {
>> 348 G4Exception("G4NeutronInelasticXS::Initialise(..)","had013",
>> 349 FatalException,
>> 350 "Environment variable G4PARTICLEXSDATA is not defined");
324 } 351 }
325 } 352 }
326 temp.resize(nIso, 0.0); << 353 return gDataDirectory[index];
>> 354 }
>> 355
>> 356 void G4ParticleInelasticXS::InitialiseOnFly(G4int Z)
>> 357 {
>> 358 #ifdef G4MULTITHREADED
>> 359 G4MUTEXLOCK(&particleInelasticXSMutex);
>> 360 if(nullptr == data[index]->GetElementData(Z)) {
>> 361 #endif
>> 362 Initialise(Z);
>> 363 #ifdef G4MULTITHREADED
>> 364 }
>> 365 G4MUTEXUNLOCK(&particleInelasticXSMutex);
>> 366 #endif
327 } 367 }
328 368
329 void G4ParticleInelasticXS::Initialise(G4int Z 369 void G4ParticleInelasticXS::Initialise(G4int Z)
330 { 370 {
331 if ( nullptr != (data[index])->GetElementDat << 371 if(nullptr != data[index]->GetElementData(Z)) { return; }
332 372
333 G4AutoLock l(&pInelasticXSMutex); << 373 // upload element data
334 if ( nullptr == (data[index])->GetElementDat << 374 std::ostringstream ost;
335 // upload element data << 375 ost << FindDirectoryPath() << Z ;
336 std::ostringstream ost; << 376 G4PhysicsVector* v = RetrieveVector(ost, true);
337 ost << gDataDirectory << "/" << pname[inde << 377 data[index]->InitialiseForElement(Z, v);
338 G4PhysicsVector* v = RetrieveVector(ost, t << 378 /*
339 data[index]->InitialiseForElement(Z, v); << 379 G4cout << "G4ParticleInelasticXS::Initialise for Z= " << Z
340 << 380 << " idx= " << index
341 // upload isotope data << 381 << " Amin= " << amin[Z]
342 G4bool noComp = true; << 382 << " Amax= " << amax[Z]
343 if (amin[Z] < amax[Z]) { << 383 << " " << FindDirectoryPath() << G4endl;
344 << 384 */
345 for (G4int A=amin[Z]; A<=amax[Z]; ++A) { << 385 // upload isotope data
346 std::ostringstream ost1; << 386 if(amin[Z] > 0) {
347 ost1 << gDataDirectory << "/" << pname[index << 387 size_t nmax = (size_t)(amax[Z]-amin[Z]+1);
348 G4PhysicsVector* v1 = RetrieveVector(ost1, f << 388 data[index]->InitialiseForComponent(Z, nmax);
349 if (nullptr != v1) { << 389
350 if (noComp) { << 390 for(G4int A=amin[Z]; A<=amax[Z]; ++A) {
351 G4int nmax = amax[Z] - A + 1; << 391 std::ostringstream ost1;
352 data[index]->InitialiseForComponent(Z, n << 392 ost1 << FindDirectoryPath() << Z << "_" << A;
353 noComp = false; << 393 G4PhysicsVector* v1 = RetrieveVector(ost1, false);
354 } << 394 data[index]->AddComponent(Z, A, v1);
355 data[index]->AddComponent(Z, A, v1); << 395 /*
356 } << 396 G4cout << " Isotope x-section Z= " << Z << " A= " << A
357 } << 397 << " v1= " << v1 << G4endl;
>> 398 */
358 } 399 }
359 // no components case <<
360 if (noComp) { data[index]->InitialiseForCo <<
361 <<
362 // smooth transition <<
363 G4double sig1 = (*v)[v->GetVectorLength() <<
364 G4double ehigh = v->GetMaxEnergy(); <<
365 G4double sig2 = highEnergyXsection->GetIne <<
366 particle, ehigh, Z, aeff[Z]); <<
367 coeff[Z][index] = (sig2 > 0.) ? sig1/sig2 <<
368 } 400 }
369 l.unlock(); << 401 // smooth transition
>> 402 G4double sig1 = (*v)[v->GetVectorLength()-1];
>> 403 G4double ehigh = v->GetMaxEnergy();
>> 404 G4double sig2 = highEnergyXsection->GetInelasticElementCrossSection(
>> 405 particle, ehigh, Z, aeff[Z]);
>> 406 coeff[Z][index] = (sig2 > 0.) ? sig1/sig2 : 1.0;
>> 407 /*
>> 408 G4cout << "G4ParticleInelasticXS: index= " << index
>> 409 << " Z= " << Z << " coeff= " << coeff[Z][index] << G4endl;
>> 410 */
370 } 411 }
371 412
372 G4PhysicsVector* 413 G4PhysicsVector*
373 G4ParticleInelasticXS::RetrieveVector(std::ost 414 G4ParticleInelasticXS::RetrieveVector(std::ostringstream& ost, G4bool warn)
374 { 415 {
375 G4PhysicsLogVector* v = nullptr; 416 G4PhysicsLogVector* v = nullptr;
376 std::ifstream filein(ost.str().c_str()); 417 std::ifstream filein(ost.str().c_str());
377 if (!filein.is_open()) { << 418 if (!(filein)) {
378 if (warn) { << 419 if(warn) {
379 G4ExceptionDescription ed; 420 G4ExceptionDescription ed;
380 ed << "Data file <" << ost.str().c_str() 421 ed << "Data file <" << ost.str().c_str()
381 << "> is not opened! index=" << index << 422 << "> is not opened!";
382 << " dir: <" << gDataDirectory << ">. "; <<
383 G4Exception("G4ParticleInelasticXS::Retr 423 G4Exception("G4ParticleInelasticXS::RetrieveVector(..)","had014",
384 FatalException, ed, "Check G4PARTICLEXSD 424 FatalException, ed, "Check G4PARTICLEXSDATA");
385 } 425 }
386 } else { 426 } else {
387 if(verboseLevel > 1) { 427 if(verboseLevel > 1) {
388 G4cout << "File " << ost.str() 428 G4cout << "File " << ost.str()
389 << " is opened by G4ParticleInelasticXS 429 << " is opened by G4ParticleInelasticXS" << G4endl;
390 } 430 }
391 // retrieve data from DB 431 // retrieve data from DB
392 v = new G4PhysicsLogVector(); 432 v = new G4PhysicsLogVector();
393 if(!v->Retrieve(filein, true)) { 433 if(!v->Retrieve(filein, true)) {
394 G4ExceptionDescription ed; 434 G4ExceptionDescription ed;
395 ed << "Data file <" << ost.str().c_str() 435 ed << "Data file <" << ost.str().c_str()
396 << "> is not retrieved!"; 436 << "> is not retrieved!";
397 G4Exception("G4ParticleInelasticXS::Retr 437 G4Exception("G4ParticleInelasticXS::RetrieveVector(..)","had015",
398 FatalException, ed, "Check G4PARTICLEXSD 438 FatalException, ed, "Check G4PARTICLEXSDATA");
399 } 439 }
400 } 440 }
401 return v; 441 return v;
402 } 442 }
403 443