| 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: G4GammaNuclearXS 31 // File name: G4GammaNuclearXS
32 // 32 //
33 // Authors V.Ivantchenko, Geant4, 20 October 33 // Authors V.Ivantchenko, Geant4, 20 October 2020
34 // B.Kutsenko, BINP/NSU, 10 August 20 34 // B.Kutsenko, BINP/NSU, 10 August 2021
35 // 35 //
36 // Modifications: 36 // Modifications:
37 // 37 //
38 38
39 #include "G4GammaNuclearXS.hh" 39 #include "G4GammaNuclearXS.hh"
40 #include "G4DynamicParticle.hh" 40 #include "G4DynamicParticle.hh"
41 #include "G4ThreeVector.hh" 41 #include "G4ThreeVector.hh"
42 #include "G4ElementTable.hh" 42 #include "G4ElementTable.hh"
43 #include "G4Material.hh" 43 #include "G4Material.hh"
44 #include "G4Element.hh" 44 #include "G4Element.hh"
45 #include "G4ElementData.hh" <<
46 #include "G4PhysicsVector.hh" <<
47 #include "G4PhysicsLinearVector.hh" 45 #include "G4PhysicsLinearVector.hh"
48 #include "G4PhysicsFreeVector.hh" <<
49 #include "G4CrossSectionDataSetRegistry.hh" 46 #include "G4CrossSectionDataSetRegistry.hh"
50 #include "G4PhotoNuclearCrossSection.hh" 47 #include "G4PhotoNuclearCrossSection.hh"
51 #include "G4HadronicParameters.hh" <<
52 #include "Randomize.hh" 48 #include "Randomize.hh"
53 #include "G4SystemOfUnits.hh" 49 #include "G4SystemOfUnits.hh"
54 #include "G4Gamma.hh" 50 #include "G4Gamma.hh"
55 #include "G4IsotopeList.hh" 51 #include "G4IsotopeList.hh"
56 #include "G4AutoLock.hh" <<
57 52
58 #include <fstream> 53 #include <fstream>
59 #include <sstream> 54 #include <sstream>
60 #include <vector> 55 #include <vector>
61 56
62 G4ElementData* G4GammaNuclearXS::data = nullpt 57 G4ElementData* G4GammaNuclearXS::data = nullptr;
63 58
64 G4double G4GammaNuclearXS::coeff[3][3]; 59 G4double G4GammaNuclearXS::coeff[3][3];
65 G4double G4GammaNuclearXS::xs150[] = {0.0}; << 60 G4double G4GammaNuclearXS::xs150[MAXZGAMMAXS] = {0.0};
66 const G4int G4GammaNuclearXS::freeVectorExcept <<
67 4, 6, 7, 8, 27, 39, 45, 65, 67, 69, 73}; <<
68 G4String G4GammaNuclearXS::gDataDirectory = "" 61 G4String G4GammaNuclearXS::gDataDirectory = "";
69 62
70 namespace << 63 #ifdef G4MULTITHREADED
71 { << 64 G4Mutex G4GammaNuclearXS::gNuclearXSMutex = G4MUTEX_INITIALIZER;
72 // Upper limit of the linear transition betw << 65 #endif
73 const G4double eTransitionBound = 150.*CLHEP <<
74 // A limit energy to correct CHIPS parameter <<
75 const G4double ehigh = 10*CLHEP::GeV; <<
76 } <<
77 66
78 G4GammaNuclearXS::G4GammaNuclearXS() 67 G4GammaNuclearXS::G4GammaNuclearXS()
79 : G4VCrossSectionDataSet(Default_Name()), ga << 68 : G4VCrossSectionDataSet(Default_Name()),
>> 69 gamma(G4Gamma::Gamma())
80 { 70 {
81 verboseLevel = 0; << 71 // verboseLevel = 0;
82 if (verboseLevel > 0) { << 72 if(verboseLevel > 0){
83 G4cout << "G4GammaNuclearXS::G4GammaNuclea << 73 G4cout << "G4GammaNuclearXS::G4GammaNuclearXS Initialise for Z < "
84 << MAXZGAMMAXS << G4endl; << 74 << MAXZGAMMAXS << G4endl;
85 } <<
86 ggXsection = dynamic_cast<G4PhotoNuclearCros <<
87 (G4CrossSectionDataSetRegistry::Instance() <<
88 if (ggXsection == nullptr) { <<
89 ggXsection = new G4PhotoNuclearCrossSectio <<
90 } 75 }
>> 76 ggXsection = G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet("PhotoNuclearXS");
>> 77 if(ggXsection == nullptr) ggXsection = new G4PhotoNuclearCrossSection();
91 SetForAllAtomsAndEnergies(true); 78 SetForAllAtomsAndEnergies(true);
>> 79 }
92 80
93 // full data set is uploaded once << 81 G4GammaNuclearXS::~G4GammaNuclearXS()
94 if (nullptr == data) { << 82 {
95 data = new G4ElementData(MAXZGAMMAXS); << 83 if(isMaster) {
96 data->SetName("gNuclear"); << 84 delete data;
97 for (G4int Z=1; Z<MAXZGAMMAXS; ++Z) { << 85 data = nullptr;
98 Initialise(Z); <<
99 } <<
100 } 86 }
101 } 87 }
102 88
103 void G4GammaNuclearXS::CrossSectionDescription 89 void G4GammaNuclearXS::CrossSectionDescription(std::ostream& outFile) const
104 { 90 {
105 outFile << "G4GammaNuclearXS calculates the 91 outFile << "G4GammaNuclearXS calculates the gamma nuclear\n"
106 << "cross-section for GDR energy reg << 92 << "cross-section for GDR energy region on nuclei using data from the high precision\n"
107 << "data from the high precision\n" <<
108 << "IAEA photonuclear database (2019 93 << "IAEA photonuclear database (2019). Then liniear connection\n"
109 << "implemented with previous CHIPS photon << 94 <<"implemented with previous CHIPS photonuclear model\n";
110 } 95 }
111 96
112 G4bool 97 G4bool
113 G4GammaNuclearXS::IsElementApplicable(const G4 98 G4GammaNuclearXS::IsElementApplicable(const G4DynamicParticle*,
114 G4int, c 99 G4int, const G4Material*)
115 { 100 {
116 return true; 101 return true;
117 } 102 }
118 103
119 G4bool G4GammaNuclearXS::IsIsoApplicable(const 104 G4bool G4GammaNuclearXS::IsIsoApplicable(const G4DynamicParticle*,
120 G4int 105 G4int, G4int,
121 const 106 const G4Element*, const G4Material*)
122 { 107 {
123 return true; 108 return true;
124 } 109 }
125 110
126 G4double 111 G4double
127 G4GammaNuclearXS::GetElementCrossSection(const 112 G4GammaNuclearXS::GetElementCrossSection(const G4DynamicParticle* aParticle,
128 G4int << 113 G4int ZZ, const G4Material* mat)
129 { 114 {
130 return ElementCrossSection(aParticle->GetKin << 115 const G4int Z = (ZZ >= MAXZGAMMAXS) ? MAXZGAMMAXS - 1 : ZZ;
131 } << 116 auto pv = GetPhysicsVector(Z);
132 117
133 G4double << 118 if(pv == nullptr) {
134 G4GammaNuclearXS::ElementCrossSection(const G4 << 119 return ggXsection->GetElementCrossSection(aParticle, Z, mat);
135 { <<
136 // check cache <<
137 const G4int Z = (ZZ < MAXZGAMMAXS) ? ZZ : MA <<
138 if(Z == fZ && ekin == fEkin) { return fXS; } <<
139 fZ = Z; <<
140 fEkin = ekin; <<
141 <<
142 auto pv = data->GetElementData(Z); <<
143 const G4double limCHIPS1 = 25*CLHEP::MeV; <<
144 const G4double limCHIPS2 = 16*CLHEP::MeV; <<
145 if (pv == nullptr || 1 == Z || Z == 40 || Z <<
146 (Z == 24 && ekin >= limCHIPS1) || <<
147 (Z == 39 && ekin >= limCHIPS1) || <<
148 (Z == 50 && ekin >= limCHIPS2) || <<
149 (Z == 64 && ekin >= limCHIPS2) <<
150 ) { <<
151 fXS = ggXsection->ComputeElementXSection(e <<
152 return fXS; <<
153 } 120 }
154 const G4double emax = pv->GetMaxEnergy(); 121 const G4double emax = pv->GetMaxEnergy();
155 << 122 const G4double ekin = aParticle->GetKineticEnergy();
156 // low energy based on data << 123 G4double xs = 0.0;
157 if(ekin <= emax) { 124 if(ekin <= emax) {
158 fXS = pv->Value(ekin); << 125 xs = pv->Value(ekin);
159 // high energy CHIPS parameterisation << 126 } else if(ekin >= rTransitionBound){
160 } else if(ekin >= eTransitionBound) { << 127 xs = ggXsection->GetElementCrossSection(aParticle, Z, mat);
161 fXS = ggXsection->ComputeElementXSection(e <<
162 // linear interpolation <<
163 } else { 128 } else {
164 const G4double rxs = xs150[Z]; 129 const G4double rxs = xs150[Z];
165 const G4double lxs = pv->Value(emax); 130 const G4double lxs = pv->Value(emax);
166 fXS = lxs + (ekin - emax)*(rxs - lxs)/(eTr << 131 xs = lxs + (ekin - emax)*(rxs - lxs)/(rTransitionBound-emax);
167 } 132 }
168 133
169 #ifdef G4VERBOSE 134 #ifdef G4VERBOSE
170 if(verboseLevel > 1) { 135 if(verboseLevel > 1) {
171 G4cout << "Z= " << Z << " Ekin(MeV)= " << 136 G4cout << "Z= " << Z << " Ekin(MeV)= " << ekin/CLHEP::MeV
172 << ", nElmXS(b)= " << fXS/CLHEP::barn << 137 << ", nElmXS(b)= " << xs/CLHEP::barn
173 << G4endl; 138 << G4endl;
174 } 139 }
175 #endif 140 #endif
176 return fXS; << 141 return xs;
177 } 142 }
178 143
179 G4double G4GammaNuclearXS::LowEnergyCrossSecti << 144 G4double
>> 145 G4GammaNuclearXS::ElementCrossSection(G4double ekin, G4int ZZ)
180 { 146 {
181 const G4int Z = (ZZ < MAXZGAMMAXS) ? ZZ : MA << 147 G4DynamicParticle theGamma(gamma, G4ThreeVector(1,0,0), ekin);
182 auto pv = data->GetElementData(Z); << 148 return GetElementCrossSection(&theGamma, ZZ);
183 return pv->Value(ekin); <<
184 } 149 }
185 150
186 G4double G4GammaNuclearXS::GetIsoCrossSection( << 151 G4double
187 const G4DynamicParticle* aParticle, << 152 G4GammaNuclearXS::IsoCrossSection(G4double ekin, G4int ZZ, G4int A)
188 G4int Z, G4int A, <<
189 const G4Isotope*, const G4Element*, const G <<
190 { 153 {
191 return IsoCrossSection(aParticle->GetKinetic << 154 G4DynamicParticle theGamma(gamma, G4ThreeVector(1,0,0), ekin);
>> 155 return GetIsoCrossSection(&theGamma, ZZ, A);
192 } 156 }
193 157
194 G4double << 158 G4double G4GammaNuclearXS::GetIsoCrossSection(
195 G4GammaNuclearXS::IsoCrossSection(const G4doub << 159 const G4DynamicParticle* aParticle,
>> 160 G4int ZZ, G4int A,
>> 161 const G4Isotope*, const G4Element*,
>> 162 const G4Material*)
196 { 163 {
197 const G4int Z = (ZZ < MAXZGAMMAXS) ? ZZ : MA << 164 const G4int Z = (ZZ >= MAXZGAMMAXS) ? MAXZGAMMAXS - 1 : ZZ;
198 // cross section per element << 165 /*
199 G4double xs = ElementCrossSection(ekin, Z); << 166 G4cout << "IsoCrossSection Z= " << Z << " A= " << A
>> 167 << " Amin= " << amin[Z] << " Amax= " << amax[Z]
>> 168 << " E(MeV)= " << ekin << G4endl;
>> 169 */
>> 170 auto pv = GetPhysicsVector(Z);
>> 171 if(pv == nullptr) {
>> 172 return ggXsection->GetIsoCrossSection(aParticle, Z, A);
>> 173 }
>> 174 const G4double ekin = aParticle->GetKineticEnergy();
>> 175 const G4double emax = pv->GetMaxEnergy();
>> 176 G4double xs = 0.0;
200 177
201 if (Z > 2) { << 178 // compute isotope cross section if applicable
202 xs *= A/aeff[Z]; << 179 if(amin[Z] < amax[Z] && A >= amin[Z] && A <= amax[Z] &&
203 } else { << 180 ekin < rTransitionBound) {
204 G4int AA = A - amin[Z]; << 181 auto pviso = data->GetComponentDataByIndex(Z, A - amin[Z]);
205 if(ekin >= ehigh && AA >=0 && AA <=2) { << 182 // isotope file exists
206 xs *= coeff[Z][AA]; << 183 if(nullptr != pviso) {
207 } else { << 184 const G4double emaxiso = pviso->GetMaxEnergy();
208 xs = ggXsection->ComputeIsoXSection(ekin << 185 if(ekin <= emaxiso) {
>> 186 xs = pviso->Value(ekin);
>> 187 } else {
>> 188 G4DynamicParticle
>> 189 theGamma(gamma, G4ThreeVector(0,0,1.), rTransitionBound);
>> 190 const G4double rxs = ggXsection->GetIsoCrossSection(&theGamma, Z, A);
>> 191 const G4double lxs = pviso->Value(emaxiso);
>> 192 xs = lxs + (ekin - emaxiso)*(rxs - lxs)/(rTransitionBound-emaxiso);
>> 193 }
>> 194 #ifdef G4VERBOSE
>> 195 if(verboseLevel > 1) {
>> 196 G4cout << "G4GammaNuclearXS::IsoXS: Z= " << Z << " A= " << A
>> 197 << " Ekin(MeV)= " << ekin/CLHEP::MeV
>> 198 << ", ElmXS(b)= " << xs/CLHEP::barn << G4endl;
>> 199 }
>> 200 #endif
>> 201 return xs;
209 } 202 }
210 } 203 }
211 204
>> 205 // use element x-section
>> 206 // for the hydrogen target there is no element data
>> 207 if(ekin <= emax && Z != 1) {
>> 208 xs = pv->Value(ekin)*A/aeff[Z];
>> 209
>> 210 // CHIPS for high energy and for the hydrogen target
>> 211 } else if(ekin >= rTransitionBound || Z == 1) {
>> 212 if(Z <= 2 && ekin > 10.*GeV) {
>> 213 xs = coeff[Z][A - amin[Z]]*
>> 214 ggXsection->GetElementCrossSection(aParticle, Z, 0);
>> 215 } else {
>> 216 xs = ggXsection->GetIsoCrossSection(aParticle, Z, A);
>> 217 }
>> 218
>> 219 // transition GDR to CHIPS
>> 220 } else {
>> 221 const G4double rxs = xs150[Z];
>> 222 const G4double lxs = pv->Value(emax)*A/aeff[Z];
>> 223 xs = lxs + (ekin - emax)*(rxs - lxs)/(rTransitionBound-emax);
>> 224 }
212 #ifdef G4VERBOSE 225 #ifdef G4VERBOSE
213 if(verboseLevel > 1) { 226 if(verboseLevel > 1) {
214 G4cout << "G4GammaNuclearXS::IsoXS: Z= " 227 G4cout << "G4GammaNuclearXS::IsoXS: Z= " << Z << " A= " << A
215 << " Ekin(MeV)= " << ekin/CLHEP::MeV 228 << " Ekin(MeV)= " << ekin/CLHEP::MeV
216 << ", ElmXS(b)= " << xs/CLHEP::barn << G 229 << ", ElmXS(b)= " << xs/CLHEP::barn << G4endl;
217 } 230 }
218 #endif 231 #endif
219 return xs; 232 return xs;
220 } 233 }
221 234
222 const G4Isotope* G4GammaNuclearXS::SelectIsoto 235 const G4Isotope* G4GammaNuclearXS::SelectIsotope(
223 const G4Element* anElement, G4double ki 236 const G4Element* anElement, G4double kinEnergy, G4double)
224 { 237 {
225 std::size_t nIso = anElement->GetNumberOfIso << 238 size_t nIso = anElement->GetNumberOfIsotopes();
226 const G4Isotope* iso = anElement->GetIsotope 239 const G4Isotope* iso = anElement->GetIsotope(0);
227 240
228 if(1 == nIso) { return iso; } 241 if(1 == nIso) { return iso; }
229 242
230 const G4double* abundVector = anElement->Get 243 const G4double* abundVector = anElement->GetRelativeAbundanceVector();
>> 244 G4double q = G4UniformRand();
231 G4double sum = 0.0; 245 G4double sum = 0.0;
>> 246 size_t j;
232 G4int Z = anElement->GetZasInt(); 247 G4int Z = anElement->GetZasInt();
233 248
>> 249 // condition to use only isotope abundance
>> 250 if(amax[Z] == amin[Z] || kinEnergy > rTransitionBound || Z >= MAXZGAMMAXS ) {
>> 251 for (j=0; j<nIso; ++j) {
>> 252 sum += abundVector[j];
>> 253 if(q <= sum) {
>> 254 iso = anElement->GetIsotope(j);
>> 255 break;
>> 256 }
>> 257 }
>> 258 return iso;
>> 259 }
234 // use isotope cross sections 260 // use isotope cross sections
235 std::size_t nn = temp.size(); << 261 size_t nn = temp.size();
236 if(nn < nIso) { temp.resize(nIso, 0.); } 262 if(nn < nIso) { temp.resize(nIso, 0.); }
237 263
238 for (std::size_t j=0; j<nIso; ++j) { << 264 for (j=0; j<nIso; ++j) {
239 //G4cout << j << "-th isotope " << (*isoVe 265 //G4cout << j << "-th isotope " << (*isoVector)[j]->GetN()
240 // << " abund= " << abundVector[j] 266 // << " abund= " << abundVector[j] << G4endl;
241 sum += abundVector[j]* 267 sum += abundVector[j]*
242 IsoCrossSection(kinEnergy, Z, anElement- << 268 IsoCrossSection(kinEnergy, Z, anElement->GetIsotope(j)->GetN());
243 temp[j] = sum; 269 temp[j] = sum;
244 } 270 }
245 sum *= G4UniformRand(); << 271 sum *= q;
246 for (std::size_t j = 0; j<nIso; ++j) { << 272 for (j = 0; j<nIso; ++j) {
247 if(temp[j] >= sum) { 273 if(temp[j] >= sum) {
248 iso = anElement->GetIsotope((G4int)j); << 274 iso = anElement->GetIsotope(j);
249 break; 275 break;
250 } 276 }
251 } 277 }
252 return iso; 278 return iso;
253 } 279 }
254 280
255 void G4GammaNuclearXS::BuildPhysicsTable(const << 281 void
>> 282 G4GammaNuclearXS::BuildPhysicsTable(const G4ParticleDefinition& p)
256 { 283 {
257 if(verboseLevel > 1) { << 284 if(verboseLevel > 0){
258 G4cout << "G4GammaNuclearXS::BuildPhysicsT 285 G4cout << "G4GammaNuclearXS::BuildPhysicsTable for "
259 << p.GetParticleName() << G4endl; 286 << p.GetParticleName() << G4endl;
260 } 287 }
261 if(p.GetParticleName() != "gamma") { 288 if(p.GetParticleName() != "gamma") {
262 G4ExceptionDescription ed; 289 G4ExceptionDescription ed;
263 ed << p.GetParticleName() << " is a wrong 290 ed << p.GetParticleName() << " is a wrong particle type -"
264 << " only gamma is allowed"; 291 << " only gamma is allowed";
265 G4Exception("G4GammaNuclearXS::BuildPhysic 292 G4Exception("G4GammaNuclearXS::BuildPhysicsTable(..)","had012",
266 FatalException, ed, ""); 293 FatalException, ed, "");
267 return; 294 return;
268 } 295 }
269 296
270 // prepare isotope selection << 297 if(nullptr == data) {
>> 298 #ifdef G4MULTITHREADED
>> 299 G4MUTEXLOCK(&gNuclearXSMutex);
>> 300 if(nullptr == data) {
>> 301 #endif
>> 302 isMaster = true;
>> 303 data = new G4ElementData();
>> 304 data->SetName("PhotoNuclear");
>> 305 FindDirectoryPath();
>> 306 #ifdef G4MULTITHREADED
>> 307 }
>> 308 G4MUTEXUNLOCK(&gNuclearXSMutex);
>> 309 #endif
>> 310 }
>> 311
>> 312
>> 313 // it is possible re-initialisation for the second run
>> 314 // Upload data for elements used in geometry
271 const G4ElementTable* table = G4Element::Get 315 const G4ElementTable* table = G4Element::GetElementTable();
272 std::size_t nIso = temp.size(); << 316 if(isMaster) {
273 for (auto const & elm : *table ) { << 317 for ( auto & elm : *table ) {
274 std::size_t n = elm->GetNumberOfIsotopes() << 318 G4int Z = std::max( 1, std::min( elm->GetZasInt(), MAXZGAMMAXS-1) );
275 if (n > nIso) { nIso = n; } << 319 if ( nullptr == data->GetElementData(Z) ) { Initialise(Z); }
>> 320 }
>> 321 }
>> 322
>> 323 // prepare isotope selection
>> 324 size_t nIso = temp.size();
>> 325 for ( auto & elm : *table ) {
>> 326 size_t n = elm->GetNumberOfIsotopes();
>> 327 if(n > nIso) { nIso = n; }
276 } 328 }
277 temp.resize(nIso, 0.0); 329 temp.resize(nIso, 0.0);
278 } 330 }
279 331
280 const G4String& G4GammaNuclearXS::FindDirector 332 const G4String& G4GammaNuclearXS::FindDirectoryPath()
281 { 333 {
>> 334 // check environment variable
282 // build the complete string identifying the 335 // build the complete string identifying the file with the data set
283 if(gDataDirectory.empty()) { 336 if(gDataDirectory.empty()) {
284 std::ostringstream ost; << 337 char* path = std::getenv("G4PARTICLEXSDATA");
285 ost << G4HadronicParameters::Instance()->G << 338 if (nullptr != path) {
286 gDataDirectory = ost.str(); << 339 std::ostringstream ost;
>> 340 ost << path << "/gamma/inel";
>> 341 gDataDirectory = ost.str();
>> 342 } else {
>> 343 G4Exception("G4GammaNuclearXS::Initialise(..)","had013",
>> 344 FatalException,
>> 345 "Environment variable G4PARTICLEXSDATA is not defined");
>> 346 }
287 } 347 }
288 return gDataDirectory; 348 return gDataDirectory;
289 } 349 }
290 350
>> 351 void G4GammaNuclearXS::InitialiseOnFly(G4int Z)
>> 352 {
>> 353 #ifdef G4MULTITHREADED
>> 354 G4MUTEXLOCK(&gNuclearXSMutex);
>> 355 if(nullptr == data->GetElementData(Z)) {
>> 356 #endif
>> 357 Initialise(Z);
>> 358 #ifdef G4MULTITHREADED
>> 359 }
>> 360 G4MUTEXUNLOCK(&gNuclearXSMutex);
>> 361 #endif
>> 362 }
>> 363
291 void G4GammaNuclearXS::Initialise(G4int Z) 364 void G4GammaNuclearXS::Initialise(G4int Z)
292 { 365 {
>> 366 if(nullptr != data->GetElementData(Z)) { return; }
>> 367
293 // upload data from file 368 // upload data from file
294 std::ostringstream ost; 369 std::ostringstream ost;
295 ost << FindDirectoryPath() << Z ; 370 ost << FindDirectoryPath() << Z ;
296 G4PhysicsVector* v = RetrieveVector(ost, tru 371 G4PhysicsVector* v = RetrieveVector(ost, true, Z);
297 372
298 data->InitialiseForElement(Z, v); 373 data->InitialiseForElement(Z, v);
299 /* 374 /*
300 G4cout << "G4GammaNuclearXS::Initialise for 375 G4cout << "G4GammaNuclearXS::Initialise for Z= " << Z
301 << " A= " << Amean << " Amin= " << amin[Z] 376 << " A= " << Amean << " Amin= " << amin[Z]
302 << " Amax= " << amax[Z] << G4endl; 377 << " Amax= " << amax[Z] << G4endl;
303 */ 378 */
304 xs150[Z] = ggXsection->ComputeElementXSectio << 379 // upload isotope data
305 << 380 G4DynamicParticle theGamma(gamma, G4ThreeVector(1,0,0), rTransitionBound);
306 // compute corrections for low Z data << 381 xs150[Z] = ggXsection->GetElementCrossSection(&theGamma, Z, 0);
307 if(Z <= 2){ << 382 if(amax[Z] > amin[Z]) {
308 if(amax[Z] > amin[Z]) { << 383 size_t nmax = (size_t)(amax[Z]-amin[Z]+1);
309 for(G4int A=amin[Z]; A<=amax[Z]; ++A) { << 384 data->InitialiseForComponent(Z, nmax);
310 G4int AA = A - amin[Z]; << 385 for(G4int A=amin[Z]; A<=amax[Z]; ++A) {
311 if(AA >= 0 && AA <= 2) { << 386 std::ostringstream ost1;
312 G4double sig1 = ggXsection->ComputeIsoXSec << 387 ost1 << gDataDirectory << Z << "_" << A;
313 G4double sig2 = ggXsection->ComputeElement << 388 G4PhysicsVector* v1 = RetrieveVector(ost1, false, Z);
314 if(sig2 > 0.) { coeff[Z][AA] = (sig1/sig2) << 389 data->AddComponent(Z, A, v1);
315 else { coeff[Z][AA] = 1.; } << 390 if(Z<=2){
316 } << 391 theGamma.SetKineticEnergy(10.*GeV);
>> 392 G4double sig1 = ggXsection->GetIsoCrossSection(&theGamma, Z, A);
>> 393 G4double sig2 = ggXsection->GetElementCrossSection(&theGamma, Z, 0);
>> 394 if(sig2 > 0.) coeff[Z][A-amin[Z]]=(sig1/sig2);
>> 395 else coeff[Z][A-amin[Z]]=1.;
317 } 396 }
318 } 397 }
319 } 398 }
320 } 399 }
321 400
322 G4PhysicsVector* 401 G4PhysicsVector*
323 G4GammaNuclearXS::RetrieveVector(std::ostrings << 402 G4GammaNuclearXS::RetrieveVector(std::ostringstream& ost, G4bool isElement, G4int Z)
324 { 403 {
325 G4PhysicsVector* v = nullptr; 404 G4PhysicsVector* v = nullptr;
326 405
327 std::ifstream filein(ost.str().c_str()); 406 std::ifstream filein(ost.str().c_str());
328 if (!filein.is_open()) { 407 if (!filein.is_open()) {
329 if(warn) { << 408 if(isElement) {
330 G4ExceptionDescription ed; 409 G4ExceptionDescription ed;
331 ed << "Data file <" << ost.str().c_str() 410 ed << "Data file <" << ost.str().c_str()
332 << "> is not opened!"; 411 << "> is not opened!";
333 G4Exception("G4GammaNuclearXS::RetrieveV 412 G4Exception("G4GammaNuclearXS::RetrieveVector(..)","had014",
334 FatalException, ed, "Check G4PARTICLEXSD 413 FatalException, ed, "Check G4PARTICLEXSDATA");
335 } 414 }
336 } else { 415 } else {
337 if(verboseLevel > 1) { 416 if(verboseLevel > 1) {
338 G4cout << "File " << ost.str() 417 G4cout << "File " << ost.str()
339 << " is opened by G4GammaNuclearXS" << 418 << " is opened by G4GammaNuclearXS" << G4endl;
340 } 419 }
341 // retrieve data from DB 420 // retrieve data from DB
342 if(std::find(std::begin(freeVectorExceptio << 421 if(std::find(std::begin(freeVectorException), std::end(freeVectorException), Z ) == std::end(freeVectorException) && isElement) {
343 == std::end(freeVectorException)) { << 422 v = new G4PhysicsLinearVector();
344 v = new G4PhysicsLinearVector(false); <<
345 } else { 423 } else {
346 v = new G4PhysicsFreeVector(false); << 424 v = new G4PhysicsVector();
347 } 425 }
348 if(!v->Retrieve(filein, true)) { 426 if(!v->Retrieve(filein, true)) {
349 G4ExceptionDescription ed; 427 G4ExceptionDescription ed;
350 ed << "Data file <" << ost.str().c_str() 428 ed << "Data file <" << ost.str().c_str()
351 << "> is not retrieved!"; 429 << "> is not retrieved!";
352 G4Exception("G4GammaNuclearXS::RetrieveV 430 G4Exception("G4GammaNuclearXS::RetrieveVector(..)","had015",
353 FatalException, ed, "Check G4PARTICLEXSD 431 FatalException, ed, "Check G4PARTICLEXSDATA");
354 } 432 }
355 } 433 }
356 return v; << 434 return v;
357 } 435 }
358 436
359 437