| 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,v 1.1 2008/10/30 14:16:35 sincerti Exp $
>> 27 // GEANT4 tag $Name: geant4-09-02-patch-01 $
26 // 28 //
27 // Author: Sebastien Incerti <<
28 // 30 October 2008 <<
29 // on base of G4LowEnergyPhotoElectric <<
30 // <<
31 // 22 Oct 2012 A & V Ivanchenko Migration da <<
32 // 1 June 2017 M Bandieramonte: New model ba <<
33 // evaluated data - parameteriza <<
34 29
35 #include "G4LivermorePhotoElectricModel.hh" 30 #include "G4LivermorePhotoElectricModel.hh"
36 31
37 #include "G4AtomicShell.hh" <<
38 #include "G4AutoLock.hh" <<
39 #include "G4CrossSectionHandler.hh" <<
40 #include "G4Electron.hh" <<
41 #include "G4Gamma.hh" <<
42 #include "G4LossTableManager.hh" <<
43 #include "G4ParticleChangeForGamma.hh" <<
44 #include "G4PhysicsFreeVector.hh" <<
45 #include "G4SauterGavrilaAngularDistribution.h <<
46 #include "G4SystemOfUnits.hh" <<
47 #include "G4VAtomDeexcitation.hh" <<
48 #include "G4EmParameters.hh" <<
49 #include <thread> <<
50 <<
51 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
52 33
53 G4ElementData* G4LivermorePhotoElectricModel:: << 34 using namespace std;
54 G4ElementData* G4LivermorePhotoElectricModel:: <<
55 std::vector<G4double>* G4LivermorePhotoElectri <<
56 std::vector<G4double>* G4LivermorePhotoElectri <<
57 G4int G4LivermorePhotoElectricModel::fNShells[ <<
58 G4int G4LivermorePhotoElectricModel::fNShellsU <<
59 G4Material* G4LivermorePhotoElectricModel::fWa <<
60 G4double G4LivermorePhotoElectricModel::fWater <<
61 G4String G4LivermorePhotoElectricModel::fDataD <<
62 <<
63 static std::once_flag applyOnce; <<
64 <<
65 namespace <<
66 { <<
67 G4Mutex livPhotoeffMutex = G4MUTEX_INITIALIZ <<
68 } <<
69 35
70 //....oooOO0OOooo........oooOO0OOooo........oo 36 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 37
72 G4LivermorePhotoElectricModel::G4LivermorePhot << 38 G4LivermorePhotoElectricModel::G4LivermorePhotoElectricModel(const G4ParticleDefinition*,
73 { << 39 const G4String& nam)
74 verboseLevel = 0; << 40 :G4VEmModel(nam),isInitialised(false)
>> 41 {
>> 42 lowEnergyLimit = 250 * eV; // SI - Could be 10 eV ?
>> 43 highEnergyLimit = 100 * GeV;
>> 44 SetLowEnergyLimit(lowEnergyLimit);
>> 45 SetHighEnergyLimit(highEnergyLimit);
>> 46
>> 47 G4cout << "Livermore Compton is constructed " << G4endl
>> 48 << "Energy range: "
>> 49 << lowEnergyLimit / keV << " keV - "
>> 50 << highEnergyLimit / GeV << " GeV"
>> 51 << G4endl;
>> 52
>> 53 verboseLevel= 0;
75 // Verbosity scale: 54 // Verbosity scale:
76 // 0 = nothing << 55 // 0 = nothing
77 // 1 = warning for energy non-conservation << 56 // 1 = warning for energy non-conservation
78 // 2 = details of energy budget 57 // 2 = details of energy budget
79 // 3 = calculation of cross sections, file o 58 // 3 = calculation of cross sections, file openings, sampling of atoms
80 // 4 = entering in methods 59 // 4 = entering in methods
81 <<
82 theGamma = G4Gamma::Gamma(); <<
83 theElectron = G4Electron::Electron(); <<
84 <<
85 // default generator <<
86 SetAngularDistribution(new G4SauterGavrilaAn <<
87 <<
88 if (verboseLevel > 0) { <<
89 G4cout << "Livermore PhotoElectric is cons <<
90 << " nShellLimit= " << nShellLimit <<
91 } <<
92 <<
93 // Mark this model as "applicable" for atomi <<
94 SetDeexcitationFlag(true); <<
95 <<
96 // For water <<
97 fSandiaCof.resize(4, 0.0); <<
98 } 60 }
99 61
100 //....oooOO0OOooo........oooOO0OOooo........oo 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 63
102 G4LivermorePhotoElectricModel::~G4LivermorePho 64 G4LivermorePhotoElectricModel::~G4LivermorePhotoElectricModel()
103 { << 65 {
104 if (isInitializer) { << 66 delete meanFreePathTable;
105 for (G4int i = 0; i < ZMAXPE; ++i) { << 67 delete crossSectionHandler;
106 if (fParamHigh[i]) { << 68 delete shellCrossSectionHandler;
107 delete fParamHigh[i]; << 69 delete ElectronAngularGenerator;
108 fParamHigh[i] = nullptr; <<
109 } <<
110 if (fParamLow[i]) { <<
111 delete fParamLow[i]; <<
112 fParamLow[i] = nullptr; <<
113 } <<
114 } <<
115 } <<
116 } 70 }
117 71
118 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119 73
120 void G4LivermorePhotoElectricModel::Initialise << 74 void G4LivermorePhotoElectricModel::Initialise(const G4ParticleDefinition* particle,
121 const G4DataVector&) << 75 const G4DataVector& cuts)
122 { << 76 {
123 if (verboseLevel > 1) { << 77 if (verboseLevel > 3)
124 G4cout << "Calling G4LivermorePhotoElectri << 78 G4cout << "Calling G4LivermorePhotoElectricModel::Initialise()" << G4endl;
125 } << 79
126 << 80 InitialiseElementSelectors(particle,cuts);
127 // initialise static tables only once << 81
128 std::call_once(applyOnce, [this]() { isIniti << 82 // Energy limits
129 << 83
130 if (isInitializer) { << 84 if (LowEnergyLimit() < lowEnergyLimit)
131 G4AutoLock l(&livPhotoeffMutex); << 85 {
132 FindDirectoryPath(); << 86 G4cout << "G4LivermorePhotoElectricModel: low energy limit increased from " <<
133 if (fWater == nullptr) { << 87 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" <<
134 fWater = G4Material::GetMaterial("G4_WAT << 88 G4endl;
135 if (fWater == nullptr) { << 89 SetLowEnergyLimit(lowEnergyLimit);
136 fWater = G4Material::GetMaterial("Wate << 90 }
137 } << 91
138 if (fWater != nullptr) { << 92 if (HighEnergyLimit() > highEnergyLimit)
139 fWaterEnergyLimit = 13.6 * CLHEP::eV; << 93 {
140 } << 94 G4cout << "G4LivermorePhotoElectricModel: high energy limit decreased from " <<
141 } << 95 HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV"
142 << 96 << G4endl;
143 if (fCrossSection == nullptr) { << 97 SetHighEnergyLimit(highEnergyLimit);
144 fCrossSection = new G4ElementData(ZMAXPE << 98 }
145 fCrossSection->SetName("PhotoEffXS"); << 99
146 fCrossSectionLE = new G4ElementData(ZMAX << 100 // Read data tables for all materials
147 fCrossSectionLE->SetName("PhotoEffLowXS" << 101
148 } << 102 crossSectionHandler = new G4CrossSectionHandler();
149 << 103 crossSectionHandler->Clear();
150 const G4ElementTable* elemTable = G4Elemen << 104 G4String crossSectionFile = "phot/pe-cs-";
151 std::size_t numElems = (*elemTable).size() << 105 crossSectionHandler->LoadData(crossSectionFile);
152 for (std::size_t ie = 0; ie < numElems; ++ << 106
153 const G4Element* elem = (*elemTable)[ie] << 107 meanFreePathTable = 0;
154 G4int Z = elem->GetZasInt(); << 108 meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
155 if (Z < ZMAXPE) { << 109
156 if (fCrossSection->GetElementData(Z) == null << 110 shellCrossSectionHandler = new G4CrossSectionHandler();
157 ReadData(Z); << 111 shellCrossSectionHandler->Clear();
158 } << 112 G4String shellCrossSectionFile = "phot/pe-ss-cs-";
159 } << 113 shellCrossSectionHandler->LoadShellData(shellCrossSectionFile);
160 } << 114
161 l.unlock(); << 115 // SI - Buggy default ?
162 } << 116 //generatorName = "geant4.6.2";
>> 117 //ElectronAngularGenerator = new G4PhotoElectricAngularGeneratorSimple("GEANTSimpleGenerator"); // default generator
>> 118
>> 119 //
>> 120
>> 121 if (verboseLevel > 2)
>> 122 G4cout << "Loaded cross section files for Livermore PhotoElectric model" << G4endl;
>> 123
>> 124 G4cout << "Livermore PhotoElectric model is initialized " << G4endl
>> 125 << "Energy range: "
>> 126 << LowEnergyLimit() / keV << " keV - "
>> 127 << HighEnergyLimit() / GeV << " GeV"
>> 128 << G4endl;
163 129
164 if (verboseLevel > 1) { << 130 if(isInitialised) return;
165 G4cout << "Loaded cross section files for <<
166 } <<
167 if (nullptr == fParticleChange) { <<
168 fParticleChange = GetParticleChangeForGamm <<
169 fAtomDeexcitation = G4LossTableManager::In <<
170 } <<
171 131
172 fDeexcitationActive = false; << 132 if(pParticleChange)
173 if (nullptr != fAtomDeexcitation) { << 133 fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
174 fDeexcitationActive = fAtomDeexcitation->I << 134 else
175 } << 135 fParticleChange = new G4ParticleChangeForGamma();
176 136
177 if (verboseLevel > 1) { << 137 isInitialised = true;
178 G4cout << "LivermorePhotoElectric model is <<
179 } <<
180 } 138 }
181 139
182 //....oooOO0OOooo........oooOO0OOooo........oo 140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
183 141
184 G4double << 142 G4double G4LivermorePhotoElectricModel::ComputeCrossSectionPerAtom(
185 G4LivermorePhotoElectricModel::CrossSectionPer << 143 const G4ParticleDefinition*,
186 << 144 G4double GammaEnergy,
187 << 145 G4double Z, G4double,
188 G4double, G4double) << 146 G4double, G4double)
189 { 147 {
190 fCurrSection = 0.0; << 148 if (verboseLevel > 3)
191 if (fWater && (material == fWater || materia << 149 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePhotoElectricModel" << G4endl;
192 if (energy <= fWaterEnergyLimit) { <<
193 fWater->GetSandiaTable()->GetSandiaCofWa <<
194 <<
195 G4double energy2 = energy * energy; <<
196 G4double energy3 = energy * energy2; <<
197 G4double energy4 = energy2 * energy2; <<
198 <<
199 fCurrSection = material->GetDensity() <<
200 * (fSandiaCof[0] / energy <<
201 fSandiaCof[2] / energy3 + fSandiaCof[3] <<
202 } <<
203 } <<
204 if (0.0 == fCurrSection) { <<
205 fCurrSection = G4VEmModel::CrossSectionPer <<
206 } <<
207 return fCurrSection; <<
208 } <<
209 <<
210 //....oooOO0OOooo........oooOO0OOooo........oo <<
211 150
212 G4double << 151 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
213 G4LivermorePhotoElectricModel::ComputeCrossSec <<
214 <<
215 G4double ZZ, G4double, <<
216 G4double, G4double) <<
217 { <<
218 if (verboseLevel > 3) { <<
219 G4cout << "\n G4LivermorePhotoElectricMode <<
220 << " Z= " << ZZ << " R(keV)= " << <<
221 } <<
222 G4double cs = 0.0; <<
223 G4int Z = G4lrint(ZZ); <<
224 if (Z >= ZMAXPE || Z <= 0) { <<
225 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; 152 return cs;
276 } 153 }
277 154
278 //....oooOO0OOooo........oooOO0OOooo........oo 155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
279 156
280 void G4LivermorePhotoElectricModel::SampleSeco 157 void G4LivermorePhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
281 << 158 const G4MaterialCutsCouple* couple,
282 << 159 const G4DynamicParticle* aDynamicGamma,
283 << 160 G4double,
284 { << 161 G4double)
285 G4double gammaEnergy = aDynamicGamma->GetKin << 162 {
286 if (verboseLevel > 3) { << 163
287 G4cout << "G4LivermorePhotoElectricModel:: << 164 // Fluorescence generated according to:
288 << gammaEnergy / keV << G4endl; << 165 // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic
289 } << 166 // subshell ionisation by a particle or due to deexcitation or decay of radionuclides",
290 << 167 // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
291 // kill incident photon << 168
292 fParticleChange->ProposeTrackStatus(fStopAnd << 169 if (verboseLevel > 3)
293 fParticleChange->SetProposedKineticEnergy(0. << 170 G4cout << "Calling SampleSecondaries() of G4LivermorePhotoElectricModel" << G4endl;
294 << 171
295 // low-energy photo-effect in water - full a << 172 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
296 const G4Material* material = couple->GetMate << 173
297 if (fWater && (material == fWater || materia << 174 if (photonEnergy <= lowEnergyLimit)
298 if (gammaEnergy <= fWaterEnergyLimit) { << 175 {
299 fParticleChange->ProposeLocalEnergyDepos << 176 fParticleChange->ProposeTrackStatus(fStopAndKill);
300 return; << 177 fParticleChange->SetProposedKineticEnergy(0.);
>> 178 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
>> 179 // SI - IS THE FOLLOWING RETURN NECESSARY ?
>> 180 return ;
301 } 181 }
302 } << 182
303 << 183 G4ThreeVector photonDirection = aDynamicGamma->GetMomentumDirection(); // Returns the normalized direction of the momentum
304 // Returns the normalized direction of the m <<
305 G4ThreeVector photonDirection = aDynamicGamm <<
306 184
307 // Select randomly one element in the curren 185 // Select randomly one element in the current material
308 const G4Element* elm = SelectRandomAtom(mate << 186 G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
309 G4int Z = elm->GetZasInt(); <<
310 187
311 // Select the ionised shell in the current a << 188 // Select the ionised shell in the current atom according to shell cross sections
312 // cross sections << 189 size_t shellIndex = shellCrossSectionHandler->SelectRandomShell(Z,photonEnergy);
313 190
314 // If element was not initialised gamma shou << 191 // Retrieve the corresponding identifier and binding energy of the selected shell
315 if (Z >= ZMAXPE || fCrossSection->GetElement << 192 const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
316 fParticleChange->ProposeLocalEnergyDeposit << 193 const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
317 return; << 194 G4double bindingEnergy = shell->BindingEnergy();
318 } << 195 G4int shellId = shell->ShellId();
>> 196
>> 197 // Create lists of pointers to DynamicParticles (photons and electrons)
>> 198 // (Is the electron vector necessary? To be checked)
>> 199 std::vector<G4DynamicParticle*>* photonVector = 0;
>> 200 std::vector<G4DynamicParticle*> electronVector;
319 201
320 // SAMPLING OF THE SHELL INDEX << 202 G4double energyDeposit = 0.0;
321 std::size_t shellIdx = 0; << 203
322 std::size_t nn = fNShellsUsed[Z]; << 204 // Primary outcoming electron
323 if (nn > 1) { << 205 G4double eKineticEnergy = photonEnergy - bindingEnergy;
324 if (gammaEnergy >= (*(fParamHigh[Z]))[0]) << 206
325 G4double x1 = 1.0 / gammaEnergy; << 207 // There may be cases where the binding energy of the selected shell is > photon energy
326 G4double x2 = x1 * x1; << 208 // In such cases do not generate secondaries
327 G4double x3 = x2 * x1; << 209 if (eKineticEnergy > 0.)
328 G4double x4 = x3 * x1; << 210 {
329 G4double x5 = x4 * x1; << 211 // SI - Removed safety
330 std::size_t idx = nn * 7 - 5; << 212
331 // when do sampling common factors are n << 213 // Generate the electron only if with large enough range w.r.t. cuts and safety
332 // so cross section is not real << 214 //G4double safety = aStep.GetPostStepPoint()->GetSafety();
333 << 215
334 G4double rand = G4UniformRand(); << 216 //if (rangeTest->Escape(G4Electron::Electron(),couple,eKineticEnergy,safety))
335 G4double cs0 = rand * << 217 {
336 ((*(fParamHigh[Z]))[idx] + x1 * (*(fPa << 218
337 + x2 * (*(fParamHigh[Z]))[idx + 2] + x << 219 // Calculate direction of the photoelectron
338 + x4 * (*(fParamHigh[Z]))[idx + 4] + x << 220 G4ThreeVector gammaPolarization = aDynamicGamma->GetPolarization();
339 << 221 G4ThreeVector electronDirection = ElectronAngularGenerator->GetPhotoElectronDirection(photonDirection,eKineticEnergy,gammaPolarization,shellId);
340 for (shellIdx = 0; shellIdx < nn; ++shel << 222
341 idx = shellIdx * 7 + 2; << 223 // The electron is created ...
342 if (gammaEnergy > (*(fParamHigh[Z]))[i << 224 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
343 G4double cs = (*(fParamHigh[Z]))[idx << 225 electronDirection,
344 + x2 * (*(fParamHigh[Z]))[idx + 2] << 226 eKineticEnergy);
345 + x4 * (*(fParamHigh[Z]))[idx + 4] << 227 electronVector.push_back(electron);
346 << 228 }
347 if (cs >= cs0) { << 229 /*else
348 break; << 230 {
349 } << 231 energyDeposit += eKineticEnergy;
350 } << 232 }*/
351 } << 233 }
352 if (shellIdx >= nn) { << 234 else
353 shellIdx = nn - 1; << 235 {
354 } << 236 bindingEnergy = photonEnergy;
>> 237 }
>> 238
>> 239 G4int nElectrons = electronVector.size();
>> 240 size_t nTotPhotons = 0;
>> 241 G4int nPhotons=0;
>> 242 const G4ProductionCutsTable* theCoupleTable=
>> 243 G4ProductionCutsTable::GetProductionCutsTable();
>> 244
>> 245 size_t index = couple->GetIndex();
>> 246 G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
>> 247 cutg = std::min(cutForLowEnergySecondaryPhotons,cutg);
>> 248
>> 249 G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
>> 250 cute = std::min(cutForLowEnergySecondaryPhotons,cute);
>> 251
>> 252 G4DynamicParticle* aPhoton;
>> 253
>> 254 // Generation of fluorescence
>> 255 // Data in EADL are available only for Z > 5
>> 256 // Protection to avoid generating photons in the unphysical case of
>> 257 // shell binding energy > photon energy
>> 258 if (Z > 5 && (bindingEnergy > cutg || bindingEnergy > cute))
>> 259 {
>> 260 photonVector = deexcitationManager.GenerateParticles(Z,shellId);
>> 261 nTotPhotons = photonVector->size();
>> 262 for (size_t k=0; k<nTotPhotons; k++)
>> 263 {
>> 264 aPhoton = (*photonVector)[k];
>> 265 if (aPhoton)
>> 266 {
>> 267 G4double itsCut = cutg;
>> 268 if(aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cute;
>> 269 G4double itsEnergy = aPhoton->GetKineticEnergy();
>> 270
>> 271 if (itsEnergy > itsCut && itsEnergy <= bindingEnergy)
>> 272 {
>> 273 nPhotons++;
>> 274 // Local energy deposit is given as the sum of the
>> 275 // energies of incident photons minus the energies
>> 276 // of the outcoming fluorescence photons
>> 277 bindingEnergy -= itsEnergy;
>> 278
>> 279 }
>> 280 else
>> 281 {
>> 282 delete aPhoton;
>> 283 (*photonVector)[k] = 0;
>> 284 }
>> 285 }
>> 286 }
355 } 287 }
356 else if (gammaEnergy >= (*(fParamLow[Z]))[ << 288
357 G4double x1 = 1.0 / gammaEnergy; << 289 energyDeposit += bindingEnergy;
358 G4double x2 = x1 * x1; << 290
359 G4double x3 = x2 * x1; << 291 // Final state
360 G4double x4 = x3 * x1; << 292
361 G4double x5 = x4 * x1; << 293 for (G4int l = 0; l<nElectrons; l++ )
362 std::size_t idx = nn * 7 - 5; << 294 {
363 // when do sampling common factors are n << 295 aPhoton = electronVector[l];
364 // so cross section is not real << 296 if(aPhoton) {
365 G4double cs0 = G4UniformRand() * << 297 fvect->push_back(aPhoton);
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 } 298 }
383 } 299 }
384 else { << 300 for ( size_t ll = 0; ll < nTotPhotons; ll++)
385 // when do sampling common factors are n << 301 {
386 // so cross section is not real << 302 aPhoton = (*photonVector)[ll];
387 G4double cs = G4UniformRand(); << 303 if(aPhoton) {
388 << 304 fvect->push_back(aPhoton);
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 } <<
403 if (cs <= 0.0 || j + 1 == (G4int)nn) { <<
404 break; <<
405 } <<
406 } 305 }
407 } 306 }
408 } <<
409 // END: SAMPLING OF THE SHELL <<
410 <<
411 G4double bindingEnergy = (*(fParamHigh[Z]))[ <<
412 const G4AtomicShell* shell = nullptr; <<
413 <<
414 // no de-excitation from the last shell <<
415 if (fDeexcitationActive && shellIdx + 1 < nn <<
416 auto as = G4AtomicShellEnumerator(shellIdx <<
417 shell = fAtomDeexcitation->GetAtomicShell( <<
418 } <<
419 307
420 // If binding energy of the selected shell i << 308 delete photonVector;
421 // do not generate secondaries <<
422 if (gammaEnergy < bindingEnergy) { <<
423 fParticleChange->ProposeLocalEnergyDeposit <<
424 return; <<
425 } <<
426 309
427 // Primary outcoming electron << 310 if (energyDeposit < 0)
428 G4double eKineticEnergy = gammaEnergy - bind << 311 {
429 G4double edep = bindingEnergy; << 312 G4cout << "WARNING - "
430 << 313 << "G4LowEnergyPhotoElectric::PostStepDoIt - Negative energy deposit"
431 // Calculate direction of the photoelectron << 314 << G4endl;
432 G4ThreeVector electronDirection = GetAngular << 315 energyDeposit = 0;
433 aDynamicGamma, eKineticEnergy, (G4int)shel <<
434 <<
435 // The electron is created <<
436 auto electron = <<
437 new G4DynamicParticle(theElectron, electro <<
438 fvect->push_back(electron); <<
439 <<
440 // Sample deexcitation <<
441 if (nullptr != shell) { <<
442 G4int index = couple->GetIndex(); <<
443 if (fAtomDeexcitation->CheckDeexcitationAc <<
444 std::size_t nbefore = fvect->size(); <<
445 <<
446 fAtomDeexcitation->GenerateParticles(fve <<
447 std::size_t nafter = fvect->size(); <<
448 if (nafter > nbefore) { <<
449 G4double esec = 0.0; <<
450 for (std::size_t j = nbefore; j < naft <<
451 G4double e = ((*fvect)[j])->GetKinet <<
452 if (esec + e > edep) { <<
453 // correct energy in order to have <<
454 e = edep - esec; <<
455 ((*fvect)[j])->SetKineticEnergy(e) <<
456 esec += e; <<
457 // delete the rest of secondaries <<
458 for (std::size_t jj = nafter - 1; <<
459 delete (*fvect)[jj]; <<
460 fvect->pop_back(); <<
461 } <<
462 break; <<
463 } <<
464 esec += e; <<
465 } <<
466 edep -= esec; <<
467 } <<
468 } 316 }
469 } << 317
470 // energy balance - excitation energy left << 318 // kill incident photon
471 if (edep > 0.0) { << 319 fParticleChange->ProposeMomentumDirection( 0., 0., 0. );
472 fParticleChange->ProposeLocalEnergyDeposit << 320 fParticleChange->SetProposedKineticEnergy(0.);
473 } << 321 fParticleChange->ProposeTrackStatus(fStopAndKill);
>> 322 fParticleChange->ProposeLocalEnergyDeposit(energyDeposit);
474 } 323 }
475 324
476 //....oooOO0OOooo........oooOO0OOooo........oo 325 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
477 326
478 const G4String& G4LivermorePhotoElectricModel: << 327 void G4LivermorePhotoElectricModel::SetCutForLowEnSecPhotons(G4double cut)
479 { 328 {
480 // no check in this method - environment var << 329 cutForLowEnergySecondaryPhotons = cut;
481 if (fDataDirectory.empty()) { << 330 deexcitationManager.SetCutForSecondaryPhotons(cut);
482 auto param = G4EmParameters::Instance(); <<
483 std::ostringstream ost; <<
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 } 331 }
494 332
495 //....oooOO0OOooo........oooOO0OOooo........oo 333 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
496 334
497 void G4LivermorePhotoElectricModel::ReadData(G << 335 void G4LivermorePhotoElectricModel::SetCutForLowEnSecElectrons(G4double cut)
498 { 336 {
499 if (Z <= 0 || Z >= ZMAXPE) { << 337 cutForLowEnergySecondaryElectrons = cut;
500 G4cout << "G4LivermorePhotoElectricModel:: << 338 deexcitationManager.SetCutForAugerElectrons(cut);
501 << Z << " is out of range - request ignor << 339 }
502 return; <<
503 } <<
504 if (verboseLevel > 1) { <<
505 G4cout << "G4LivermorePhotoElectricModel:: <<
506 } <<
507 340
508 if (fCrossSection->GetElementData(Z) != null << 341 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
509 return; <<
510 } <<
511 342
512 // spline for photoeffect total x-section ab << 343 void G4LivermorePhotoElectricModel::ActivateAuger(G4bool val)
513 // but below the parameterized ones << 344 {
>> 345 deexcitationManager.ActivateAugerElectronProduction(val);
>> 346 }
514 347
515 G4bool spline = (G4EmParameters::Instance()- << 348 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
516 G4int number = G4EmParameters::Instance()->N <<
517 auto pv = new G4PhysicsFreeVector(spline); <<
518 <<
519 // fDataDirectory will be defined after thes <<
520 std::ostringstream ost; <<
521 ost << FindDirectoryPath() << "pe-cs-" << Z <<
522 std::ifstream fin(ost.str().c_str()); <<
523 if (!fin.is_open()) { <<
524 G4ExceptionDescription ed; <<
525 ed << "G4LivermorePhotoElectricModel data <<
526 << ost.str().c_str() << "> is not opene <<
527 G4Exception("G4LivermorePhotoElectricModel <<
528 "G4LEDATA version should be G4 <<
529 return; <<
530 } <<
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 <<
542 // read high-energy fit parameters <<
543 fParamHigh[Z] = new std::vector<G4double>; <<
544 G4int n1 = 0; <<
545 G4int n2 = 0; <<
546 G4double x; <<
547 std::ostringstream ost1; <<
548 ost1 << fDataDirectory << "pe-high-" << Z << <<
549 std::ifstream fin1(ost1.str().c_str()); <<
550 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 349
570 fin1 >> n2; << 350 void G4LivermorePhotoElectricModel::SetAngularGenerator(G4VPhotoElectricAngularDistribution* distribution)
571 if (fin1.fail()) { << 351 {
572 return; << 352 ElectronAngularGenerator = distribution;
573 } << 353 ElectronAngularGenerator->PrintGeneratorInformation();
574 if (0 > n2 || n2 >= INT_MAX) { << 354 }
575 n2 = 0; <<
576 } <<
577 355
578 fin1 >> x; << 356 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
579 if (fin1.fail()) { <<
580 return; <<
581 } <<
582 357
583 fNShells[Z] = n1; << 358 void G4LivermorePhotoElectricModel::SetAngularGenerator(const G4String& name)
584 fParamHigh[Z]->reserve(7 * n1 + 1); << 359 {
585 fParamHigh[Z]->push_back(x * MeV); << 360 if (name == "default")
586 for (G4int i = 0; i < n1; ++i) { << 361 {
587 for (G4int j = 0; j < 7; ++j) { << 362 delete ElectronAngularGenerator;
588 fin1 >> x; << 363 ElectronAngularGenerator = new G4PhotoElectricAngularGeneratorSimple("GEANT4LowEnergySimpleGenerator");
589 if (0 == j) { << 364 generatorName = name;
590 x *= MeV; <<
591 } <<
592 else { <<
593 x *= barn; <<
594 } <<
595 fParamHigh[Z]->push_back(x); <<
596 } 365 }
597 } << 366 else if (name == "standard")
598 fin1.close(); << 367 {
599 << 368 delete ElectronAngularGenerator;
600 // read low-energy fit parameters << 369 ElectronAngularGenerator = new G4PhotoElectricAngularGeneratorSauterGavrila("GEANT4SauterGavrilaGenerator");
601 fParamLow[Z] = new std::vector<G4double>; << 370 generatorName = name;
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 } 371 }
655 } << 372 else if (name == "polarized")
656 fin1_low.close(); << 373 {
657 << 374 delete ElectronAngularGenerator;
658 // there is a possibility to use only main s << 375 ElectronAngularGenerator = new G4PhotoElectricAngularGeneratorPolarized("GEANT4LowEnergyPolarizedGenerator");
659 if (nShellLimit < n2) { << 376 generatorName = name;
660 n2 = nShellLimit; << 377 }
661 } << 378 else
662 fCrossSection->InitialiseForComponent(Z, n2) << 379 {
663 fNShellsUsed[Z] = n2; << 380 G4Exception("G4LowEnergyPhotoElectric::SetAngularGenerator - generator does not exist");
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 } <<
681 <<
682 G4int n3, n4; <<
683 G4double y; <<
684 for (G4int i = 0; i < n2; ++i) { <<
685 fin2 >> x >> y >> n3 >> n4; <<
686 auto v = new G4PhysicsFreeVector(n3, x, <<
687 for (G4int j = 0; j < n3; ++j) { <<
688 fin2 >> x >> y; <<
689 v->PutValues(j, x * CLHEP::MeV, y * CL <<
690 } <<
691 v->EnableLogBinSearch(number); <<
692 fCrossSection->AddComponent(Z, n4, v); <<
693 } 381 }
694 fin2.close(); <<
695 } <<
696 382
697 // no spline for photoeffect total x-section << 383 ElectronAngularGenerator->PrintGeneratorInformation();
698 if (1 < fNShells[Z]) { <<
699 auto pv1 = new G4PhysicsFreeVector(false); <<
700 std::ostringstream ost3; <<
701 ost3 << fDataDirectory << "pe-le-cs-" << Z <<
702 std::ifstream fin3(ost3.str().c_str()); <<
703 if (!fin3.is_open()) { <<
704 G4ExceptionDescription ed; <<
705 ed << "G4LivermorePhotoElectricModel dat <<
706 << G4endl; <<
707 G4Exception("G4LivermorePhotoElectricMod <<
708 "G4LEDATA version should be <<
709 return; <<
710 } <<
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 } 384 }
722 385
723 //....oooOO0OOooo........oooOO0OOooo........oo 386 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
724 387
725 G4double G4LivermorePhotoElectricModel::GetBin << 388 G4double G4LivermorePhotoElectricModel::GetMeanFreePath(const G4Track& track,
>> 389 G4double, // previousStepSize
>> 390 G4ForceCondition*)
726 { 391 {
727 if (Z < 1 || Z >= ZMAXPE) { << 392 const G4DynamicParticle* photon = track.GetDynamicParticle();
728 return -1; << 393 G4double energy = photon->GetKineticEnergy();
729 } // If Z is out of the supported return 0 << 394 G4Material* material = track.GetMaterial();
730 << 395 // size_t materialIndex = material->GetIndex();
731 InitialiseOnFly(Z); <<
732 if (fCrossSection->GetElementData(Z) == null <<
733 shell < 0 || shell >= fNShellsUsed[Z]) { <<
734 return -1; <<
735 } <<
736 396
737 if (Z > 2) { << 397 G4double meanFreePath = DBL_MAX;
738 return fCrossSection->GetComponentDataByIn <<
739 } <<
740 else { <<
741 return fCrossSection->GetElementData(Z)->E <<
742 } <<
743 } <<
744 <<
745 //....oooOO0OOooo........oooOO0OOooo........oo <<
746 398
747 void << 399 // if (energy > highEnergyLimit)
748 G4LivermorePhotoElectricModel::InitialiseForEl << 400 // meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
749 { << 401 // else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
750 if (fCrossSection == nullptr) { << 402 // else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
751 fCrossSection = new G4ElementData(ZMAXPE); <<
752 fCrossSection->SetName("PhotoEffXS"); <<
753 fCrossSectionLE = new G4ElementData(ZMAXPE <<
754 fCrossSectionLE->SetName("PhotoEffLowXS"); <<
755 } <<
756 ReadData(Z); <<
757 } <<
758 403
759 //....oooOO0OOooo........oooOO0OOooo........oo << 404 G4double cross = shellCrossSectionHandler->ValueForMaterial(material,energy);
>> 405 if(cross > 0.0) meanFreePath = 1.0/cross;
760 406
761 void G4LivermorePhotoElectricModel::Initialise << 407 return meanFreePath;
762 { <<
763 if (fCrossSection->GetElementData(Z) == null <<
764 G4AutoLock l(&livPhotoeffMutex); <<
765 if (fCrossSection->GetElementData(Z) == nu <<
766 ReadData(Z); <<
767 } <<
768 l.unlock(); <<
769 } <<
770 } 408 }
771 409
772 //....oooOO0OOooo........oooOO0OOooo........oo <<
773 410