| 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 // Author: Sebastien Incerti 27 // Author: Sebastien Incerti
28 // 30 October 2008 28 // 30 October 2008
29 // on base of G4LowEnergyPhotoElectric 29 // on base of G4LowEnergyPhotoElectric developed by A.Forti and M.G.Pia
30 // 30 //
31 // 22 Oct 2012 A & V Ivanchenko Migration da 31 // 22 Oct 2012 A & V Ivanchenko Migration data structure to G4PhysicsVector
32 // 1 June 2017 M Bandieramonte: New model ba << 32 // 1 June 2017 M Bandieramonte: New model based on livermore/epics2014
33 // evaluated data - parameteriza 33 // evaluated data - parameterization fits in two ranges
34 34
35 #include "G4LivermorePhotoElectricModel.hh" 35 #include "G4LivermorePhotoElectricModel.hh"
36 << 36 #include "G4SystemOfUnits.hh"
37 #include "G4AtomicShell.hh" << 37 #include "G4PhysicalConstants.hh"
38 #include "G4AutoLock.hh" << 38 #include "G4LossTableManager.hh"
39 #include "G4CrossSectionHandler.hh" <<
40 #include "G4Electron.hh" 39 #include "G4Electron.hh"
41 #include "G4Gamma.hh" 40 #include "G4Gamma.hh"
42 #include "G4LossTableManager.hh" <<
43 #include "G4ParticleChangeForGamma.hh" 41 #include "G4ParticleChangeForGamma.hh"
44 #include "G4PhysicsFreeVector.hh" << 42 #include "G4CrossSectionHandler.hh"
45 #include "G4SauterGavrilaAngularDistribution.h << 43 #include "G4LPhysicsFreeVector.hh"
46 #include "G4SystemOfUnits.hh" <<
47 #include "G4VAtomDeexcitation.hh" 44 #include "G4VAtomDeexcitation.hh"
48 #include "G4EmParameters.hh" << 45 #include "G4SauterGavrilaAngularDistribution.hh"
49 #include <thread> << 46 #include "G4AtomicShell.hh"
50 47
51 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
52 49
53 G4ElementData* G4LivermorePhotoElectricModel:: << 50 G4LPhysicsFreeVector* G4LivermorePhotoElectricModel::fCrossSection[] = {nullptr};
54 G4ElementData* G4LivermorePhotoElectricModel:: << 51 G4LPhysicsFreeVector* G4LivermorePhotoElectricModel::fCrossSectionLE[] = {nullptr};
55 std::vector<G4double>* G4LivermorePhotoElectri 52 std::vector<G4double>* G4LivermorePhotoElectricModel::fParamHigh[] = {nullptr};
56 std::vector<G4double>* G4LivermorePhotoElectri 53 std::vector<G4double>* G4LivermorePhotoElectricModel::fParamLow[] = {nullptr};
57 G4int G4LivermorePhotoElectricModel::fNShells[ << 54 G4int G4LivermorePhotoElectricModel::fNShells[] = {0};
58 G4int G4LivermorePhotoElectricModel::fNShellsU << 55 G4int G4LivermorePhotoElectricModel::fNShellsUsed[] = {0};
59 G4Material* G4LivermorePhotoElectricModel::fWa << 56 G4ElementData* G4LivermorePhotoElectricModel::fShellCrossSection = nullptr;
60 G4double G4LivermorePhotoElectricModel::fWater << 57 G4Material* G4LivermorePhotoElectricModel::fWater = nullptr;
61 G4String G4LivermorePhotoElectricModel::fDataD << 58 G4double G4LivermorePhotoElectricModel::fWaterEnergyLimit = 0.0;
>> 59 G4String G4LivermorePhotoElectricModel::fDataDirectory = "";
>> 60
>> 61 #ifdef G4MULTITHREADED
>> 62 G4Mutex G4LivermorePhotoElectricModel::livPhotoeffMutex = G4MUTEX_INITIALIZER;
>> 63 #endif
62 64
63 static std::once_flag applyOnce; << 65 using namespace std;
64 <<
65 namespace <<
66 { <<
67 G4Mutex livPhotoeffMutex = G4MUTEX_INITIALIZ <<
68 } <<
69 66
70 //....oooOO0OOooo........oooOO0OOooo........oo 67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 68
72 G4LivermorePhotoElectricModel::G4LivermorePhot << 69 G4LivermorePhotoElectricModel::G4LivermorePhotoElectricModel(const G4String& nam)
73 { << 70 : G4VEmModel(nam),fParticleChange(nullptr),maxZ(99),
74 verboseLevel = 0; << 71 nShellLimit(100),fDeexcitationActive(false),isInitialised(false),
75 // Verbosity scale: << 72 fAtomDeexcitation(nullptr)
76 // 0 = nothing << 73 {
77 // 1 = warning for energy non-conservation << 74 verboseLevel= 0;
78 // 2 = details of energy budget << 75 // Verbosity scale:
79 // 3 = calculation of cross sections, file o << 76 // 0 = nothing
80 // 4 = entering in methods << 77 // 1 = warning for energy non-conservation
81 << 78 // 2 = details of energy budget
82 theGamma = G4Gamma::Gamma(); << 79 // 3 = calculation of cross sections, file openings, sampling of atoms
83 theElectron = G4Electron::Electron(); << 80 // 4 = entering in methods
84 << 81
85 // default generator << 82 theGamma = G4Gamma::Gamma();
86 SetAngularDistribution(new G4SauterGavrilaAn << 83 theElectron = G4Electron::Electron();
87 << 84
88 if (verboseLevel > 0) { << 85 // default generator
89 G4cout << "Livermore PhotoElectric is cons << 86 SetAngularDistribution(new G4SauterGavrilaAngularDistribution());
90 << " nShellLimit= " << nShellLimit << 87
91 } << 88 if(verboseLevel>0) {
92 << 89 G4cout << "Livermore PhotoElectric is constructed "
93 // Mark this model as "applicable" for atomi << 90 << " nShellLimit= " << nShellLimit << G4endl;
94 SetDeexcitationFlag(true); << 91 }
95 << 92
96 // For water << 93 //Mark this model as "applicable" for atomic deexcitation
97 fSandiaCof.resize(4, 0.0); << 94 SetDeexcitationFlag(true);
>> 95 fSandiaCof.resize(4,0.0);
>> 96 fCurrSection = 0.0;
98 } 97 }
99 98
100 //....oooOO0OOooo........oooOO0OOooo........oo 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 100
102 G4LivermorePhotoElectricModel::~G4LivermorePho 101 G4LivermorePhotoElectricModel::~G4LivermorePhotoElectricModel()
103 { 102 {
104 if (isInitializer) { << 103 if(IsMaster()) {
105 for (G4int i = 0; i < ZMAXPE; ++i) { << 104 delete fShellCrossSection;
106 if (fParamHigh[i]) { << 105 fShellCrossSection = nullptr;
107 delete fParamHigh[i]; << 106 for(G4int i=0; i<maxZ; ++i) {
108 fParamHigh[i] = nullptr; << 107 delete fParamHigh[i];
109 } << 108 fParamHigh[i] = nullptr;
110 if (fParamLow[i]) { << 109 delete fParamLow[i];
111 delete fParamLow[i]; << 110 fParamLow[i] = nullptr;
112 fParamLow[i] = nullptr; << 111 delete fCrossSection[i];
113 } << 112 fCrossSection[i] = nullptr;
>> 113 delete fCrossSectionLE[i];
>> 114 fCrossSectionLE[i] = nullptr;
>> 115 }
114 } 116 }
115 } <<
116 } 117 }
117 118
118 //....oooOO0OOooo........oooOO0OOooo........oo 119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119 120
120 void G4LivermorePhotoElectricModel::Initialise << 121 void
121 const G4DataVector&) << 122 G4LivermorePhotoElectricModel::Initialise(const G4ParticleDefinition*,
>> 123 const G4DataVector&)
122 { 124 {
123 if (verboseLevel > 1) { << 125 if (verboseLevel > 2) {
124 G4cout << "Calling G4LivermorePhotoElectri << 126 G4cout << "Calling G4LivermorePhotoElectricModel::Initialise() " << G4endl;
125 } << 127 }
126 << 128
127 // initialise static tables only once << 129 if(IsMaster()) {
128 std::call_once(applyOnce, [this]() { isIniti << 130
129 << 131 if(!fWater) {
130 if (isInitializer) { << 132 fWater = G4Material::GetMaterial("G4_WATER", false);
131 G4AutoLock l(&livPhotoeffMutex); << 133 if(fWater) { fWaterEnergyLimit = 13.6*eV; }
132 FindDirectoryPath(); << 134 }
133 if (fWater == nullptr) { << 135
134 fWater = G4Material::GetMaterial("G4_WAT << 136 if(!fShellCrossSection) { fShellCrossSection = new G4ElementData(); }
135 if (fWater == nullptr) { << 137
136 fWater = G4Material::GetMaterial("Wate << 138 G4ProductionCutsTable* theCoupleTable =
137 } << 139 G4ProductionCutsTable::GetProductionCutsTable();
138 if (fWater != nullptr) { << 140 G4int numOfCouples = theCoupleTable->GetTableSize();
139 fWaterEnergyLimit = 13.6 * CLHEP::eV; << 141
140 } << 142 for(G4int i=0; i<numOfCouples; ++i) {
141 } << 143 const G4MaterialCutsCouple* couple =
142 << 144 theCoupleTable->GetMaterialCutsCouple(i);
143 if (fCrossSection == nullptr) { << 145 const G4Material* material = couple->GetMaterial();
144 fCrossSection = new G4ElementData(ZMAXPE << 146 const G4ElementVector* theElementVector = material->GetElementVector();
145 fCrossSection->SetName("PhotoEffXS"); << 147 G4int nelm = material->GetNumberOfElements();
146 fCrossSectionLE = new G4ElementData(ZMAX << 148
147 fCrossSectionLE->SetName("PhotoEffLowXS" << 149 for (G4int j=0; j<nelm; ++j) {
148 } << 150 G4int Z = std::min(maxZ, (*theElementVector)[j]->GetZasInt());
149 << 151 if(!fCrossSection[Z]) { ReadData(Z); }
150 const G4ElementTable* elemTable = G4Elemen << 152 }
151 std::size_t numElems = (*elemTable).size() << 153 }
152 for (std::size_t ie = 0; ie < numElems; ++ << 154 }
153 const G4Element* elem = (*elemTable)[ie] << 155
154 G4int Z = elem->GetZasInt(); << 156 if (verboseLevel > 2) {
155 if (Z < ZMAXPE) { << 157 G4cout << "Loaded cross section files for new LivermorePhotoElectric model"
156 if (fCrossSection->GetElementData(Z) == null << 158 << G4endl;
157 ReadData(Z); << 159 }
158 } << 160 if(!isInitialised) {
159 } << 161 isInitialised = true;
>> 162 fParticleChange = GetParticleChangeForGamma();
>> 163
>> 164 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
>> 165 }
>> 166 fDeexcitationActive = false;
>> 167 if(fAtomDeexcitation) {
>> 168 fDeexcitationActive = fAtomDeexcitation->IsFluoActive();
>> 169 }
>> 170
>> 171 if (verboseLevel > 0) {
>> 172 G4cout << "LivermorePhotoElectric model is initialized " << G4endl
>> 173 << G4endl;
160 } 174 }
161 l.unlock(); <<
162 } <<
163 <<
164 if (verboseLevel > 1) { <<
165 G4cout << "Loaded cross section files for <<
166 } <<
167 if (nullptr == fParticleChange) { <<
168 fParticleChange = GetParticleChangeForGamm <<
169 fAtomDeexcitation = G4LossTableManager::In <<
170 } <<
171 <<
172 fDeexcitationActive = false; <<
173 if (nullptr != fAtomDeexcitation) { <<
174 fDeexcitationActive = fAtomDeexcitation->I <<
175 } <<
176 <<
177 if (verboseLevel > 1) { <<
178 G4cout << "LivermorePhotoElectric model is <<
179 } <<
180 } 175 }
181 176
182 //....oooOO0OOooo........oooOO0OOooo........oo 177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
183 178
184 G4double << 179 G4double G4LivermorePhotoElectricModel::CrossSectionPerVolume(
185 G4LivermorePhotoElectricModel::CrossSectionPer << 180 const G4Material* material,
186 << 181 const G4ParticleDefinition* p,
187 << 182 G4double energy,
188 G4double, G4double) << 183 G4double, G4double)
189 { << 184 {
190 fCurrSection = 0.0; << 185 fCurrSection = 0.0;
191 if (fWater && (material == fWater || materia << 186 if(fWater && (material == fWater ||
192 if (energy <= fWaterEnergyLimit) { << 187 material->GetBaseMaterial() == fWater)) {
193 fWater->GetSandiaTable()->GetSandiaCofWa << 188 if(energy <= fWaterEnergyLimit) {
194 << 189 fWater->GetSandiaTable()->GetSandiaCofWater(energy, fSandiaCof);
195 G4double energy2 = energy * energy; << 190
196 G4double energy3 = energy * energy2; << 191 G4double energy2 = energy*energy;
197 G4double energy4 = energy2 * energy2; << 192 G4double energy3 = energy*energy2;
198 << 193 G4double energy4 = energy2*energy2;
199 fCurrSection = material->GetDensity() << 194
200 * (fSandiaCof[0] / energy << 195 fCurrSection = material->GetDensity()*
201 fSandiaCof[2] / energy3 + fSandiaCof[3] << 196 (fSandiaCof[0]/energy + fSandiaCof[1]/energy2 +
>> 197 fSandiaCof[2]/energy3 + fSandiaCof[3]/energy4);
>> 198 }
202 } 199 }
203 } << 200 if(0.0 == fCurrSection) {
204 if (0.0 == fCurrSection) { << 201 fCurrSection = G4VEmModel::CrossSectionPerVolume(material, p, energy);
205 fCurrSection = G4VEmModel::CrossSectionPer << 202 }
206 } << 203 return fCurrSection;
207 return fCurrSection; <<
208 } 204 }
209 205
210 //....oooOO0OOooo........oooOO0OOooo........oo 206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
211 207
212 G4double << 208 G4double G4LivermorePhotoElectricModel::ComputeCrossSectionPerAtom(
213 G4LivermorePhotoElectricModel::ComputeCrossSec << 209 const G4ParticleDefinition*,
214 << 210 G4double energy,
215 G4double ZZ, G4double, << 211 G4double ZZ, G4double,
216 G4double, G4double) << 212 G4double, G4double)
217 { 213 {
218 if (verboseLevel > 3) { << 214 if (verboseLevel > 3) {
219 G4cout << "\n G4LivermorePhotoElectricMode << 215 G4cout << "\n G4LivermorePhotoElectricModel::ComputeCrossSectionPerAtom():"
220 << " Z= " << ZZ << " R(keV)= " << << 216 << " Z= " << ZZ << " R(keV)= " << energy/keV << G4endl;
221 } << 217 }
222 G4double cs = 0.0; << 218 G4double cs = 0.0;
223 G4int Z = G4lrint(ZZ); << 219 G4int Z = G4lrint(ZZ);
224 if (Z >= ZMAXPE || Z <= 0) { << 220 if(Z >= maxZ) { return cs; }
>> 221 // if element was not initialised
>> 222
>> 223 // do initialisation safely for MT mode
>> 224 if(!fCrossSection[Z]) { InitialiseForElement(theGamma, Z); }
>> 225
>> 226 //7: rows in the parameterization file; 5: number of parameters
>> 227 G4int idx = fNShells[Z]*7 - 5;
>> 228
>> 229 energy = std::max(energy, (*(fParamHigh[Z]))[idx-1]);
>> 230
>> 231 G4double x1 = 1.0/energy;
>> 232 G4double x2 = x1*x1;
>> 233 G4double x3 = x2*x1;
>> 234
>> 235 // high energy parameterisation
>> 236 if(energy >= (*(fParamHigh[Z]))[0]) {
>> 237
>> 238 G4double x4 = x2*x2;
>> 239 G4double x5 = x4*x1;
>> 240
>> 241 cs = x1*((*(fParamHigh[Z]))[idx] + x1*(*(fParamHigh[Z]))[idx+1]
>> 242 + x2*(*(fParamHigh[Z]))[idx+2] + x3*(*(fParamHigh[Z]))[idx+3]
>> 243 + x4*(*(fParamHigh[Z]))[idx+4]+ x5*(*(fParamHigh[Z]))[idx+5]);
>> 244
>> 245 }
>> 246 // low energy parameterisation
>> 247 else if(energy >= (*(fParamLow[Z]))[0]) {
>> 248
>> 249 G4double x4 = x2*x2;
>> 250 G4double x5 = x4*x1; //this variable usage can probably be optimized
>> 251 cs = x1*((*(fParamLow[Z]))[idx] + x1*(*(fParamLow[Z]))[idx+1]
>> 252 + x2*(*(fParamLow[Z]))[idx+2] + x3*(*(fParamLow[Z]))[idx+3]
>> 253 + x4*(*(fParamLow[Z]))[idx+4]+ x5*(*(fParamLow[Z]))[idx+5]);
>> 254
>> 255 }
>> 256 // Tabulated values above k-shell ionization energy
>> 257 else if(energy >= (*(fParamHigh[Z]))[1]) {
>> 258 cs = x3*(fCrossSection[Z])->Value(energy);
>> 259 }
>> 260 // Tabulated values below k-shell ionization energy
>> 261 else
>> 262 {
>> 263 cs = x3*(fCrossSectionLE[Z])->Value(energy);
>> 264 }
>> 265 if (verboseLevel > 1) {
>> 266 G4cout << "G4LivermorePhotoElectricModel: E(keV)= " << energy/keV
>> 267 << " Z= " << Z << " cross(barn)= " << cs/barn << G4endl;
>> 268 }
225 return cs; 269 return cs;
226 } <<
227 // if element was not initialised <<
228 <<
229 // do initialisation safely for MT mode <<
230 if (fCrossSection->GetElementData(Z) == null <<
231 InitialiseOnFly(Z); <<
232 if (fCrossSection->GetElementData(Z) == nu <<
233 } <<
234 <<
235 // 7: rows in the parameterization file; 5: <<
236 G4int idx = fNShells[Z] * 7 - 5; <<
237 <<
238 energy = std::max(energy, (*(fParamHigh[Z])) <<
239 <<
240 G4double x1 = 1.0 / energy; <<
241 G4double x2 = x1 * x1; <<
242 G4double x3 = x2 * x1; <<
243 <<
244 // high energy parameterisation <<
245 if (energy >= (*(fParamHigh[Z]))[0]) { <<
246 G4double x4 = x2 * x2; <<
247 G4double x5 = x4 * x1; <<
248 <<
249 cs = x1 * <<
250 ((*(fParamHigh[Z]))[idx] + x1 * (*(fPara <<
251 + x2 * (*(fParamHigh[Z]))[idx + 2] + x3 <<
252 + x4 * (*(fParamHigh[Z]))[idx + 4] + x5 <<
253 } <<
254 // low energy parameterisation <<
255 else if (energy >= (*(fParamLow[Z]))[0]) { <<
256 G4double x4 = x2 * x2; <<
257 G4double x5 = x4 * x1; // this variable u <<
258 cs = x1 * <<
259 ((*(fParamLow[Z]))[idx] + x1 * (*(fParam <<
260 + x2 * (*(fParamLow[Z]))[idx + 2] + x3 <<
261 + x4 * (*(fParamLow[Z]))[idx + 4] + x5 <<
262 } <<
263 // Tabulated values above k-shell ionization <<
264 else if (energy >= (*(fParamHigh[Z]))[1]) { <<
265 cs = x3 * (fCrossSection->GetElementData(Z <<
266 } <<
267 // Tabulated values below k-shell ionization <<
268 else { <<
269 cs = x3 * (fCrossSectionLE->GetElementData <<
270 } <<
271 if (verboseLevel > 1) { <<
272 G4cout << "G4LivermorePhotoElectricModel: <<
273 << " cross(barn)= " << cs / barn << <<
274 } <<
275 return cs; <<
276 } 270 }
277 271
278 //....oooOO0OOooo........oooOO0OOooo........oo 272 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
279 273
280 void G4LivermorePhotoElectricModel::SampleSeco << 274 void
281 << 275 G4LivermorePhotoElectricModel::SampleSecondaries(
282 << 276 std::vector<G4DynamicParticle*>* fvect,
283 << 277 const G4MaterialCutsCouple* couple,
>> 278 const G4DynamicParticle* aDynamicGamma,
>> 279 G4double, G4double)
284 { 280 {
285 G4double gammaEnergy = aDynamicGamma->GetKin << 281 G4double gammaEnergy = aDynamicGamma->GetKineticEnergy();
286 if (verboseLevel > 3) { << 282 if (verboseLevel > 3) {
287 G4cout << "G4LivermorePhotoElectricModel:: << 283 G4cout << "G4LivermorePhotoElectricModel::SampleSecondaries() Egamma(keV)= "
288 << gammaEnergy / keV << G4endl; << 284 << gammaEnergy/keV << G4endl;
289 } <<
290 <<
291 // kill incident photon <<
292 fParticleChange->ProposeTrackStatus(fStopAnd <<
293 fParticleChange->SetProposedKineticEnergy(0. <<
294 <<
295 // low-energy photo-effect in water - full a <<
296 const G4Material* material = couple->GetMate <<
297 if (fWater && (material == fWater || materia <<
298 if (gammaEnergy <= fWaterEnergyLimit) { <<
299 fParticleChange->ProposeLocalEnergyDepos <<
300 return; <<
301 } 285 }
302 } << 286
303 << 287 // kill incident photon
304 // Returns the normalized direction of the m << 288 fParticleChange->ProposeTrackStatus(fStopAndKill);
305 G4ThreeVector photonDirection = aDynamicGamm << 289 fParticleChange->SetProposedKineticEnergy(0.);
306 << 290
307 // Select randomly one element in the curren << 291 // low-energy photo-effect in water - full absorption
308 const G4Element* elm = SelectRandomAtom(mate << 292 const G4Material* material = couple->GetMaterial();
309 G4int Z = elm->GetZasInt(); << 293 if(fWater && (material == fWater ||
310 << 294 material->GetBaseMaterial() == fWater)) {
311 // Select the ionised shell in the current a << 295 if(gammaEnergy <= fWaterEnergyLimit) {
312 // cross sections << 296 fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy);
313 << 297 return;
314 // If element was not initialised gamma shou <<
315 if (Z >= ZMAXPE || fCrossSection->GetElement <<
316 fParticleChange->ProposeLocalEnergyDeposit <<
317 return; <<
318 } <<
319 <<
320 // SAMPLING OF THE SHELL INDEX <<
321 std::size_t shellIdx = 0; <<
322 std::size_t nn = fNShellsUsed[Z]; <<
323 if (nn > 1) { <<
324 if (gammaEnergy >= (*(fParamHigh[Z]))[0]) <<
325 G4double x1 = 1.0 / gammaEnergy; <<
326 G4double x2 = x1 * x1; <<
327 G4double x3 = x2 * x1; <<
328 G4double x4 = x3 * x1; <<
329 G4double x5 = x4 * x1; <<
330 std::size_t idx = nn * 7 - 5; <<
331 // when do sampling common factors are n <<
332 // so cross section is not real <<
333 <<
334 G4double rand = G4UniformRand(); <<
335 G4double cs0 = rand * <<
336 ((*(fParamHigh[Z]))[idx] + x1 * (*(fPa <<
337 + x2 * (*(fParamHigh[Z]))[idx + 2] + x <<
338 + x4 * (*(fParamHigh[Z]))[idx + 4] + x <<
339 <<
340 for (shellIdx = 0; shellIdx < nn; ++shel <<
341 idx = shellIdx * 7 + 2; <<
342 if (gammaEnergy > (*(fParamHigh[Z]))[i <<
343 G4double cs = (*(fParamHigh[Z]))[idx <<
344 + x2 * (*(fParamHigh[Z]))[idx + 2] <<
345 + x4 * (*(fParamHigh[Z]))[idx + 4] <<
346 <<
347 if (cs >= cs0) { <<
348 break; <<
349 } <<
350 } <<
351 } <<
352 if (shellIdx >= nn) { <<
353 shellIdx = nn - 1; <<
354 } <<
355 } <<
356 else if (gammaEnergy >= (*(fParamLow[Z]))[ <<
357 G4double x1 = 1.0 / gammaEnergy; <<
358 G4double x2 = x1 * x1; <<
359 G4double x3 = x2 * x1; <<
360 G4double x4 = x3 * x1; <<
361 G4double x5 = x4 * x1; <<
362 std::size_t idx = nn * 7 - 5; <<
363 // when do sampling common factors are n <<
364 // so cross section is not real <<
365 G4double cs0 = G4UniformRand() * <<
366 ((*(fParamLow[Z]))[idx] + x1 * (*(fPar <<
367 + x2 * (*(fParamLow[Z]))[idx + 2] + x3 <<
368 + x4 * (*(fParamLow[Z]))[idx + 4] + x5 <<
369 for (shellIdx = 0; shellIdx < nn; ++shel <<
370 idx = shellIdx * 7 + 2; <<
371 if (gammaEnergy > (*(fParamLow[Z]))[id <<
372 G4double cs = (*(fParamLow[Z]))[idx] <<
373 + x2 * (*(fParamLow[Z]))[idx + 2] <<
374 + x4 * (*(fParamLow[Z]))[idx + 4] <<
375 if (cs >= cs0) { <<
376 break; <<
377 } <<
378 } <<
379 } <<
380 if (shellIdx >= nn) { <<
381 shellIdx = nn - 1; <<
382 } <<
383 } <<
384 else { <<
385 // when do sampling common factors are n <<
386 // so cross section is not real <<
387 G4double cs = G4UniformRand(); <<
388 <<
389 if (gammaEnergy >= (*(fParamHigh[Z]))[1] <<
390 // above K-shell binding energy <<
391 cs *= fCrossSection->GetElementData(Z) <<
392 } <<
393 else { <<
394 // below K-shell binding energy <<
395 cs *= fCrossSectionLE->GetElementData( <<
396 } <<
397 <<
398 for (G4int j = 0; j < (G4int)nn; ++j) { <<
399 shellIdx = (std::size_t)fCrossSection- <<
400 if (gammaEnergy > (*(fParamLow[Z]))[7 <<
401 cs -= fCrossSection->GetValueForComp <<
402 } 298 }
403 if (cs <= 0.0 || j + 1 == (G4int)nn) { << 299 }
404 break; << 300
>> 301 // Returns the normalized direction of the momentum
>> 302 G4ThreeVector photonDirection = aDynamicGamma->GetMomentumDirection();
>> 303
>> 304 // Select randomly one element in the current material
>> 305 //G4cout << "Select random atom Egamma(keV)= " << gammaEnergy/keV << G4endl;
>> 306 const G4Element* elm = SelectRandomAtom(material, theGamma, gammaEnergy);
>> 307 G4int Z = elm->GetZasInt();
>> 308
>> 309 // Select the ionised shell in the current atom according to shell
>> 310 // cross sections
>> 311 // G4cout << "Select random shell Z= " << Z << G4endl;
>> 312
>> 313 if(Z >= maxZ) { Z = maxZ-1; }
>> 314
>> 315 // element was not initialised gamma should be absorbed
>> 316 if(!fCrossSection[Z]) {
>> 317 fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy);
>> 318 return;
>> 319 }
>> 320
>> 321 // SAMPLING OF THE SHELL INDEX
>> 322 size_t shellIdx = 0;
>> 323 size_t nn = fNShellsUsed[Z];
>> 324 if(nn > 1)
>> 325 {
>> 326 if(gammaEnergy >= (*(fParamHigh[Z]))[0])
>> 327 {
>> 328 G4double x1 = 1.0/gammaEnergy;
>> 329 G4double x2 = x1*x1;
>> 330 G4double x3 = x2*x1;
>> 331 G4double x4 = x3*x1;
>> 332 G4double x5 = x4*x1;
>> 333 G4int idx = nn*7 - 5;
>> 334 // when do sampling common factors are not taken into account
>> 335 // so cross section is not real
>> 336
>> 337 G4double rand=G4UniformRand();
>> 338 G4double cs0 = rand*( (*(fParamHigh[Z]))[idx]
>> 339 + x1*(*(fParamHigh[Z]))[idx+1]
>> 340 + x2*(*(fParamHigh[Z]))[idx+2]
>> 341 + x3*(*(fParamHigh[Z]))[idx+3]
>> 342 + x4*(*(fParamHigh[Z]))[idx+4]
>> 343 + x5*(*(fParamHigh[Z]))[idx+5]);
>> 344
>> 345 for(shellIdx=0; shellIdx<nn; ++shellIdx)
>> 346 {
>> 347 idx = shellIdx*7 + 2;
>> 348 if(gammaEnergy > (*(fParamHigh[Z]))[idx-1])
>> 349 {
>> 350 G4double cs =
>> 351 (*(fParamHigh[Z]))[idx]
>> 352 + x1*(*(fParamHigh[Z]))[idx+1]
>> 353 + x2*(*(fParamHigh[Z]))[idx+2]
>> 354 + x3*(*(fParamHigh[Z]))[idx+3]
>> 355 + x4*(*(fParamHigh[Z]))[idx+4]
>> 356 + x5*(*(fParamHigh[Z]))[idx+5];
>> 357
>> 358 if(cs >= cs0) { break; }
>> 359 }
>> 360 }
>> 361 if(shellIdx >= nn) { shellIdx = nn-1; }
>> 362 }
>> 363 else if(gammaEnergy >= (*(fParamLow[Z]))[0])
>> 364 {
>> 365 G4double x1 = 1.0/gammaEnergy;
>> 366 G4double x2 = x1*x1;
>> 367 G4double x3 = x2*x1;
>> 368 G4double x4 = x3*x1;
>> 369 G4double x5 = x4*x1;
>> 370 G4int idx = nn*7 - 5;
>> 371 // when do sampling common factors are not taken into account
>> 372 // so cross section is not real
>> 373 G4double cs0 = G4UniformRand()*((*(fParamLow[Z]))[idx]
>> 374 + x1*(*(fParamLow[Z]))[idx+1]
>> 375 + x2*(*(fParamLow[Z]))[idx+2]
>> 376 + x3*(*(fParamLow[Z]))[idx+3]
>> 377 + x4*(*(fParamLow[Z]))[idx+4]
>> 378 + x5*(*(fParamLow[Z]))[idx+5]);
>> 379 for(shellIdx=0; shellIdx<nn; ++shellIdx)
>> 380 {
>> 381 idx = shellIdx*7 + 2;
>> 382 if(gammaEnergy > (*(fParamLow[Z]))[idx-1])
>> 383 {
>> 384 G4double cs = (*(fParamLow[Z]))[idx] + x1*(*(fParamLow[Z]))[idx+1]
>> 385 + x2*(*(fParamLow[Z]))[idx+2] + x3*(*(fParamLow[Z]))[idx+3]
>> 386 + x4*(*(fParamLow[Z]))[idx+4]+ x5*(*(fParamLow[Z]))[idx+5];
>> 387 if(cs >= cs0) { break; }
>> 388 }
>> 389 }
>> 390 if(shellIdx >= nn) {shellIdx = nn-1;}
>> 391 }
>> 392 else
>> 393 {
>> 394 // when do sampling common factors are not taken into account
>> 395 // so cross section is not real
>> 396 G4double cs = G4UniformRand();
>> 397
>> 398 if(gammaEnergy >= (*(fParamHigh[Z]))[1]) {
>> 399 //above K-shell binding energy
>> 400 cs*= (fCrossSection[Z])->Value(gammaEnergy);
>> 401 }
>> 402 else
>> 403 {
>> 404 //below K-shell binding energy
>> 405 cs *= (fCrossSectionLE[Z])->Value(gammaEnergy);
>> 406 }
>> 407
>> 408 for(size_t j=0; j<nn; ++j) {
>> 409
>> 410 shellIdx = (size_t)fShellCrossSection->GetComponentID(Z, j);
>> 411 if(gammaEnergy > (*(fParamLow[Z]))[7*shellIdx+1]) {
>> 412 cs -= fShellCrossSection->GetValueForComponent(Z, j, gammaEnergy);
>> 413 }
>> 414 if(cs <= 0.0 || j+1 == nn) {break;}
>> 415 }
405 } 416 }
406 } <<
407 } 417 }
408 } << 418 // END: SAMPLING OF THE SHELL
409 // END: SAMPLING OF THE SHELL << 419
410 << 420 G4double bindingEnergy = (*(fParamHigh[Z]))[shellIdx*7 + 1];
411 G4double bindingEnergy = (*(fParamHigh[Z]))[ << 421 //G4cout << "Z= " << Z << " shellIdx= " << shellIdx
412 const G4AtomicShell* shell = nullptr; << 422 // << " nShells= " << fNShells[Z]
413 << 423 // << " Ebind(keV)= " << bindingEnergy/keV
414 // no de-excitation from the last shell << 424 // << " Egamma(keV)= " << gammaEnergy/keV << G4endl;
415 if (fDeexcitationActive && shellIdx + 1 < nn << 425
416 auto as = G4AtomicShellEnumerator(shellIdx << 426 const G4AtomicShell* shell = 0;
417 shell = fAtomDeexcitation->GetAtomicShell( << 427
418 } << 428 // no de-excitation from the last shell
419 << 429 if(fDeexcitationActive && shellIdx + 1 < nn) {
420 // If binding energy of the selected shell i << 430 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIdx);
421 // do not generate secondaries << 431 shell = fAtomDeexcitation->GetAtomicShell(Z, as);
422 if (gammaEnergy < bindingEnergy) { << 432 }
423 fParticleChange->ProposeLocalEnergyDeposit << 433
424 return; << 434 // If binding energy of the selected shell is larger than photon energy
425 } << 435 // do not generate secondaries
426 << 436 if(gammaEnergy < bindingEnergy) {
427 // Primary outcoming electron << 437 fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy);
428 G4double eKineticEnergy = gammaEnergy - bind << 438 return;
429 G4double edep = bindingEnergy; << 439 }
430 << 440
431 // Calculate direction of the photoelectron << 441 // Primary outcoming electron
432 G4ThreeVector electronDirection = GetAngular << 442 G4double eKineticEnergy = gammaEnergy - bindingEnergy;
433 aDynamicGamma, eKineticEnergy, (G4int)shel << 443 G4double edep = bindingEnergy;
434 << 444
435 // The electron is created << 445 // Calculate direction of the photoelectron
436 auto electron = << 446 G4ThreeVector electronDirection =
437 new G4DynamicParticle(theElectron, electro << 447 GetAngularDistribution()->SampleDirection(aDynamicGamma,
438 fvect->push_back(electron); << 448 eKineticEnergy,
439 << 449 shellIdx,
440 // Sample deexcitation << 450 couple->GetMaterial());
441 if (nullptr != shell) { << 451
442 G4int index = couple->GetIndex(); << 452 // The electron is created
443 if (fAtomDeexcitation->CheckDeexcitationAc << 453 G4DynamicParticle* electron = new G4DynamicParticle (theElectron,
444 std::size_t nbefore = fvect->size(); << 454 electronDirection,
445 << 455 eKineticEnergy);
446 fAtomDeexcitation->GenerateParticles(fve << 456 fvect->push_back(electron);
447 std::size_t nafter = fvect->size(); << 457
448 if (nafter > nbefore) { << 458 // Sample deexcitation
449 G4double esec = 0.0; << 459 if(shell) {
450 for (std::size_t j = nbefore; j < naft << 460 G4int index = couple->GetIndex();
451 G4double e = ((*fvect)[j])->GetKinet << 461 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
452 if (esec + e > edep) { << 462 G4int nbefore = fvect->size();
453 // correct energy in order to have << 463
454 e = edep - esec; << 464 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
455 ((*fvect)[j])->SetKineticEnergy(e) << 465 G4int nafter = fvect->size();
456 esec += e; << 466 if(nafter > nbefore) {
457 // delete the rest of secondaries << 467 G4double esec = 0.0;
458 for (std::size_t jj = nafter - 1; << 468 for (G4int j=nbefore; j<nafter; ++j) {
459 delete (*fvect)[jj]; << 469
460 fvect->pop_back(); << 470 G4double e = ((*fvect)[j])->GetKineticEnergy();
>> 471 if(esec + e > edep) {
>> 472 // correct energy in order to have energy balance
>> 473 e = edep - esec;
>> 474 ((*fvect)[j])->SetKineticEnergy(e);
>> 475 esec += e;
>> 476 // delete the rest of secondaries (should not happens)
>> 477 for (G4int jj=nafter-1; jj>j; --jj) {
>> 478 delete (*fvect)[jj];
>> 479 fvect->pop_back();
>> 480 }
>> 481 break;
>> 482 }
>> 483 esec += e;
>> 484 }
>> 485 edep -= esec;
461 } 486 }
462 break; <<
463 } <<
464 esec += e; <<
465 } 487 }
466 edep -= esec; <<
467 } <<
468 } 488 }
469 } << 489 // energy balance - excitation energy left
470 // energy balance - excitation energy left << 490 if(edep > 0.0) {
471 if (edep > 0.0) { << 491 fParticleChange->ProposeLocalEnergyDeposit(edep);
472 fParticleChange->ProposeLocalEnergyDeposit << 492 }
473 } <<
474 } 493 }
475 494
476 //....oooOO0OOooo........oooOO0OOooo........oo 495 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
477 496
478 const G4String& G4LivermorePhotoElectricModel: 497 const G4String& G4LivermorePhotoElectricModel::FindDirectoryPath()
479 { 498 {
480 // no check in this method - environment var << 499 // check environment variable
481 if (fDataDirectory.empty()) { << 500 // build the complete string identifying the file with the data set
482 auto param = G4EmParameters::Instance(); << 501 if(fDataDirectory.empty()) {
483 std::ostringstream ost; << 502 const char* path = std::getenv("G4LEDATA");
484 if (param->LivermoreDataDir() == "livermor << 503 if (path) {
485 ost << param->GetDirLEDATA() << "/liverm << 504 std::ostringstream ost;
486 } << 505 ost << path << "/livermore/phot_epics2014/";
487 else { << 506 fDataDirectory = ost.str();
488 ost << param->GetDirLEDATA() << "/epics2 << 507 } else {
>> 508 G4Exception("G4SeltzerBergerModel::FindDirectoryPath()","em0006",
>> 509 FatalException,
>> 510 "Environment variable G4LEDATA not defined");
489 } 511 }
490 fDataDirectory = ost.str(); <<
491 } 512 }
492 return fDataDirectory; 513 return fDataDirectory;
493 } 514 }
494 515
495 //....oooOO0OOooo........oooOO0OOooo........oo 516 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
496 517
497 void G4LivermorePhotoElectricModel::ReadData(G 518 void G4LivermorePhotoElectricModel::ReadData(G4int Z)
498 { 519 {
499 if (Z <= 0 || Z >= ZMAXPE) { << 520 if (verboseLevel > 1)
500 G4cout << "G4LivermorePhotoElectricModel:: << 521 {
501 << Z << " is out of range - request ignor << 522 G4cout << "Calling ReadData() of G4LivermorePhotoElectricModel"
502 return; << 523 << G4endl;
503 } << 524 }
504 if (verboseLevel > 1) { << 525
505 G4cout << "G4LivermorePhotoElectricModel:: << 526 if(fCrossSection[Z]) { return; }
506 } << 527
507 << 528 // spline for photoeffect total x-section above K-shell
508 if (fCrossSection->GetElementData(Z) != null << 529 // but below the parameterized ones
509 return; << 530 fCrossSection[Z] = new G4LPhysicsFreeVector();
510 } << 531 fCrossSection[Z]->SetSpline(true);
511 << 532 std::ostringstream ost;
512 // spline for photoeffect total x-section ab << 533 ost << FindDirectoryPath() << "pe-cs-" << Z <<".dat";
513 // but below the parameterized ones << 534 std::ifstream fin(ost.str().c_str());
514 << 535 if( !fin.is_open()) {
515 G4bool spline = (G4EmParameters::Instance()- << 536 G4ExceptionDescription ed;
516 G4int number = G4EmParameters::Instance()->N << 537 ed << "G4LivermorePhotoElectricModel data file <" << ost.str().c_str()
517 auto pv = new G4PhysicsFreeVector(spline); << 538 << "> is not opened!" << G4endl;
518 << 539 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
519 // fDataDirectory will be defined after thes << 540 "em0003",FatalException,
520 std::ostringstream ost; << 541 ed,"G4LEDATA version should be G4EMLOW7.2 or later.");
521 ost << FindDirectoryPath() << "pe-cs-" << Z << 542 return;
522 std::ifstream fin(ost.str().c_str()); << 543 }
523 if (!fin.is_open()) { << 544 if(verboseLevel > 3) {
524 G4ExceptionDescription ed; << 545 G4cout << "File " << ost.str().c_str()
525 ed << "G4LivermorePhotoElectricModel data << 546 << " is opened by G4LivermorePhotoElectricModel" << G4endl;
526 << ost.str().c_str() << "> is not opene << 547 }
527 G4Exception("G4LivermorePhotoElectricModel << 548 fCrossSection[Z]->Retrieve(fin, true);
528 "G4LEDATA version should be G4 << 549 fCrossSection[Z]->ScaleVector(MeV, barn);
529 return; << 550 fin.close();
530 } << 551
531 if (verboseLevel > 3) { << 552 // read high-energy fit parameters
532 G4cout << "File " << ost.str().c_str() << << 553 fParamHigh[Z] = new std::vector<G4double>;
533 << G4endl; << 554 G4int n1 = 0;
534 } << 555 G4int n2 = 0;
535 pv->Retrieve(fin, true); << 556 G4double x;
536 pv->ScaleVector(MeV, barn); << 557 std::ostringstream ost1;
537 pv->FillSecondDerivatives(); << 558 ost1 << fDataDirectory << "pe-high-" << Z <<".dat";
538 pv->EnableLogBinSearch(number); << 559 std::ifstream fin1(ost1.str().c_str());
539 fCrossSection->InitialiseForElement(Z, pv); << 560 if( !fin1.is_open()) {
540 fin.close(); << 561 G4ExceptionDescription ed;
541 << 562 ed << "G4LivermorePhotoElectricModel data file <" << ost1.str().c_str()
542 // read high-energy fit parameters << 563 << "> is not opened!" << G4endl;
543 fParamHigh[Z] = new std::vector<G4double>; << 564 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
544 G4int n1 = 0; << 565 "em0003",FatalException,
545 G4int n2 = 0; << 566 ed,"G4LEDATA version should be G4EMLOW7.2 or later.");
546 G4double x; << 567 return;
547 std::ostringstream ost1; << 568 }
548 ost1 << fDataDirectory << "pe-high-" << Z << << 569 if(verboseLevel > 3) {
549 std::ifstream fin1(ost1.str().c_str()); << 570 G4cout << "File " << ost1.str().c_str()
550 if (!fin1.is_open()) { << 571 << " is opened by G4LivermorePhotoElectricModel" << G4endl;
551 G4ExceptionDescription ed; << 572 }
552 ed << "G4LivermorePhotoElectricModel data << 573 fin1 >> n1;
553 << ost1.str().c_str() << "> is not open << 574 if(fin1.fail()) { return; }
554 G4Exception("G4LivermorePhotoElectricModel << 575 if(0 > n1 || n1 >= INT_MAX) { n1 = 0; }
555 "G4LEDATA version should be G4 << 576
556 return; << 577 fin1 >> n2;
557 } << 578 if(fin1.fail()) { return; }
558 if (verboseLevel > 3) { << 579 if(0 > n2 || n2 >= INT_MAX) { n2 = 0; }
559 G4cout << "File " << ost1.str().c_str() << << 580
560 << G4endl; << 581 fin1 >> x;
561 } << 582 if(fin1.fail()) { return; }
562 fin1 >> n1; << 583
563 if (fin1.fail()) { << 584 fNShells[Z] = n1;
564 return; << 585 fParamHigh[Z]->reserve(7*n1+1);
565 } << 586 fParamHigh[Z]->push_back(x*MeV);
566 if (0 > n1 || n1 >= INT_MAX) { << 587 for(G4int i=0; i<n1; ++i) {
567 n1 = 0; << 588 for(G4int j=0; j<7; ++j) {
568 } << 589 fin1 >> x;
569 << 590 if(0 == j) { x *= MeV; }
570 fin1 >> n2; << 591 else { x *= barn; }
571 if (fin1.fail()) { << 592 fParamHigh[Z]->push_back(x);
572 return; << 593 }
573 } <<
574 if (0 > n2 || n2 >= INT_MAX) { <<
575 n2 = 0; <<
576 } <<
577 <<
578 fin1 >> x; <<
579 if (fin1.fail()) { <<
580 return; <<
581 } <<
582 <<
583 fNShells[Z] = n1; <<
584 fParamHigh[Z]->reserve(7 * n1 + 1); <<
585 fParamHigh[Z]->push_back(x * MeV); <<
586 for (G4int i = 0; i < n1; ++i) { <<
587 for (G4int j = 0; j < 7; ++j) { <<
588 fin1 >> x; <<
589 if (0 == j) { <<
590 x *= MeV; <<
591 } <<
592 else { <<
593 x *= barn; <<
594 } <<
595 fParamHigh[Z]->push_back(x); <<
596 } <<
597 } <<
598 fin1.close(); <<
599 <<
600 // read low-energy fit parameters <<
601 fParamLow[Z] = new std::vector<G4double>; <<
602 G4int n1_low = 0; <<
603 G4int n2_low = 0; <<
604 G4double x_low; <<
605 std::ostringstream ost1_low; <<
606 ost1_low << fDataDirectory << "pe-low-" << Z <<
607 std::ifstream fin1_low(ost1_low.str().c_str( <<
608 if (!fin1_low.is_open()) { <<
609 G4ExceptionDescription ed; <<
610 ed << "G4LivermorePhotoElectricModel data <<
611 << "> is not opened!" << G4endl; <<
612 G4Exception("G4LivermorePhotoElectricModel <<
613 "G4LEDATA version should be G4 <<
614 return; <<
615 } <<
616 if (verboseLevel > 3) { <<
617 G4cout << "File " << ost1_low.str().c_str( <<
618 << G4endl; <<
619 } <<
620 fin1_low >> n1_low; <<
621 if (fin1_low.fail()) { <<
622 return; <<
623 } <<
624 if (0 > n1_low || n1_low >= INT_MAX) { <<
625 n1_low = 0; <<
626 } <<
627 <<
628 fin1_low >> n2_low; <<
629 if (fin1_low.fail()) { <<
630 return; <<
631 } <<
632 if (0 > n2_low || n2_low >= INT_MAX) { <<
633 n2_low = 0; <<
634 } <<
635 <<
636 fin1_low >> x_low; <<
637 if (fin1_low.fail()) { <<
638 return; <<
639 } <<
640 <<
641 fNShells[Z] = n1_low; <<
642 fParamLow[Z]->reserve(7 * n1_low + 1); <<
643 fParamLow[Z]->push_back(x_low * MeV); <<
644 for (G4int i = 0; i < n1_low; ++i) { <<
645 for (G4int j = 0; j < 7; ++j) { <<
646 fin1_low >> x_low; <<
647 if (0 == j) { <<
648 x_low *= MeV; <<
649 } <<
650 else { <<
651 x_low *= barn; <<
652 } <<
653 fParamLow[Z]->push_back(x_low); <<
654 } <<
655 } <<
656 fin1_low.close(); <<
657 <<
658 // there is a possibility to use only main s <<
659 if (nShellLimit < n2) { <<
660 n2 = nShellLimit; <<
661 } <<
662 fCrossSection->InitialiseForComponent(Z, n2) <<
663 fNShellsUsed[Z] = n2; <<
664 <<
665 if (1 < n2) { <<
666 std::ostringstream ost2; <<
667 ost2 << fDataDirectory << "pe-ss-cs-" << Z <<
668 std::ifstream fin2(ost2.str().c_str()); <<
669 if (!fin2.is_open()) { <<
670 G4ExceptionDescription ed; <<
671 ed << "G4LivermorePhotoElectricModel dat <<
672 << G4endl; <<
673 G4Exception("G4LivermorePhotoElectricMod <<
674 "G4LEDATA version should be <<
675 return; <<
676 } <<
677 if (verboseLevel > 3) { <<
678 G4cout << "File " << ost2.str().c_str() <<
679 << G4endl; <<
680 } 594 }
681 << 595 fin1.close();
682 G4int n3, n4; << 596
683 G4double y; << 597 // read low-energy fit parameters
684 for (G4int i = 0; i < n2; ++i) { << 598 fParamLow[Z] = new std::vector<G4double>;
685 fin2 >> x >> y >> n3 >> n4; << 599 G4int n1_low = 0;
686 auto v = new G4PhysicsFreeVector(n3, x, << 600 G4int n2_low = 0;
687 for (G4int j = 0; j < n3; ++j) { << 601 G4double x_low;
688 fin2 >> x >> y; << 602 std::ostringstream ost1_low;
689 v->PutValues(j, x * CLHEP::MeV, y * CL << 603 ost1_low << fDataDirectory << "pe-low-" << Z <<".dat";
690 } << 604 std::ifstream fin1_low(ost1_low.str().c_str());
691 v->EnableLogBinSearch(number); << 605 if( !fin1_low.is_open()) {
692 fCrossSection->AddComponent(Z, n4, v); << 606 G4ExceptionDescription ed;
>> 607 ed << "G4LivermorePhotoElectricModel data file <" << ost1_low.str().c_str()
>> 608 << "> is not opened!" << G4endl;
>> 609 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
>> 610 "em0003",FatalException,
>> 611 ed,"G4LEDATA version should be G4EMLOW7.2 or later.");
>> 612 return;
>> 613 }
>> 614 if(verboseLevel > 3) {
>> 615 G4cout << "File " << ost1_low.str().c_str()
>> 616 << " is opened by G4LivermorePhotoElectricModel" << G4endl;
>> 617 }
>> 618 fin1_low >> n1_low;
>> 619 if(fin1_low.fail()) { return; }
>> 620 if(0 > n1_low || n1_low >= INT_MAX) { n1_low = 0; }
>> 621
>> 622 fin1_low >> n2_low;
>> 623 if(fin1_low.fail()) { return; }
>> 624 if(0 > n2_low || n2_low >= INT_MAX) { n2_low = 0; }
>> 625
>> 626 fin1_low >> x_low;
>> 627 if(fin1_low.fail()) { return; }
>> 628
>> 629 fNShells[Z] = n1_low;
>> 630 fParamLow[Z]->reserve(7*n1_low+1);
>> 631 fParamLow[Z]->push_back(x_low*MeV);
>> 632 for(G4int i=0; i<n1_low; ++i) {
>> 633 for(G4int j=0; j<7; ++j) {
>> 634 fin1_low >> x_low;
>> 635 if(0 == j) { x_low *= MeV; }
>> 636 else { x_low *= barn; }
>> 637 fParamLow[Z]->push_back(x_low);
>> 638 }
693 } 639 }
694 fin2.close(); << 640 fin1_low.close();
695 } << 641
696 << 642 // there is a possibility to use only main shells
697 // no spline for photoeffect total x-section << 643 if(nShellLimit < n2) { n2 = nShellLimit; }
698 if (1 < fNShells[Z]) { << 644 fShellCrossSection->InitialiseForComponent(Z, n2);//number of shells
699 auto pv1 = new G4PhysicsFreeVector(false); << 645 fNShellsUsed[Z] = n2;
700 std::ostringstream ost3; << 646
701 ost3 << fDataDirectory << "pe-le-cs-" << Z << 647 if(1 < n2) {
702 std::ifstream fin3(ost3.str().c_str()); << 648 std::ostringstream ost2;
703 if (!fin3.is_open()) { << 649 ost2 << fDataDirectory << "pe-ss-cs-" << Z <<".dat";
704 G4ExceptionDescription ed; << 650 std::ifstream fin2(ost2.str().c_str());
705 ed << "G4LivermorePhotoElectricModel dat << 651 if( !fin2.is_open()) {
706 << G4endl; << 652 G4ExceptionDescription ed;
707 G4Exception("G4LivermorePhotoElectricMod << 653 ed << "G4LivermorePhotoElectricModel data file <" << ost2.str().c_str()
708 "G4LEDATA version should be << 654 << "> is not opened!" << G4endl;
709 return; << 655 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
>> 656 "em0003",FatalException,
>> 657 ed,"G4LEDATA version should be G4EMLOW7.2 or later.");
>> 658 return;
>> 659 }
>> 660 if(verboseLevel > 3) {
>> 661 G4cout << "File " << ost2.str().c_str()
>> 662 << " is opened by G4LivermorePhotoElectricModel" << G4endl;
>> 663 }
>> 664
>> 665 G4int n3, n4;
>> 666 G4double y;
>> 667 for(G4int i=0; i<n2; ++i) {
>> 668 fin2 >> x >> y >> n3 >> n4;
>> 669 G4LPhysicsFreeVector* v = new G4LPhysicsFreeVector(n3, x, y);
>> 670 for(G4int j=0; j<n3; ++j) {
>> 671 fin2 >> x >> y;
>> 672 v->PutValues(j, x*MeV, y*barn);
>> 673 }
>> 674 fShellCrossSection->AddComponent(Z, n4, v);
>> 675 }
>> 676 fin2.close();
710 } 677 }
711 if (verboseLevel > 3) { << 678
712 G4cout << "File " << ost3.str().c_str() << 679 // no spline for photoeffect total x-section below K-shell
713 << G4endl; << 680 if(1 < fNShells[Z]) {
>> 681 fCrossSectionLE[Z] = new G4LPhysicsFreeVector();
>> 682 std::ostringstream ost3;
>> 683 ost3 << fDataDirectory << "pe-le-cs-" << Z <<".dat";
>> 684 std::ifstream fin3(ost3.str().c_str());
>> 685 if( !fin3.is_open()) {
>> 686 G4ExceptionDescription ed;
>> 687 ed << "G4LivermorePhotoElectricModel data file <" << ost3.str().c_str()
>> 688 << "> is not opened!" << G4endl;
>> 689 G4Exception("G4LivermorePhotoElectricModel::ReadData()",
>> 690 "em0003",FatalException,
>> 691 ed,"G4LEDATA version should be G4EMLOW7.2 or later.");
>> 692 return;
>> 693 }
>> 694 if(verboseLevel > 3) {
>> 695 G4cout << "File " << ost3.str().c_str()
>> 696 << " is opened by G4LivermorePhotoElectricModel" << G4endl;
>> 697 }
>> 698 fCrossSectionLE[Z]->Retrieve(fin3, true);
>> 699 fCrossSectionLE[Z]->ScaleVector(MeV, barn);
>> 700 fin3.close();
714 } 701 }
715 pv1->Retrieve(fin3, true); <<
716 pv1->ScaleVector(CLHEP::MeV, CLHEP::barn); <<
717 pv1->EnableLogBinSearch(number); <<
718 fCrossSectionLE->InitialiseForElement(Z, p <<
719 fin3.close(); <<
720 } <<
721 } 702 }
722 703
723 //....oooOO0OOooo........oooOO0OOooo........oo 704 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
724 705
725 G4double G4LivermorePhotoElectricModel::GetBin 706 G4double G4LivermorePhotoElectricModel::GetBindingEnergy(G4int Z, G4int shell)
726 { 707 {
727 if (Z < 1 || Z >= ZMAXPE) { << 708 if(Z < 1 || Z >= maxZ) { return -1;} //If Z is out of the supported return 0
728 return -1; << 709 //If necessary load data for Z
729 } // If Z is out of the supported return 0 << 710 InitialiseForElement(theGamma, Z);
730 << 711 if(!fCrossSection[Z] || shell < 0 || shell >= fNShellsUsed[Z]) { return -1; }
731 InitialiseOnFly(Z); << 712
732 if (fCrossSection->GetElementData(Z) == null << 713 if(Z>2)
733 shell < 0 || shell >= fNShellsUsed[Z]) { << 714 return fShellCrossSection->GetComponentDataByIndex(Z, shell)->Energy(0);
734 return -1; << 715 else
735 } << 716 return fCrossSection[Z]->Energy(0);
736 <<
737 if (Z > 2) { <<
738 return fCrossSection->GetComponentDataByIn <<
739 } <<
740 else { <<
741 return fCrossSection->GetElementData(Z)->E <<
742 } <<
743 } <<
744 <<
745 //....oooOO0OOooo........oooOO0OOooo........oo <<
746 <<
747 void <<
748 G4LivermorePhotoElectricModel::InitialiseForEl <<
749 { <<
750 if (fCrossSection == nullptr) { <<
751 fCrossSection = new G4ElementData(ZMAXPE); <<
752 fCrossSection->SetName("PhotoEffXS"); <<
753 fCrossSectionLE = new G4ElementData(ZMAXPE <<
754 fCrossSectionLE->SetName("PhotoEffLowXS"); <<
755 } <<
756 ReadData(Z); <<
757 } 717 }
758 718
759 //....oooOO0OOooo........oooOO0OOooo........oo 719 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
760 720
761 void G4LivermorePhotoElectricModel::Initialise << 721 void G4LivermorePhotoElectricModel::InitialiseForElement(
>> 722 const G4ParticleDefinition*, G4int Z)
762 { 723 {
763 if (fCrossSection->GetElementData(Z) == null << 724 if (!fCrossSection[Z]) {
764 G4AutoLock l(&livPhotoeffMutex); << 725 #ifdef G4MULTITHREADED
765 if (fCrossSection->GetElementData(Z) == nu << 726 G4MUTEXLOCK(&livPhotoeffMutex);
>> 727 if (!fCrossSection[Z]) {
>> 728 #endif
766 ReadData(Z); 729 ReadData(Z);
>> 730 #ifdef G4MULTITHREADED
767 } 731 }
768 l.unlock(); << 732 G4MUTEXUNLOCK(&livPhotoeffMutex);
769 } << 733 #endif
>> 734 }
770 } 735 }
771 736
772 //....oooOO0OOooo........oooOO0OOooo........oo 737 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
773 738