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