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