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