| 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: G4PenelopePhotoElectricModel.cc,v 1.2.2.1 2009/03/05 08:50:19 gcosmo Exp $
>> 27 // GEANT4 tag $Name: geant4-09-02-patch-01 $
26 // 28 //
27 // Author: Luciano Pandola 29 // Author: Luciano Pandola
28 // 30 //
29 // History: 31 // History:
30 // -------- 32 // --------
31 // 08 Jan 2010 L Pandola First implementati << 33 // 08 Oct 2008 L Pandola Migration from process to model
32 // 01 Feb 2011 L Pandola Suppress fake ener << 34 // 08 Jan 2009 L. Pandola Check shell index to avoid mismatch between
33 // is active. << 35 // the Penelope cross section database and the
34 // Make sure that flu << 36 // G4AtomicTransitionManager database. It suppresses
35 // only if above thre << 37 // a warning from G4AtomicTransitionManager only.
36 // 25 May 2011 L Pandola Renamed (make v200 << 38 // Results are unchanged.
37 // 10 Jun 2011 L Pandola Migrate atomic dee <<
38 // 07 Oct 2011 L Pandola Bug fix (potential <<
39 // 27 Sep 2013 L Pandola Migrate to MT para <<
40 // tables. <<
41 // 02 Oct 2013 L Pandola Rewrite sampling a <<
42 // to improve CPU per <<
43 // 39 //
44 40
45 #include "G4PenelopePhotoElectricModel.hh" 41 #include "G4PenelopePhotoElectricModel.hh"
46 #include "G4PhysicalConstants.hh" <<
47 #include "G4SystemOfUnits.hh" <<
48 #include "G4ParticleDefinition.hh" 42 #include "G4ParticleDefinition.hh"
49 #include "G4MaterialCutsCouple.hh" 43 #include "G4MaterialCutsCouple.hh"
>> 44 #include "G4ProductionCutsTable.hh"
50 #include "G4DynamicParticle.hh" 45 #include "G4DynamicParticle.hh"
51 #include "G4PhysicsTable.hh" 46 #include "G4PhysicsTable.hh"
52 #include "G4PhysicsFreeVector.hh" <<
53 #include "G4ElementTable.hh" 47 #include "G4ElementTable.hh"
54 #include "G4Element.hh" 48 #include "G4Element.hh"
>> 49 #include "G4CrossSectionHandler.hh"
55 #include "G4AtomicTransitionManager.hh" 50 #include "G4AtomicTransitionManager.hh"
56 #include "G4AtomicShell.hh" 51 #include "G4AtomicShell.hh"
57 #include "G4Gamma.hh" 52 #include "G4Gamma.hh"
58 #include "G4Electron.hh" 53 #include "G4Electron.hh"
59 #include "G4AutoLock.hh" <<
60 #include "G4LossTableManager.hh" <<
61 #include "G4Exp.hh" <<
62 54
63 //....oooOO0OOooo........oooOO0OOooo........oo 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64 56
65 const G4int G4PenelopePhotoElectricModel::fMax <<
66 G4PhysicsTable* G4PenelopePhotoElectricModel:: <<
67 57
68 //....oooOO0OOooo........oooOO0OOooo........oo << 58 G4PenelopePhotoElectricModel::G4PenelopePhotoElectricModel(const G4ParticleDefinition*,
69 << 59 const G4String& nam)
70 G4PenelopePhotoElectricModel::G4PenelopePhotoE << 60 :G4VEmModel(nam),isInitialised(false),crossSectionHandler(0),
71 const G4String& nam) << 61 shellCrossSectionHandler(0)
72 :G4VEmModel(nam),fParticleChange(nullptr),fP <<
73 fAtomDeexcitation(nullptr),fIsInitialised(f <<
74 { 62 {
75 fIntrinsicLowEnergyLimit = 100.0*eV; 63 fIntrinsicLowEnergyLimit = 100.0*eV;
76 fIntrinsicHighEnergyLimit = 100.0*GeV; 64 fIntrinsicHighEnergyLimit = 100.0*GeV;
77 // SetLowEnergyLimit(fIntrinsicLowEnergyLim << 65 SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
78 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 66 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
79 // 67 //
80 << 68 fUseAtomicDeexcitation = true;
81 if (part) << 69 verboseLevel= 0;
82 SetParticle(part); <<
83 <<
84 fVerboseLevel= 0; <<
85 // Verbosity scale: 70 // Verbosity scale:
86 // 0 = nothing << 71 // 0 = nothing
87 // 1 = warning for energy non-conservation << 72 // 1 = warning for energy non-conservation
88 // 2 = details of energy budget 73 // 2 = details of energy budget
89 // 3 = calculation of cross sections, file o 74 // 3 = calculation of cross sections, file openings, sampling of atoms
90 // 4 = entering in methods 75 // 4 = entering in methods
91 <<
92 //Mark this model as "applicable" for atomic <<
93 SetDeexcitationFlag(true); <<
94 <<
95 fTransitionManager = G4AtomicTransitionManag <<
96 } 76 }
97 77
98 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
99 79
100 G4PenelopePhotoElectricModel::~G4PenelopePhoto 80 G4PenelopePhotoElectricModel::~G4PenelopePhotoElectricModel()
101 { << 81 {
102 if (IsMaster() || fLocalTable) << 82 if (crossSectionHandler) delete crossSectionHandler;
103 { << 83 if (shellCrossSectionHandler) delete shellCrossSectionHandler;
104 for(G4int i=0; i<=fMaxZ; ++i) <<
105 { <<
106 if(fLogAtomicShellXS[i]) { <<
107 fLogAtomicShellXS[i]->clearAndDestroy(); <<
108 delete fLogAtomicShellXS[i]; <<
109 fLogAtomicShellXS[i] = nullptr; <<
110 } <<
111 } <<
112 } <<
113 } 84 }
114 85
115 //....oooOO0OOooo........oooOO0OOooo........oo 86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
116 87
117 void G4PenelopePhotoElectricModel::Initialise( << 88 void G4PenelopePhotoElectricModel::Initialise(const G4ParticleDefinition*,
118 const G4DataVector& cuts) << 89 const G4DataVector& )
119 { 90 {
120 if (fVerboseLevel > 3) << 91 if (verboseLevel > 3)
121 G4cout << "Calling G4PenelopePhotoElectri 92 G4cout << "Calling G4PenelopePhotoElectricModel::Initialise()" << G4endl;
122 << 93 if (crossSectionHandler)
123 fAtomDeexcitation = G4LossTableManager::Inst <<
124 //Issue warning if the AtomicDeexcitation ha <<
125 if (!fAtomDeexcitation) <<
126 { 94 {
127 G4cout << G4endl; << 95 crossSectionHandler->Clear();
128 G4cout << "WARNING from G4PenelopePhotoE << 96 delete crossSectionHandler;
129 G4cout << "Atomic de-excitation module i <<
130 G4cout << "any fluorescence/Auger emissi <<
131 G4cout << "Please make sure this is inte <<
132 } 97 }
133 << 98 if (shellCrossSectionHandler)
134 SetParticle(particle); <<
135 <<
136 //Only the master model creates/fills/destro <<
137 if (IsMaster() && particle == fParticle) <<
138 { 99 {
139 G4ProductionCutsTable* theCoupleTable = << 100 shellCrossSectionHandler->Clear();
140 G4ProductionCutsTable::GetProductionCutsTabl << 101 delete shellCrossSectionHandler;
141 << 102 }
142 for (G4int i=0;i<(G4int)theCoupleTable-> <<
143 { <<
144 const G4Material* material = <<
145 theCoupleTable->GetMaterialCutsCouple(i) <<
146 const G4ElementVector* theElementVector = <<
147 <<
148 for (std::size_t j=0;j<material->GetNumber <<
149 { <<
150 G4int iZ = theElementVector->at(j)->Ge <<
151 //read data files only in the master <<
152 if (!fLogAtomicShellXS[iZ]) <<
153 ReadDataFile(iZ); <<
154 } <<
155 } <<
156 <<
157 InitialiseElementSelectors(particle,cuts <<
158 103
159 if (fVerboseLevel > 0) { << 104 //Check energy limits
160 G4cout << "Penelope Photo-Electric model v20 << 105 if (LowEnergyLimit() < fIntrinsicLowEnergyLimit)
161 << "Energy range: " << 106 {
162 << LowEnergyLimit() / MeV << " MeV - << 107 G4cout << "G4PenelopePhotoElectricModel: low energy limit increased from " <<
163 << HighEnergyLimit() / GeV << " GeV"; << 108 LowEnergyLimit()/eV << " eV to " << fIntrinsicLowEnergyLimit/eV << " eV" <<
164 } << 109 G4endl;
>> 110 SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
>> 111 }
>> 112 if (HighEnergyLimit() > fIntrinsicHighEnergyLimit)
>> 113 {
>> 114 G4cout << "G4PenelopePhotoElectricModel: high energy limit decreased from " <<
>> 115 HighEnergyLimit()/GeV << " GeV to " << fIntrinsicHighEnergyLimit/GeV << " GeV"
>> 116 << G4endl;
>> 117 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
165 } 118 }
166 119
167 if(fIsInitialised) return; <<
168 fParticleChange = GetParticleChangeForGamma( <<
169 fIsInitialised = true; <<
170 120
171 } << 121 //Re-initialize cross section handlers
>> 122 crossSectionHandler = new G4CrossSectionHandler();
>> 123 crossSectionHandler->Clear();
>> 124 G4String crossSectionFile = "penelope/ph-cs-pen-";
>> 125 crossSectionHandler->LoadData(crossSectionFile);
>> 126 shellCrossSectionHandler = new G4CrossSectionHandler();
>> 127 shellCrossSectionHandler->Clear();
>> 128 crossSectionFile = "penelope/ph-ss-cs-pen-";
>> 129 shellCrossSectionHandler->LoadShellData(crossSectionFile);
>> 130 //This is used to retrieve cross section values later on
>> 131 crossSectionHandler->BuildMeanFreePathForMaterials();
172 132
173 void G4PenelopePhotoElectricModel::InitialiseL << 133 if (verboseLevel > 2)
174 G4VEmModel *masterModel) << 134 G4cout << "Loaded cross section files for PenelopePhotoElectric" << G4endl;
175 { <<
176 if (fVerboseLevel > 3) <<
177 G4cout << "Calling G4PenelopePhotoElectri <<
178 // <<
179 //Check that particle matches: one might hav <<
180 //for e+ and e-). <<
181 // <<
182 if (part == fParticle) <<
183 { <<
184 SetElementSelectors(masterModel->GetElem <<
185 135
186 //Get the const table pointers from the << 136 G4cout << "Penelope Photo-Electric model is initialized " << G4endl
187 const G4PenelopePhotoElectricModel* theM << 137 << "Energy range: "
188 static_cast<G4PenelopePhotoElectricModel*> ( << 138 << LowEnergyLimit() / MeV << " MeV - "
189 for(G4int i=0; i<=fMaxZ; ++i) << 139 << HighEnergyLimit() / GeV << " GeV"
190 fLogAtomicShellXS[i] = theModel->fLogAtomicS << 140 << G4endl;
191 //Same verbosity for all workers, as the << 141
192 fVerboseLevel = theModel->fVerboseLevel; << 142 if(isInitialised) return;
193 } <<
194 143
195 return; << 144 if(pParticleChange)
>> 145 fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
>> 146 else
>> 147 fParticleChange = new G4ParticleChangeForGamma();
>> 148 isInitialised = true;
196 } 149 }
197 150
198 //....oooOO0OOooo........oooOO0OOooo........oo 151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
199 namespace { G4Mutex PenelopePhotoElectricMode << 152
200 G4double G4PenelopePhotoElectricModel::Compute 153 G4double G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom(
201 const G4ParticleDefinition*, << 154 const G4ParticleDefinition*,
202 G4double energy, << 155 G4double energy,
203 G4double Z, G4double, << 156 G4double Z, G4double,
204 G4double, G4double) << 157 G4double, G4double)
205 { 158 {
206 // 159 //
207 // Penelope model v2008 << 160 // Penelope model. Use data-driven approach for cross section estimate (and
>> 161 // also shell sampling from a given atom). Data are from the Livermore database
>> 162 // D.E. Cullen et al., Report UCRL-50400 (1989)
208 // 163 //
209 if (fVerboseLevel > 3) <<
210 G4cout << "Calling ComputeCrossSectionPerA <<
211 164
212 G4int iZ = G4int(Z); << 165 if (verboseLevel > 3)
>> 166 G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopePhotoElectricModel" << G4endl;
213 167
214 if (!fLogAtomicShellXS[iZ]) << 168 G4int iZ = (G4int) Z;
>> 169 if (!crossSectionHandler)
215 { 170 {
216 //If we are here, it means that Initiali << 171 G4cout << "G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom" << G4endl;
217 //not filled up. This can happen in a Un << 172 G4cout << "The cross section handler is not correctly initialized" << G4endl;
218 if (fVerboseLevel > 0) << 173 G4Exception();
219 { << 174 }
220 //Issue a G4Exception (warning) only in ve << 175 G4double cs = crossSectionHandler->FindValue(iZ,energy);
221 G4ExceptionDescription ed; << 176
222 ed << "Unable to retrieve the shell cross << 177 if (verboseLevel > 2)
223 ed << "This can happen only in Unit Tests <<
224 G4Exception("G4PenelopePhotoElectricModel: <<
225 "em2038",JustWarning,ed); <<
226 } <<
227 //protect file reading via autolock <<
228 G4AutoLock lock(&PenelopePhotoElectricMo <<
229 ReadDataFile(iZ); <<
230 lock.unlock(); <<
231 } <<
232 <<
233 G4double cross = 0; <<
234 G4PhysicsTable* theTable = fLogAtomicShellX <<
235 G4PhysicsFreeVector* totalXSLog = (G4Physics <<
236 <<
237 if (!totalXSLog) <<
238 { <<
239 G4Exception("G4PenelopePhotoElectricMod <<
240 "em2039",FatalException, <<
241 "Unable to retrieve the total cross sec <<
242 return 0; <<
243 } <<
244 G4double logene = G4Log(energy); <<
245 G4double logXS = totalXSLog->Value(logene); <<
246 cross = G4Exp(logXS); <<
247 <<
248 if (fVerboseLevel > 2) <<
249 G4cout << "Photoelectric cross section at 178 G4cout << "Photoelectric cross section at " << energy/MeV << " MeV for Z=" << Z <<
250 " = " << cross/barn << " barn" << G4endl << 179 " = " << cs/barn << " barn" << G4endl;
251 return cross; << 180 return cs;
252 } 181 }
253 182
254 //....oooOO0OOooo........oooOO0OOooo........oo 183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
255 184
256 void G4PenelopePhotoElectricModel::SampleSecon 185 void G4PenelopePhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
257 const G4MaterialCutsCouple* c << 186 const G4MaterialCutsCouple* couple,
258 const G4DynamicParticle* aDyn << 187 const G4DynamicParticle* aDynamicGamma,
259 G4double, << 188 G4double,
260 G4double) << 189 G4double)
261 { 190 {
262 // 191 //
263 // Photoelectric effect, Penelope model v200 << 192 // Photoelectric effect, Penelope model
264 // 193 //
265 // The target atom and the target shell are << 194 // The target atom and the target shell are sampled according to the Livermore
266 // database << 195 // database
267 // D.E. Cullen et al., Report UCRL-50400 (1 196 // D.E. Cullen et al., Report UCRL-50400 (1989)
268 // The angular distribution of the electron << 197 // The angular distribution of the electron in the final state is sampled
269 // according to the Sauter distribution from << 198 // according to the Sauter distribution from
270 // F. Sauter, Ann. Phys. 11 (1931) 454 199 // F. Sauter, Ann. Phys. 11 (1931) 454
271 // The energy of the final electron is given << 200 // The energy of the final electron is given by the initial photon energy minus
272 // the binding energy. Fluorescence de-excit << 201 // the binding energy. Fluorescence de-excitation is subsequently produced
273 // (to fill the vacancy) according to the ge 202 // (to fill the vacancy) according to the general Geant4 G4DeexcitationManager:
274 // J. Stepanek, Comp. Phys. Comm. 1206 pp 1 203 // J. Stepanek, Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
275 204
276 if (fVerboseLevel > 3) << 205 if (verboseLevel > 3)
277 G4cout << "Calling SamplingSecondaries() o 206 G4cout << "Calling SamplingSecondaries() of G4PenelopePhotoElectricModel" << G4endl;
278 207
279 G4double photonEnergy = aDynamicGamma->GetKi 208 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
280 209
281 // always kill primary << 210 if (photonEnergy <= LowEnergyLimit())
282 fParticleChange->ProposeTrackStatus(fStopAnd << 211 {
283 fParticleChange->SetProposedKineticEnergy(0. << 212 fParticleChange->ProposeTrackStatus(fStopAndKill);
284 << 213 fParticleChange->SetProposedKineticEnergy(0.);
285 if (photonEnergy <= fIntrinsicLowEnergyLimit <<
286 { <<
287 fParticleChange->ProposeLocalEnergyDepos 214 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
288 return ; 215 return ;
289 } << 216 }
290 217
291 G4ParticleMomentum photonDirection = aDynami 218 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
292 219
293 // Select randomly one element in the curren 220 // Select randomly one element in the current material
294 if (fVerboseLevel > 2) << 221 if (verboseLevel > 2)
295 G4cout << "Going to select element in " << 222 G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
296 << 223 //use crossSectionHandler instead of G4EmElementSelector because in this case
297 // atom can be selected efficiently if eleme << 224 //the dimension of the table is equal to the dimension of the database (less interpolation errors)
298 const G4Element* anElement = << 225 G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
299 SelectRandomAtom(couple,G4Gamma::GammaDefi << 226 if (verboseLevel > 2)
300 G4int Z = anElement->GetZasInt(); << 227 G4cout << "Selected Z = " << Z << G4endl;
301 if (fVerboseLevel > 2) <<
302 G4cout << "Selected " << anElement->GetNam <<
303 228
304 // Select the ionised shell in the current a 229 // Select the ionised shell in the current atom according to shell cross sections
305 //shellIndex = 0 --> K shell << 230 size_t shellIndex = shellCrossSectionHandler->SelectRandomShell(Z,photonEnergy);
306 // 1-3 --> L shells <<
307 // 4-8 --> M shells <<
308 // 9 --> outer shells cumulative <<
309 // <<
310 std::size_t shellIndex = SelectRandomShell(Z <<
311 <<
312 if (fVerboseLevel > 2) <<
313 G4cout << "Selected shell " << shellIndex <<
314 231
315 // Retrieve the corresponding identifier and 232 // Retrieve the corresponding identifier and binding energy of the selected shell
316 const G4AtomicTransitionManager* transitionM 233 const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
317 234
318 //The number of shell cross section possibly << 235 //The number of shell cross section possibly reported in the Penelope database
319 //might be different from the number of shel 236 //might be different from the number of shells in the G4AtomicTransitionManager
320 //(namely, Penelope may contain more shell, 237 //(namely, Penelope may contain more shell, especially for very light elements).
321 //In order to avoid a warning message from t << 238 //In order to avoid a warning message from the G4AtomicTransitionManager, I
322 //add this protection. Results are anyway ch 239 //add this protection. Results are anyway changed, because when G4AtomicTransitionManager
323 //has a shellID>maxID, it sets the shellID t << 240 //has a shellID>maxID, it sets the shellID to the last valid shell.
324 std::size_t numberOfShells = (std::size_t) t << 241 size_t numberOfShells = (size_t) transitionManager->NumberOfShells(Z);
325 if (shellIndex >= numberOfShells) 242 if (shellIndex >= numberOfShells)
326 shellIndex = numberOfShells-1; 243 shellIndex = numberOfShells-1;
327 244
328 const G4AtomicShell* shell = fTransitionMana << 245 const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
329 G4double bindingEnergy = shell->BindingEnerg 246 G4double bindingEnergy = shell->BindingEnergy();
330 << 247 G4int shellId = shell->ShellId();
331 //Penelope considers only K, L and M shells. <<
332 //not included in the Penelope database. If <<
333 //shellIndex = 9, it means that an outer she <<
334 //Penelope recipe is to set bindingEnergy = <<
335 //to the electron) and to disregard fluoresc <<
336 if (shellIndex == 9) <<
337 bindingEnergy = 0.*eV; <<
338 248
339 G4double localEnergyDeposit = 0.0; 249 G4double localEnergyDeposit = 0.0;
340 G4double cosTheta = 1.0; <<
341 250
342 // Primary outcoming electron 251 // Primary outcoming electron
343 G4double eKineticEnergy = photonEnergy - bin 252 G4double eKineticEnergy = photonEnergy - bindingEnergy;
>> 253
>> 254 G4double cutForLowEnergySecondaryParticles = 250.0*eV;
>> 255 const G4ProductionCutsTable* theCoupleTable=
>> 256 G4ProductionCutsTable::GetProductionCutsTable();
>> 257 size_t indx = couple->GetIndex();
>> 258 G4double cutE = (*(theCoupleTable->GetEnergyCutsVector(1)))[indx];
>> 259 cutE = std::max(cutForLowEnergySecondaryParticles,cutE);
344 260
345 // There may be cases where the binding ener 261 // There may be cases where the binding energy of the selected shell is > photon energy
346 // In such cases do not generate secondaries 262 // In such cases do not generate secondaries
347 if (eKineticEnergy > 0.) 263 if (eKineticEnergy > 0.)
348 { 264 {
349 // The electron is created << 265 //Now check if the electron is above cuts: if so, it is created explicitely
350 // Direction sampled from the Sauter dis << 266 if (eKineticEnergy > cutE)
351 cosTheta = SampleElectronDirection(eKine << 267 {
352 G4double sinTheta = std::sqrt(1-cosTheta << 268 // The electron is created
353 G4double phi = twopi * G4UniformRand() ; << 269 // Direction sampled from the Sauter distribution
354 G4double dirx = sinTheta * std::cos(phi) << 270 G4double cosTheta = SampleElectronDirection(eKineticEnergy);
355 G4double diry = sinTheta * std::sin(phi) << 271 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
356 G4double dirz = cosTheta ; << 272 G4double phi = twopi * G4UniformRand() ;
357 G4ThreeVector electronDirection(dirx,dir << 273 G4double dirx = sinTheta * std::cos(phi);
358 electronDirection.rotateUz(photonDirecti << 274 G4double diry = sinTheta * std::sin(phi);
359 G4DynamicParticle* electron = new G4Dyna << 275 G4double dirz = cosTheta ;
360 electronDirection, << 276 G4ThreeVector electronDirection(dirx,diry,dirz); //electron direction
361 eKineticEnergy); << 277 electronDirection.rotateUz(photonDirection);
362 fvect->push_back(electron); << 278 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
>> 279 electronDirection,
>> 280 eKineticEnergy);
>> 281 fvect->push_back(electron);
>> 282 }
>> 283 else
>> 284 localEnergyDeposit += eKineticEnergy;
363 } 285 }
364 else 286 else
365 bindingEnergy = photonEnergy; << 287 bindingEnergy = photonEnergy;
366 << 288
>> 289
367 G4double energyInFluorescence = 0; //testing 290 G4double energyInFluorescence = 0; //testing purposes
368 G4double energyInAuger = 0; //testing purpos <<
369 291
370 //Now, take care of fluorescence, if require << 292 //Now, take care of fluorescence, if required
371 //recipe, I have to skip fluoresence complet << 293 if (fUseAtomicDeexcitation)
372 //(= sampling of a shell outer than K,L,M) <<
373 if (fAtomDeexcitation && shellIndex<9) <<
374 { 294 {
375 G4int index = couple->GetIndex(); << 295 G4double cutG = (*(theCoupleTable->GetEnergyCutsVector(0)))[indx];
376 if (fAtomDeexcitation->CheckDeexcitation << 296 cutG = std::min(cutForLowEnergySecondaryParticles,cutG);
>> 297
>> 298 std::vector<G4DynamicParticle*>* photonVector = 0;
>> 299
>> 300 // Protection to avoid generating photons in the unphysical case of
>> 301 // shell binding energy > photon energy
>> 302 if (Z > 5 && (bindingEnergy > cutG || bindingEnergy > cutE))
377 { 303 {
378 std::size_t nBefore = fvect->size(); << 304 photonVector = deexcitationManager.GenerateParticles(Z,shellId);
379 fAtomDeexcitation->GenerateParticles(fvect << 305 //Check for single photons (if they are above threshold)
380 std::size_t nAfter = fvect->size(); << 306 for (size_t k=0; k< photonVector->size(); k++)
381 <<
382 if (nAfter > nBefore) //actual production <<
383 { 307 {
384 for (std::size_t j=nBefore;j<nAfter;++ << 308 G4DynamicParticle* aPhoton = (*photonVector)[k];
>> 309 if (aPhoton)
385 { 310 {
386 G4double itsEnergy = ((*fvect)[j])->GetK << 311 G4double itsCut = cutG;
387 if (itsEnergy < bindingEnergy) // valid << 312 if(aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cutE;
>> 313 G4double itsEnergy = aPhoton->GetKineticEnergy();
>> 314
>> 315 if (itsEnergy > itsCut && itsEnergy <= bindingEnergy)
388 { 316 {
>> 317 // Local energy deposit is given as the sum of the
>> 318 // energies of incident photons minus the energies
>> 319 // of the outcoming fluorescence photons
389 bindingEnergy -= itsEnergy; 320 bindingEnergy -= itsEnergy;
390 if (((*fvect)[j])->GetParticleDefini << 321
391 energyInFluorescence += itsEnergy; << 322 }
392 else if (((*fvect)[j])->GetParticleD << 323 else
393 energyInAuger += itsEnergy; << 324 {
>> 325 (*photonVector)[k] = 0;
394 } 326 }
395 else //invalid secondary: takes more tha <<
396 { <<
397 delete (*fvect)[j]; <<
398 (*fvect)[j] = nullptr; <<
399 } <<
400 } 327 }
401 } 328 }
402 } 329 }
>> 330 //Register valid secondaries
>> 331 if (photonVector)
>> 332 {
>> 333 for ( size_t ll = 0; ll <photonVector->size(); ll++)
>> 334 {
>> 335 G4DynamicParticle* aPhoton = (*photonVector)[ll];
>> 336 if (aPhoton)
>> 337 {
>> 338 energyInFluorescence += aPhoton->GetKineticEnergy();
>> 339 fvect->push_back(aPhoton);
>> 340 }
>> 341 }
>> 342 delete photonVector;
>> 343 }
403 } 344 }
404 <<
405 //Residual energy is deposited locally 345 //Residual energy is deposited locally
406 localEnergyDeposit += bindingEnergy; 346 localEnergyDeposit += bindingEnergy;
>> 347
407 348
408 if (localEnergyDeposit < 0) //Should not be: << 349 if (localEnergyDeposit < 0)
409 { 350 {
410 G4Exception("G4PenelopePhotoElectricMode << 351 G4cout << "WARNING - "
411 "em2099",JustWarning,"WARNING: Negative << 352 << "G4PenelopePhotoElectric::PostStepDoIt - Negative energy deposit"
>> 353 << G4endl;
412 localEnergyDeposit = 0; 354 localEnergyDeposit = 0;
413 } 355 }
414 356
>> 357 //Update the status of the primary gamma (kill it)
>> 358 fParticleChange->SetProposedKineticEnergy(0.);
415 fParticleChange->ProposeLocalEnergyDeposit(l 359 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
>> 360 fParticleChange->ProposeTrackStatus(fStopAndKill);
416 361
417 if (fVerboseLevel > 1) << 362 if (verboseLevel > 1)
418 { 363 {
419 G4cout << "----------------------------- 364 G4cout << "-----------------------------------------------------------" << G4endl;
420 G4cout << "Energy balance from G4Penelop 365 G4cout << "Energy balance from G4PenelopePhotoElectric" << G4endl;
421 G4cout << "Selected shell: " << WriteTar <<
422 anElement->GetName() << G4endl; <<
423 G4cout << "Incoming photon energy: " << 366 G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
424 G4cout << "----------------------------- 367 G4cout << "-----------------------------------------------------------" << G4endl;
425 if (eKineticEnergy) 368 if (eKineticEnergy)
426 G4cout << "Outgoing electron " << eKineticEn 369 G4cout << "Outgoing electron " << eKineticEnergy/keV << " keV" << G4endl;
427 if (energyInFluorescence) << 370 G4cout << "Fluorescence: " << energyInFluorescence/keV << " keV" << G4endl;
428 G4cout << "Fluorescence x-rays: " << energyI <<
429 if (energyInAuger) <<
430 G4cout << "Auger electrons: " << energyInAug <<
431 G4cout << "Local energy deposit " << loc 371 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
432 G4cout << "Total final state: " << << 372 G4cout << "Total final state: " << (eKineticEnergy+energyInFluorescence+localEnergyDeposit)/keV <<
433 (eKineticEnergy+energyInFluorescence+localEn <<
434 " keV" << G4endl; 373 " keV" << G4endl;
435 G4cout << "----------------------------- 374 G4cout << "-----------------------------------------------------------" << G4endl;
436 } 375 }
437 if (fVerboseLevel > 0) << 376 if (verboseLevel > 0)
438 { 377 {
439 G4double energyDiff = << 378 G4double energyDiff = std::fabs(eKineticEnergy+energyInFluorescence+localEnergyDeposit-photonEnergy);
440 std::fabs(eKineticEnergy+energyInFluorescenc <<
441 if (energyDiff > 0.05*keV) 379 if (energyDiff > 0.05*keV)
442 { << 380 G4cout << "Warning from G4PenelopePhotoElectric: problem with energy conservation: " <<
443 G4cout << "Warning from G4PenelopePhotoEle << 381 (eKineticEnergy+energyInFluorescence+localEnergyDeposit)/keV << " keV (final) vs. " <<
444 (eKineticEnergy+energyInFluorescence+loc << 382 photonEnergy/keV << " keV (initial)" << G4endl;
445 << " keV (final) vs. " << <<
446 photonEnergy/keV << " keV (initial)" << <<
447 G4cout << "------------------------------- <<
448 G4cout << "Energy balance from G4PenelopeP <<
449 G4cout << "Selected shell: " << WriteTarge <<
450 anElement->GetName() << G4endl; <<
451 G4cout << "Incoming photon energy: " << ph <<
452 G4cout << "------------------------------- <<
453 if (eKineticEnergy) <<
454 G4cout << "Outgoing electron " << eKinet <<
455 if (energyInFluorescence) <<
456 G4cout << "Fluorescence x-rays: " << ene <<
457 if (energyInAuger) <<
458 G4cout << "Auger electrons: " << energyI <<
459 G4cout << "Local energy deposit " << local <<
460 G4cout << "Total final state: " << <<
461 (eKineticEnergy+energyInFluorescence+loc <<
462 " keV" << G4endl; <<
463 G4cout << "------------------------------- <<
464 } <<
465 } 383 }
466 } 384 }
467 385
468 //....oooOO0OOooo........oooOO0OOooo........oo 386 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
469 387
>> 388 void G4PenelopePhotoElectricModel::ActivateAuger(G4bool augerbool)
>> 389 {
>> 390 if (!fUseAtomicDeexcitation)
>> 391 {
>> 392 G4cout << "WARNING - G4PenelopePhotoElectricModel" << G4endl;
>> 393 G4cout << "The use of the Atomic Deexcitation Manager is set to false " << G4endl;
>> 394 G4cout << "Therefore, Auger electrons will be not generated anyway" << G4endl;
>> 395 }
>> 396 deexcitationManager.ActivateAugerElectronProduction(augerbool);
>> 397 if (verboseLevel > 1)
>> 398 G4cout << "Auger production set to " << augerbool << G4endl;
>> 399 }
>> 400
>> 401 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 402
470 G4double G4PenelopePhotoElectricModel::SampleE 403 G4double G4PenelopePhotoElectricModel::SampleElectronDirection(G4double energy)
471 { 404 {
472 G4double costheta = 1.0; 405 G4double costheta = 1.0;
473 if (energy>1*GeV) return costheta; 406 if (energy>1*GeV) return costheta;
474 << 407
475 //1) initialize energy-dependent variables 408 //1) initialize energy-dependent variables
476 // Variable naming according to Eq. (2.24) o 409 // Variable naming according to Eq. (2.24) of Penelope Manual
477 // (pag. 44) 410 // (pag. 44)
478 G4double gamma = 1.0 + energy/electron_mass_ 411 G4double gamma = 1.0 + energy/electron_mass_c2;
479 G4double gamma2 = gamma*gamma; 412 G4double gamma2 = gamma*gamma;
480 G4double beta = std::sqrt((gamma2-1.0)/gamma 413 G4double beta = std::sqrt((gamma2-1.0)/gamma2);
481 << 414
482 // ac corresponds to "A" of Eq. (2.31) 415 // ac corresponds to "A" of Eq. (2.31)
483 // 416 //
484 G4double ac = (1.0/beta) - 1.0; 417 G4double ac = (1.0/beta) - 1.0;
485 G4double a1 = 0.5*beta*gamma*(gamma-1.0)*(ga 418 G4double a1 = 0.5*beta*gamma*(gamma-1.0)*(gamma-2.0);
486 G4double a2 = ac + 2.0; 419 G4double a2 = ac + 2.0;
487 G4double gtmax = 2.0*(a1 + 1.0/ac); 420 G4double gtmax = 2.0*(a1 + 1.0/ac);
488 << 421
489 G4double tsam = 0; 422 G4double tsam = 0;
490 G4double gtr = 0; 423 G4double gtr = 0;
491 424
492 //2) sampling. Eq. (2.31) of Penelope Manual 425 //2) sampling. Eq. (2.31) of Penelope Manual
493 // tsam = 1-std::cos(theta) 426 // tsam = 1-std::cos(theta)
494 // gtr = rejection function according to Eq. 427 // gtr = rejection function according to Eq. (2.28)
495 do{ 428 do{
496 G4double rand = G4UniformRand(); 429 G4double rand = G4UniformRand();
497 tsam = 2.0*ac * (2.0*rand + a2*std::sqrt(r 430 tsam = 2.0*ac * (2.0*rand + a2*std::sqrt(rand)) / (a2*a2 - 4.0*rand);
498 gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam)); 431 gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam));
499 }while(G4UniformRand()*gtmax > gtr); 432 }while(G4UniformRand()*gtmax > gtr);
500 costheta = 1.0-tsam; 433 costheta = 1.0-tsam;
501 <<
502 return costheta; 434 return costheta;
503 } 435 }
504 436
505 //....oooOO0OOooo........oooOO0OOooo........oo <<
506 <<
507 void G4PenelopePhotoElectricModel::ReadDataFil <<
508 { <<
509 if (!IsMaster()) <<
510 //Should not be here! <<
511 G4Exception("G4PenelopePhotoElectricModel: <<
512 "em0100",FatalException,"Worker thread in <<
513 <<
514 if (fVerboseLevel > 2) <<
515 { <<
516 G4cout << "G4PenelopePhotoElectricModel: <<
517 G4cout << "Going to read PhotoElectric d <<
518 } <<
519 <<
520 const char* path = G4FindDataDir("G4LEDATA <<
521 if(!path) <<
522 { <<
523 G4String excep = "G4PenelopePhotoElectri <<
524 G4Exception("G4PenelopePhotoElectricMode <<
525 "em0006",FatalException,excep); <<
526 return; <<
527 } <<
528 <<
529 /* <<
530 Read the cross section file <<
531 */ <<
532 std::ostringstream ost; <<
533 if (Z>9) <<
534 ost << path << "/penelope/photoelectric/pd <<
535 else <<
536 ost << path << "/penelope/photoelectric/pd <<
537 std::ifstream file(ost.str().c_str()); <<
538 if (!file.is_open()) <<
539 { <<
540 G4String excep = "G4PenelopePhotoElectri <<
541 G4Exception("G4PenelopePhotoElectricMode <<
542 "em0003",FatalException,excep); <<
543 } <<
544 //I have to know in advance how many points <<
545 //to initialize the G4PhysicsFreeVector() <<
546 std::size_t ndata=0; <<
547 G4String line; <<
548 while( getline(file, line) ) <<
549 ndata++; <<
550 ndata -= 1; <<
551 //G4cout << "Found: " << ndata << " lines" < <<
552 <<
553 file.clear(); <<
554 file.close(); <<
555 file.open(ost.str().c_str()); <<
556 <<
557 G4int readZ =0; <<
558 std::size_t nShells= 0; <<
559 file >> readZ >> nShells; <<
560 <<
561 if (fVerboseLevel > 3) <<
562 G4cout << "Element Z=" << Z << " , nShells <<
563 <<
564 //check the right file is opened. <<
565 if (readZ != Z || nShells <= 0 || nShells > <<
566 { <<
567 G4ExceptionDescription ed; <<
568 ed << "Corrupted data file for Z=" << Z <<
569 G4Exception("G4PenelopePhotoElectricMode <<
570 "em0005",FatalException,ed); <<
571 return; <<
572 } <<
573 G4PhysicsTable* thePhysicsTable = new G4Phys <<
574 <<
575 //the table has to contain nShell+1 G4Physic <<
576 //(theTable)[0] --> total cross section <<
577 //(theTable)[ishell] --> cross section for s <<
578 <<
579 //reserve space for the vectors <<
580 //everything is log-log <<
581 for (std::size_t i=0;i<nShells+1;++i) <<
582 thePhysicsTable->push_back(new G4PhysicsFr <<
583 <<
584 std::size_t k =0; <<
585 for (k=0;k<ndata && !file.eof();++k) <<
586 { <<
587 G4double energy = 0; <<
588 G4double aValue = 0; <<
589 file >> energy ; <<
590 energy *= eV; <<
591 G4double logene = G4Log(energy); <<
592 //loop on the columns <<
593 for (std::size_t i=0;i<nShells+1;++i) <<
594 { <<
595 file >> aValue; <<
596 aValue *= barn; <<
597 G4PhysicsFreeVector* theVec = (G4PhysicsFr <<
598 if (aValue < 1e-40*cm2) //protection again <<
599 aValue = 1e-40*cm2; <<
600 theVec->PutValue(k,logene,G4Log(aValue)); <<
601 } <<
602 } <<
603 <<
604 if (fVerboseLevel > 2) <<
605 { <<
606 G4cout << "G4PenelopePhotoElectricModel: <<
607 << Z << G4endl; <<
608 } <<
609 <<
610 fLogAtomicShellXS[Z] = thePhysicsTable; <<
611 <<
612 file.close(); <<
613 return; <<
614 } <<
615 <<
616 //....oooOO0OOooo........oooOO0OOooo........oo <<
617 <<
618 std::size_t G4PenelopePhotoElectricModel::GetN <<
619 { <<
620 if (!IsMaster()) <<
621 //Should not be here! <<
622 G4Exception("G4PenelopePhotoElectricModel: <<
623 "em0100",FatalException,"Worker thread in <<
624 <<
625 //read data files <<
626 if (!fLogAtomicShellXS[Z]) <<
627 ReadDataFile(Z); <<
628 //now it should be ok <<
629 if (!fLogAtomicShellXS[Z]) <<
630 { <<
631 G4ExceptionDescription ed; <<
632 ed << "Cannot find shell cross section <<
633 G4Exception("G4PenelopePhotoElectricMod <<
634 "em2038",FatalException,ed); <<
635 } <<
636 //one vector is allocated for the _total_ cr <<
637 std::size_t nEntries = fLogAtomicShellXS[Z]- <<
638 return (nEntries-1); <<
639 } <<
640 <<
641 //....oooOO0OOooo........oooOO0OOooo........oo <<
642 <<
643 G4double G4PenelopePhotoElectricModel::GetShel <<
644 { <<
645 //this forces also the loading of the data <<
646 std::size_t entries = GetNumberOfShellXS(Z); <<
647 <<
648 if (shellID >= entries) <<
649 { <<
650 G4cout << "Element Z=" << Z << " has dat <<
651 G4cout << "so shellID should be from 0 t <<
652 return 0; <<
653 } <<
654 <<
655 G4PhysicsTable* theTable = fLogAtomicShellX <<
656 //[0] is the total XS, shellID is in the ele <<
657 G4PhysicsFreeVector* totalXSLog = (G4Physics <<
658 <<
659 if (!totalXSLog) <<
660 { <<
661 G4Exception("G4PenelopePhotoElectricMod <<
662 "em2039",FatalException, <<
663 "Unable to retrieve the total cross sec <<
664 return 0; <<
665 } <<
666 G4double logene = G4Log(energy); <<
667 G4double logXS = totalXSLog->Value(logene); <<
668 G4double cross = G4Exp(logXS); <<
669 if (cross < 2e-40*cm2) cross = 0; <<
670 return cross; <<
671 } <<
672 <<
673 //....oooOO0OOooo........oooOO0OOooo........oo <<
674 <<
675 G4String G4PenelopePhotoElectricModel::WriteTa <<
676 { <<
677 G4String theShell = "outer shell"; <<
678 if (shellID == 0) <<
679 theShell = "K"; <<
680 else if (shellID == 1) <<
681 theShell = "L1"; <<
682 else if (shellID == 2) <<
683 theShell = "L2"; <<
684 else if (shellID == 3) <<
685 theShell = "L3"; <<
686 else if (shellID == 4) <<
687 theShell = "M1"; <<
688 else if (shellID == 5) <<
689 theShell = "M2"; <<
690 else if (shellID == 6) <<
691 theShell = "M3"; <<
692 else if (shellID == 7) <<
693 theShell = "M4"; <<
694 else if (shellID == 8) <<
695 theShell = "M5"; <<
696 <<
697 return theShell; <<
698 } <<
699 <<
700 //....oooOO0OOooo........oooOO0OOooo........oo <<
701 <<
702 void G4PenelopePhotoElectricModel::SetParticle <<
703 { <<
704 if(!fParticle) { <<
705 fParticle = p; <<
706 } <<
707 } <<
708 <<
709 //....oooOO0OOooo........oooOO0OOooo........oo <<
710 <<
711 std::size_t G4PenelopePhotoElectricModel::Sele <<
712 { <<
713 G4double logEnergy = G4Log(energy); <<
714 <<
715 //Check if data have been read (it should be <<
716 if (!fLogAtomicShellXS[Z]) <<
717 { <<
718 G4ExceptionDescription ed; <<
719 ed << "Cannot find shell cross section <<
720 G4Exception("G4PenelopePhotoElectricMod <<
721 "em2038",FatalException,ed); <<
722 } <<
723 <<
724 G4PhysicsTable* theTable = fLogAtomicShellX <<
725 <<
726 G4double sum = 0; <<
727 G4PhysicsFreeVector* totalXSLog = (G4Physics <<
728 G4double logXS = totalXSLog->Value(logEnergy <<
729 G4double totalXS = G4Exp(logXS); <<
730 <<
731 //Notice: totalXS is the total cross section <<
732 //the sum of partialXS's, since these includ <<
733 // <<
734 // Therefore, here one have to consider the <<
735 // an outer shell. Conventionally, it is ind <<
736 // <<
737 G4double random = G4UniformRand()*totalXS; <<
738 <<
739 for (std::size_t k=1;k<theTable->entries();+ <<
740 { <<
741 //Add one shell <<
742 G4PhysicsFreeVector* partialXSLog = (G4P <<
743 G4double logXSLocal = partialXSLog->Valu <<
744 G4double partialXS = G4Exp(logXSLocal); <<
745 sum += partialXS; <<
746 if (random <= sum) <<
747 return k-1; <<
748 } <<
749 //none of the shells K, L, M: return outer s <<
750 return 9; <<
751 } <<
752 437