| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // Author: Luciano Pandola 27 // Author: Luciano Pandola
28 // 28 //
29 // History: 29 // History:
30 // -------- 30 // --------
31 // 27 Jul 2010 L Pandola First complete i 31 // 27 Jul 2010 L Pandola First complete implementation
32 // 18 Jan 2011 L.Pandola Stricter check o 32 // 18 Jan 2011 L.Pandola Stricter check on production of sub-treshold delta-rays.
33 // Should never hap 33 // Should never happen now
34 // 01 Feb 2011 L Pandola Suppress fake ener 34 // 01 Feb 2011 L Pandola Suppress fake energy-violation warning when Auger is active.
35 // Make sure that flu 35 // Make sure that fluorescence/Auger is generated only if
36 // above threshold 36 // above threshold
37 // 25 May 2011 L Pandola Renamed (make v200 37 // 25 May 2011 L Pandola Renamed (make v2008 as default Penelope)
38 // 26 Jan 2012 L Pandola Migration of Atomi 38 // 26 Jan 2012 L Pandola Migration of AtomicDeexcitation to the new interface
39 // 09 Mar 2012 L Pandola Moved the manageme 39 // 09 Mar 2012 L Pandola Moved the management and calculation of
40 // cross sections to 40 // cross sections to a separate class. Use a different method to
41 // get normalized she 41 // get normalized shell cross sections
42 // 07 Oct 2013 L. Pandola Migration to MT 42 // 07 Oct 2013 L. Pandola Migration to MT
43 // 23 Jun 2015 L. Pandola Keep track of the 43 // 23 Jun 2015 L. Pandola Keep track of the PIXE flag, to avoid double-production of
44 // atomic de-excitat 44 // atomic de-excitation (bug #1761)
45 // 29 Aug 2018 L. Pandola Fix bug causing en 45 // 29 Aug 2018 L. Pandola Fix bug causing energy non-conservation
46 // 46 //
47 47
48 #include "G4PenelopeIonisationModel.hh" 48 #include "G4PenelopeIonisationModel.hh"
49 #include "G4PhysicalConstants.hh" 49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh" 50 #include "G4SystemOfUnits.hh"
51 #include "G4ParticleDefinition.hh" 51 #include "G4ParticleDefinition.hh"
52 #include "G4MaterialCutsCouple.hh" 52 #include "G4MaterialCutsCouple.hh"
53 #include "G4ProductionCutsTable.hh" 53 #include "G4ProductionCutsTable.hh"
54 #include "G4DynamicParticle.hh" 54 #include "G4DynamicParticle.hh"
55 #include "G4AtomicTransitionManager.hh" 55 #include "G4AtomicTransitionManager.hh"
56 #include "G4AtomicShell.hh" 56 #include "G4AtomicShell.hh"
57 #include "G4Gamma.hh" 57 #include "G4Gamma.hh"
58 #include "G4Electron.hh" 58 #include "G4Electron.hh"
59 #include "G4Positron.hh" 59 #include "G4Positron.hh"
60 #include "G4PenelopeOscillatorManager.hh" 60 #include "G4PenelopeOscillatorManager.hh"
61 #include "G4PenelopeOscillator.hh" 61 #include "G4PenelopeOscillator.hh"
62 #include "G4PenelopeCrossSection.hh" 62 #include "G4PenelopeCrossSection.hh"
63 #include "G4PhysicsFreeVector.hh" 63 #include "G4PhysicsFreeVector.hh"
64 #include "G4PhysicsLogVector.hh" 64 #include "G4PhysicsLogVector.hh"
65 #include "G4LossTableManager.hh" 65 #include "G4LossTableManager.hh"
66 #include "G4PenelopeIonisationXSHandler.hh" 66 #include "G4PenelopeIonisationXSHandler.hh"
67 #include "G4EmParameters.hh" 67 #include "G4EmParameters.hh"
68 #include "G4AutoLock.hh" 68 #include "G4AutoLock.hh"
69 69
70 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 namespace { G4Mutex PenelopeIonisationModelMu 71 namespace { G4Mutex PenelopeIonisationModelMutex = G4MUTEX_INITIALIZER; }
72 72
73 G4PenelopeIonisationModel::G4PenelopeIonisatio 73 G4PenelopeIonisationModel::G4PenelopeIonisationModel(const G4ParticleDefinition* part,
74 const G4String& nam) 74 const G4String& nam)
75 :G4VEmModel(nam),fParticleChange(nullptr),fP 75 :G4VEmModel(nam),fParticleChange(nullptr),fParticle(nullptr),
76 fCrossSectionHandler(nullptr), 76 fCrossSectionHandler(nullptr),
77 fAtomDeexcitation(nullptr), fKineticEnergy1 77 fAtomDeexcitation(nullptr), fKineticEnergy1(0.*eV),
78 fCosThetaPrimary(1.0),fEnergySecondary(0.*e 78 fCosThetaPrimary(1.0),fEnergySecondary(0.*eV),
79 fCosThetaSecondary(0.0),fTargetOscillator(- 79 fCosThetaSecondary(0.0),fTargetOscillator(-1),
80 fIsInitialised(false),fPIXEflag(false),fLoc 80 fIsInitialised(false),fPIXEflag(false),fLocalTable(false)
81 { 81 {
82 fIntrinsicLowEnergyLimit = 100.0*eV; 82 fIntrinsicLowEnergyLimit = 100.0*eV;
83 fIntrinsicHighEnergyLimit = 100.0*GeV; 83 fIntrinsicHighEnergyLimit = 100.0*GeV;
84 // SetLowEnergyLimit(fIntrinsicLowEnergyLim 84 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
85 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 85 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
86 fNBins = 200; 86 fNBins = 200;
87 87
88 if (part) 88 if (part)
89 SetParticle(part); 89 SetParticle(part);
90 90
91 // 91 //
92 fOscManager = G4PenelopeOscillatorManager::G 92 fOscManager = G4PenelopeOscillatorManager::GetOscillatorManager();
93 // 93 //
94 fVerboseLevel= 0; 94 fVerboseLevel= 0;
95 // Verbosity scale: 95 // Verbosity scale:
96 // 0 = nothing 96 // 0 = nothing
97 // 1 = warning for energy non-conservation 97 // 1 = warning for energy non-conservation
98 // 2 = details of energy budget 98 // 2 = details of energy budget
99 // 3 = calculation of cross sections, file o 99 // 3 = calculation of cross sections, file openings, sampling of atoms
100 // 4 = entering in methods 100 // 4 = entering in methods
101 101
102 // Atomic deexcitation model activated by de 102 // Atomic deexcitation model activated by default
103 SetDeexcitationFlag(true); 103 SetDeexcitationFlag(true);
104 } 104 }
105 105
106 //....oooOO0OOooo........oooOO0OOooo........oo 106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107 107
108 G4PenelopeIonisationModel::~G4PenelopeIonisati 108 G4PenelopeIonisationModel::~G4PenelopeIonisationModel()
109 { 109 {
110 if (IsMaster() || fLocalTable) 110 if (IsMaster() || fLocalTable)
111 { 111 {
112 if (fCrossSectionHandler) 112 if (fCrossSectionHandler)
113 delete fCrossSectionHandler; 113 delete fCrossSectionHandler;
114 } 114 }
115 } 115 }
116 116
117 //....oooOO0OOooo........oooOO0OOooo........oo 117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
118 118
119 void G4PenelopeIonisationModel::Initialise(con 119 void G4PenelopeIonisationModel::Initialise(const G4ParticleDefinition* particle,
120 const G4DataVector& theCuts) 120 const G4DataVector& theCuts)
121 { 121 {
122 if (fVerboseLevel > 3) 122 if (fVerboseLevel > 3)
123 G4cout << "Calling G4PenelopeIonisationMod 123 G4cout << "Calling G4PenelopeIonisationModel::Initialise()" << G4endl;
124 124
125 fAtomDeexcitation = G4LossTableManager::Inst 125 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
126 //Issue warning if the AtomicDeexcitation ha 126 //Issue warning if the AtomicDeexcitation has not been declared
127 if (!fAtomDeexcitation) 127 if (!fAtomDeexcitation)
128 { 128 {
129 G4cout << G4endl; 129 G4cout << G4endl;
130 G4cout << "WARNING from G4PenelopeIonisa 130 G4cout << "WARNING from G4PenelopeIonisationModel " << G4endl;
131 G4cout << "Atomic de-excitation module i 131 G4cout << "Atomic de-excitation module is not instantiated, so there will not be ";
132 G4cout << "any fluorescence/Auger emissi 132 G4cout << "any fluorescence/Auger emission." << G4endl;
133 G4cout << "Please make sure this is inte 133 G4cout << "Please make sure this is intended" << G4endl;
134 } 134 }
135 135
136 if (fAtomDeexcitation) 136 if (fAtomDeexcitation)
137 fPIXEflag = fAtomDeexcitation->IsPIXEActiv 137 fPIXEflag = fAtomDeexcitation->IsPIXEActive();
138 138
139 //If the PIXE flag is active, the PIXE inter 139 //If the PIXE flag is active, the PIXE interface will take care of the
140 //atomic de-excitation. The model does not n 140 //atomic de-excitation. The model does not need to do that.
141 //Issue warnings here 141 //Issue warnings here
142 if (fPIXEflag && IsMaster() && particle==G4E 142 if (fPIXEflag && IsMaster() && particle==G4Electron::Electron())
143 { 143 {
144 G4String theModel = G4EmParameters::Inst 144 G4String theModel = G4EmParameters::Instance()->PIXEElectronCrossSectionModel();
145 G4cout << "============================= 145 G4cout << "======================================================================" << G4endl;
146 G4cout << "The G4PenelopeIonisationModel 146 G4cout << "The G4PenelopeIonisationModel is being used with the PIXE flag ON." << G4endl;
147 G4cout << "Atomic de-excitation will be 147 G4cout << "Atomic de-excitation will be produced statistically by the PIXE " << G4endl;
148 G4cout << "interface by using the shell 148 G4cout << "interface by using the shell cross section --> " << theModel << G4endl;
149 G4cout << "The built-in model procedure 149 G4cout << "The built-in model procedure for atomic de-excitation is disabled. " << G4endl;
150 G4cout << "*Please be sure this is inten 150 G4cout << "*Please be sure this is intended*, or disable PIXE by" << G4endl;
151 G4cout << "/process/em/pixe false" << G4 151 G4cout << "/process/em/pixe false" << G4endl;
152 G4cout << "============================= 152 G4cout << "======================================================================" << G4endl;
153 } 153 }
154 154
155 SetParticle(particle); 155 SetParticle(particle);
156 156
157 //Only the master model creates/manages the 157 //Only the master model creates/manages the tables. All workers get the
158 //pointer to the table, and use it as readon 158 //pointer to the table, and use it as readonly
159 if (IsMaster() && particle == fParticle) 159 if (IsMaster() && particle == fParticle)
160 { 160 {
161 //Set the number of bins for the tables. 161 //Set the number of bins for the tables. 20 points per decade
162 fNBins = (std::size_t) (20*std::log10(Hi 162 fNBins = (std::size_t) (20*std::log10(HighEnergyLimit()/LowEnergyLimit()));
163 fNBins = std::max(fNBins,(std::size_t)10 163 fNBins = std::max(fNBins,(std::size_t)100);
164 164
165 //Clear and re-build the tables 165 //Clear and re-build the tables
166 if (fCrossSectionHandler) 166 if (fCrossSectionHandler)
167 { 167 {
168 delete fCrossSectionHandler; 168 delete fCrossSectionHandler;
169 fCrossSectionHandler = 0; 169 fCrossSectionHandler = 0;
170 } 170 }
171 fCrossSectionHandler = new G4PenelopeIon 171 fCrossSectionHandler = new G4PenelopeIonisationXSHandler(fNBins);
172 fCrossSectionHandler->SetVerboseLevel(fV 172 fCrossSectionHandler->SetVerboseLevel(fVerboseLevel);
173 173
174 //Build tables for all materials 174 //Build tables for all materials
175 G4ProductionCutsTable* theCoupleTable = 175 G4ProductionCutsTable* theCoupleTable =
176 G4ProductionCutsTable::GetProductionCutsTabl 176 G4ProductionCutsTable::GetProductionCutsTable();
177 for (G4int i=0;i<(G4int)theCoupleTable-> 177 for (G4int i=0;i<(G4int)theCoupleTable->GetTableSize();++i)
178 { 178 {
179 const G4Material* theMat = 179 const G4Material* theMat =
180 theCoupleTable->GetMaterialCutsCouple(i) 180 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
181 //Forces the building of the cross section 181 //Forces the building of the cross section tables
182 fCrossSectionHandler->BuildXSTable(theMat, 182 fCrossSectionHandler->BuildXSTable(theMat,theCuts.at(i),particle,
183 IsMaster()); 183 IsMaster());
184 } 184 }
185 185
186 if (fVerboseLevel > 2) { 186 if (fVerboseLevel > 2) {
187 G4cout << "Penelope Ionisation model v2008 i 187 G4cout << "Penelope Ionisation model v2008 is initialized " << G4endl
188 << "Energy range: " 188 << "Energy range: "
189 << LowEnergyLimit() / keV << " keV - 189 << LowEnergyLimit() / keV << " keV - "
190 << HighEnergyLimit() / GeV << " GeV. 190 << HighEnergyLimit() / GeV << " GeV. Using "
191 << fNBins << " bins." 191 << fNBins << " bins."
192 << G4endl; 192 << G4endl;
193 } 193 }
194 } 194 }
195 195
196 if(fIsInitialised) 196 if(fIsInitialised)
197 return; 197 return;
198 fParticleChange = GetParticleChangeForLoss() 198 fParticleChange = GetParticleChangeForLoss();
199 fIsInitialised = true; 199 fIsInitialised = true;
200 200
201 return; 201 return;
202 } 202 }
203 203
204 //....oooOO0OOooo........oooOO0OOooo........oo 204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
205 205
206 void G4PenelopeIonisationModel::InitialiseLoca 206 void G4PenelopeIonisationModel::InitialiseLocal(const G4ParticleDefinition* part,
207 G4VEmModel *masterModel) 207 G4VEmModel *masterModel)
208 { 208 {
209 if (fVerboseLevel > 3) 209 if (fVerboseLevel > 3)
210 G4cout << "Calling G4PenelopeIonisationMo 210 G4cout << "Calling G4PenelopeIonisationModel::InitialiseLocal()" << G4endl;
211 // 211 //
212 //Check that particle matches: one might hav 212 //Check that particle matches: one might have multiple master models (e.g.
213 //for e+ and e-). 213 //for e+ and e-).
214 // 214 //
215 if (part == fParticle) 215 if (part == fParticle)
216 { 216 {
217 //Get the const table pointers from the 217 //Get the const table pointers from the master to the workers
218 const G4PenelopeIonisationModel* theMode 218 const G4PenelopeIonisationModel* theModel =
219 static_cast<G4PenelopeIonisationModel*> (mas 219 static_cast<G4PenelopeIonisationModel*> (masterModel);
220 220
221 //Copy pointers to the data tables 221 //Copy pointers to the data tables
222 fCrossSectionHandler = theModel->fCrossS 222 fCrossSectionHandler = theModel->fCrossSectionHandler;
223 223
224 //copy data 224 //copy data
225 fNBins = theModel->fNBins; 225 fNBins = theModel->fNBins;
226 226
227 //Same verbosity for all workers, as the 227 //Same verbosity for all workers, as the master
228 fVerboseLevel = theModel->fVerboseLevel; 228 fVerboseLevel = theModel->fVerboseLevel;
229 } 229 }
230 return; 230 return;
231 } 231 }
232 232
233 233
234 //....oooOO0OOooo........oooOO0OOooo........oo 234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
235 235
236 G4double G4PenelopeIonisationModel::CrossSecti 236 G4double G4PenelopeIonisationModel::CrossSectionPerVolume(const G4Material* material,
237 const G4ParticleDefinition* 237 const G4ParticleDefinition*
238 theParticle, 238 theParticle,
239 G4double energy, 239 G4double energy,
240 G4double cutEnergy, 240 G4double cutEnergy,
241 G4double) 241 G4double)
242 { 242 {
243 // Penelope model v2008 to calculate the cro 243 // Penelope model v2008 to calculate the cross section for inelastic collisions above the
244 // threshold. It makes use of the Generalise 244 // threshold. It makes use of the Generalised Oscillator Strength (GOS) model from
245 // D. Liljequist, J. Phys. D: Appl. Phys. 1 245 // D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567
246 // 246 //
247 // The total cross section is calculated ana 247 // The total cross section is calculated analytically by taking
248 // into account the atomic oscillators comin 248 // into account the atomic oscillators coming into the play for a given threshold.
249 // 249 //
250 // For incident e- the maximum energy allowe 250 // For incident e- the maximum energy allowed for the delta-rays is energy/2.
251 // because particles are undistinghishable. 251 // because particles are undistinghishable.
252 // 252 //
253 // The contribution is splitted in three par 253 // The contribution is splitted in three parts: distant longitudinal collisions,
254 // distant transverse collisions and close c 254 // distant transverse collisions and close collisions. Each term is described by
255 // its own analytical function. 255 // its own analytical function.
256 // Fermi density correction is calculated an 256 // Fermi density correction is calculated analytically according to
257 // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963), 257 // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1
258 // 258 //
259 if (fVerboseLevel > 3) 259 if (fVerboseLevel > 3)
260 G4cout << "Calling CrossSectionPerVolume() 260 G4cout << "Calling CrossSectionPerVolume() of G4PenelopeIonisationModel" << G4endl;
261 261
262 SetupForMaterial(theParticle, material, ener 262 SetupForMaterial(theParticle, material, energy);
263 263
264 G4double totalCross = 0.0; 264 G4double totalCross = 0.0;
265 G4double crossPerMolecule = 0.; 265 G4double crossPerMolecule = 0.;
266 266
267 //Either Initialize() was not called, or we 267 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
268 //not invoked 268 //not invoked
269 if (!fCrossSectionHandler) 269 if (!fCrossSectionHandler)
270 { 270 {
271 //create a **thread-local** version of t 271 //create a **thread-local** version of the table. Used only for G4EmCalculator and
272 //Unit Tests 272 //Unit Tests
273 fLocalTable = true; 273 fLocalTable = true;
274 fCrossSectionHandler = new G4PenelopeIon 274 fCrossSectionHandler = new G4PenelopeIonisationXSHandler(fNBins);
275 } 275 }
276 276
277 const G4PenelopeCrossSection* theXS = 277 const G4PenelopeCrossSection* theXS =
278 fCrossSectionHandler->GetCrossSectionTable 278 fCrossSectionHandler->GetCrossSectionTableForCouple(theParticle,
279 material, 279 material,
280 cutEnergy); 280 cutEnergy);
281 if (!theXS) 281 if (!theXS)
282 { 282 {
283 //If we are here, it means that Initiali 283 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
284 //not filled up. This can happen in a Un 284 //not filled up. This can happen in a UnitTest or via G4EmCalculator
285 if (fVerboseLevel > 0) 285 if (fVerboseLevel > 0)
286 { 286 {
287 //Issue a G4Exception (warning) only in ve 287 //Issue a G4Exception (warning) only in verbose mode
288 G4ExceptionDescription ed; 288 G4ExceptionDescription ed;
289 ed << "Unable to retrieve the cross sectio 289 ed << "Unable to retrieve the cross section table for " <<
290 theParticle->GetParticleName() << 290 theParticle->GetParticleName() <<
291 " in " << material->GetName() << ", cut 291 " in " << material->GetName() << ", cut = " << cutEnergy/keV << " keV " << G4endl;
292 ed << "This can happen only in Unit Tests 292 ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
293 G4Exception("G4PenelopeIonisationModel::Cr 293 G4Exception("G4PenelopeIonisationModel::CrossSectionPerVolume()",
294 "em2038",JustWarning,ed); 294 "em2038",JustWarning,ed);
295 } 295 }
296 //protect file reading via autolock 296 //protect file reading via autolock
297 G4AutoLock lock(&PenelopeIonisationModel 297 G4AutoLock lock(&PenelopeIonisationModelMutex);
298 fCrossSectionHandler->BuildXSTable(mater 298 fCrossSectionHandler->BuildXSTable(material,cutEnergy,theParticle);
299 lock.unlock(); 299 lock.unlock();
300 //now it should be ok 300 //now it should be ok
301 theXS = 301 theXS =
302 fCrossSectionHandler->GetCrossSectionTableFo 302 fCrossSectionHandler->GetCrossSectionTableForCouple(theParticle,
303 material, 303 material,
304 cutEnergy); 304 cutEnergy);
305 } 305 }
306 306
307 if (theXS) 307 if (theXS)
308 crossPerMolecule = theXS->GetHardCrossSect 308 crossPerMolecule = theXS->GetHardCrossSection(energy);
309 309
310 G4double atomDensity = material->GetTotNbOfA 310 G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
311 G4double atPerMol = fOscManager->GetAtomsPe 311 G4double atPerMol = fOscManager->GetAtomsPerMolecule(material);
312 312
313 if (fVerboseLevel > 3) 313 if (fVerboseLevel > 3)
314 G4cout << "Material " << material->GetName 314 G4cout << "Material " << material->GetName() << " has " << atPerMol <<
315 "atoms per molecule" << G4endl; 315 "atoms per molecule" << G4endl;
316 316
317 G4double moleculeDensity = 0.; 317 G4double moleculeDensity = 0.;
318 if (atPerMol) 318 if (atPerMol)
319 moleculeDensity = atomDensity/atPerMol; 319 moleculeDensity = atomDensity/atPerMol;
320 G4double crossPerVolume = crossPerMolecule*m 320 G4double crossPerVolume = crossPerMolecule*moleculeDensity;
321 321
322 if (fVerboseLevel > 2) 322 if (fVerboseLevel > 2)
323 { 323 {
324 G4cout << "G4PenelopeIonisationModel " << 324 G4cout << "G4PenelopeIonisationModel " << G4endl;
325 G4cout << "Mean free path for delta emissi 325 G4cout << "Mean free path for delta emission > " << cutEnergy/keV << " keV at " <<
326 energy/keV << " keV = " << (1./crossPerV 326 energy/keV << " keV = " << (1./crossPerVolume)/mm << " mm" << G4endl;
327 if (theXS) 327 if (theXS)
328 totalCross = (theXS->GetTotalCrossSectio 328 totalCross = (theXS->GetTotalCrossSection(energy))*moleculeDensity;
329 G4cout << "Total free path for ionisation 329 G4cout << "Total free path for ionisation (no threshold) at " <<
330 energy/keV << " keV = " << (1./totalCros 330 energy/keV << " keV = " << (1./totalCross)/mm << " mm" << G4endl;
331 } 331 }
332 return crossPerVolume; 332 return crossPerVolume;
333 } 333 }
334 334
335 //....oooOO0OOooo........oooOO0OOooo........oo 335 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
336 336
337 //This is a dummy method. Never inkoved by the 337 //This is a dummy method. Never inkoved by the tracking, it just issues
338 //a warning if one tries to get Cross Sections 338 //a warning if one tries to get Cross Sections per Atom via the
339 //G4EmCalculator. 339 //G4EmCalculator.
340 G4double G4PenelopeIonisationModel::ComputeCro 340 G4double G4PenelopeIonisationModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
341 G4double, 341 G4double,
342 G4double, 342 G4double,
343 G4double, 343 G4double,
344 G4double, 344 G4double,
345 G4double) 345 G4double)
346 { 346 {
347 G4cout << "*** G4PenelopeIonisationModel -- 347 G4cout << "*** G4PenelopeIonisationModel -- WARNING ***" << G4endl;
348 G4cout << "Penelope Ionisation model v2008 d 348 G4cout << "Penelope Ionisation model v2008 does not calculate cross section _per atom_ " << G4endl;
349 G4cout << "so the result is always zero. For 349 G4cout << "so the result is always zero. For physics values, please invoke " << G4endl;
350 G4cout << "GetCrossSectionPerVolume() or Get 350 G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl;
351 return 0; 351 return 0;
352 } 352 }
353 353
354 //....oooOO0OOooo........oooOO0OOooo........oo 354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
355 355
356 G4double G4PenelopeIonisationModel::ComputeDED 356 G4double G4PenelopeIonisationModel::ComputeDEDXPerVolume(const G4Material* material,
357 const G4ParticleDefinition* the 357 const G4ParticleDefinition* theParticle,
358 G4double kineticEnergy, 358 G4double kineticEnergy,
359 G4double cutEnergy) 359 G4double cutEnergy)
360 { 360 {
361 // Penelope model v2008 to calculate the sto 361 // Penelope model v2008 to calculate the stopping power for soft inelastic collisions
362 // below the threshold. It makes use of the 362 // below the threshold. It makes use of the Generalised Oscillator Strength (GOS)
363 // model from 363 // model from
364 // D. Liljequist, J. Phys. D: Appl. Phys. 1 364 // D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567
365 // 365 //
366 // The stopping power is calculated analytic 366 // The stopping power is calculated analytically using the dsigma/dW cross
367 // section from the GOS models, which includ 367 // section from the GOS models, which includes separate contributions from
368 // distant longitudinal collisions, distant 368 // distant longitudinal collisions, distant transverse collisions and
369 // close collisions. Only the atomic oscilla 369 // close collisions. Only the atomic oscillators that come in the play
370 // (according to the threshold) are consider 370 // (according to the threshold) are considered for the calculation.
371 // Differential cross sections have a differ 371 // Differential cross sections have a different form for e+ and e-.
372 // 372 //
373 // Fermi density correction is calculated an 373 // Fermi density correction is calculated analytically according to
374 // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963), 374 // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1
375 375
376 if (fVerboseLevel > 3) 376 if (fVerboseLevel > 3)
377 G4cout << "Calling ComputeDEDX() of G4Pene 377 G4cout << "Calling ComputeDEDX() of G4PenelopeIonisationModel" << G4endl;
378 378
379 //Either Initialize() was not called, or we 379 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
380 //not invoked 380 //not invoked
381 if (!fCrossSectionHandler) 381 if (!fCrossSectionHandler)
382 { 382 {
383 //create a **thread-local** version of t 383 //create a **thread-local** version of the table. Used only for G4EmCalculator and
384 //Unit Tests 384 //Unit Tests
385 fLocalTable = true; 385 fLocalTable = true;
386 fCrossSectionHandler = new G4PenelopeIon 386 fCrossSectionHandler = new G4PenelopeIonisationXSHandler(fNBins);
387 } 387 }
388 388
389 const G4PenelopeCrossSection* theXS = 389 const G4PenelopeCrossSection* theXS =
390 fCrossSectionHandler->GetCrossSectionTable 390 fCrossSectionHandler->GetCrossSectionTableForCouple(theParticle,material,
391 cutEnergy); 391 cutEnergy);
392 if (!theXS) 392 if (!theXS)
393 { 393 {
394 //If we are here, it means that Initiali 394 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
395 //not filled up. This can happen in a Un 395 //not filled up. This can happen in a UnitTest or via G4EmCalculator
396 if (fVerboseLevel > 0) 396 if (fVerboseLevel > 0)
397 { 397 {
398 //Issue a G4Exception (warning) only in ve 398 //Issue a G4Exception (warning) only in verbose mode
399 G4ExceptionDescription ed; 399 G4ExceptionDescription ed;
400 ed << "Unable to retrieve the cross sectio 400 ed << "Unable to retrieve the cross section table for " <<
401 theParticle->GetParticleName() << 401 theParticle->GetParticleName() <<
402 " in " << material->GetName() << ", cut 402 " in " << material->GetName() << ", cut = " << cutEnergy/keV << " keV " << G4endl;
403 ed << "This can happen only in Unit Tests 403 ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
404 G4Exception("G4PenelopeIonisationModel::Co 404 G4Exception("G4PenelopeIonisationModel::ComputeDEDXPerVolume()",
405 "em2038",JustWarning,ed); 405 "em2038",JustWarning,ed);
406 } 406 }
407 //protect file reading via autolock 407 //protect file reading via autolock
408 G4AutoLock lock(&PenelopeIonisationModel 408 G4AutoLock lock(&PenelopeIonisationModelMutex);
409 fCrossSectionHandler->BuildXSTable(mater 409 fCrossSectionHandler->BuildXSTable(material,cutEnergy,theParticle);
410 lock.unlock(); 410 lock.unlock();
411 //now it should be ok 411 //now it should be ok
412 theXS = 412 theXS =
413 fCrossSectionHandler->GetCrossSectionTableFo 413 fCrossSectionHandler->GetCrossSectionTableForCouple(theParticle,
414 material, 414 material,
415 cutEnergy); 415 cutEnergy);
416 } 416 }
417 417
418 G4double sPowerPerMolecule = 0.0; 418 G4double sPowerPerMolecule = 0.0;
419 if (theXS) 419 if (theXS)
420 sPowerPerMolecule = theXS->GetSoftStopping 420 sPowerPerMolecule = theXS->GetSoftStoppingPower(kineticEnergy);
421 421
422 G4double atomDensity = material->GetTotNbOfA 422 G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
423 G4double atPerMol = fOscManager->GetAtomsPe 423 G4double atPerMol = fOscManager->GetAtomsPerMolecule(material);
424 424
425 G4double moleculeDensity = 0.; 425 G4double moleculeDensity = 0.;
426 if (atPerMol) 426 if (atPerMol)
427 moleculeDensity = atomDensity/atPerMol; 427 moleculeDensity = atomDensity/atPerMol;
428 G4double sPowerPerVolume = sPowerPerMolecule 428 G4double sPowerPerVolume = sPowerPerMolecule*moleculeDensity;
429 429
430 if (fVerboseLevel > 2) 430 if (fVerboseLevel > 2)
431 { 431 {
432 G4cout << "G4PenelopeIonisationModel " < 432 G4cout << "G4PenelopeIonisationModel " << G4endl;
433 G4cout << "Stopping power < " << cutEner 433 G4cout << "Stopping power < " << cutEnergy/keV << " keV at " <<
434 kineticEnergy/keV << " keV = " << 434 kineticEnergy/keV << " keV = " <<
435 sPowerPerVolume/(keV/mm) << " keV/mm" << G4e 435 sPowerPerVolume/(keV/mm) << " keV/mm" << G4endl;
436 } 436 }
437 return sPowerPerVolume; 437 return sPowerPerVolume;
438 } 438 }
439 439
440 //....oooOO0OOooo........oooOO0OOooo........oo 440 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
441 441
442 G4double G4PenelopeIonisationModel::MinEnergyC 442 G4double G4PenelopeIonisationModel::MinEnergyCut(const G4ParticleDefinition*,
443 const G4MaterialCutsCouple*) 443 const G4MaterialCutsCouple*)
444 { 444 {
445 return fIntrinsicLowEnergyLimit; 445 return fIntrinsicLowEnergyLimit;
446 } 446 }
447 447
448 //....oooOO0OOooo........oooOO0OOooo........oo 448 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
449 449
450 void G4PenelopeIonisationModel::SampleSecondar 450 void G4PenelopeIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
451 const G4MaterialCutsCouple* coup 451 const G4MaterialCutsCouple* couple,
452 const G4DynamicParticle* aDynami 452 const G4DynamicParticle* aDynamicParticle,
453 G4double cutE, G4double) 453 G4double cutE, G4double)
454 { 454 {
455 // Penelope model v2008 to sample the final 455 // Penelope model v2008 to sample the final state following an hard inelastic interaction.
456 // It makes use of the Generalised Oscillato 456 // It makes use of the Generalised Oscillator Strength (GOS) model from
457 // D. Liljequist, J. Phys. D: Appl. Phys. 1 457 // D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567
458 // 458 //
459 // The GOS model is used to calculate the in 459 // The GOS model is used to calculate the individual cross sections for all
460 // the atomic oscillators coming in the play 460 // the atomic oscillators coming in the play, taking into account the three
461 // contributions (distant longitudinal colli 461 // contributions (distant longitudinal collisions, distant transverse collisions and
462 // close collisions). Then the target shell 462 // close collisions). Then the target shell and the interaction channel are
463 // sampled. Final state of the delta-ray (en 463 // sampled. Final state of the delta-ray (energy, angle) are generated according
464 // to the analytical distributions (dSigma/d 464 // to the analytical distributions (dSigma/dW) for the selected interaction
465 // channels. 465 // channels.
466 // For e-, the maximum energy for the delta- 466 // For e-, the maximum energy for the delta-ray is initialEnergy/2. (because
467 // particles are indistinghusbable), while i 467 // particles are indistinghusbable), while it is the full initialEnergy for
468 // e+. 468 // e+.
469 // The efficiency of the random sampling alg 469 // The efficiency of the random sampling algorithm (e.g. for close collisions)
470 // decreases when initial and cutoff energy 470 // decreases when initial and cutoff energy increase (e.g. 87% for 10-keV primary
471 // and 1 keV threshold, 99% for 10-MeV prima 471 // and 1 keV threshold, 99% for 10-MeV primary and 10-keV threshold).
472 // Differential cross sections have a differ 472 // Differential cross sections have a different form for e+ and e-.
473 // 473 //
474 // WARNING: The model provides an _average_ 474 // WARNING: The model provides an _average_ description of inelastic collisions.
475 // Anyway, the energy spectrum associated to 475 // Anyway, the energy spectrum associated to distant excitations of a given
476 // atomic shell is approximated as a single 476 // atomic shell is approximated as a single resonance. The simulated energy spectra
477 // show _unphysical_ narrow peaks at energie 477 // show _unphysical_ narrow peaks at energies that are multiple of the shell
478 // resonance energies. The spurious speaks a 478 // resonance energies. The spurious speaks are automatically smoothed out after
479 // multiple inelastic collisions. 479 // multiple inelastic collisions.
480 // 480 //
481 // The model determines also the original sh 481 // The model determines also the original shell from which the delta-ray is expelled,
482 // in order to produce fluorescence de-excit 482 // in order to produce fluorescence de-excitation (from G4DeexcitationManager)
483 // 483 //
484 // Fermi density correction is calculated an 484 // Fermi density correction is calculated analytically according to
485 // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963), 485 // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1
486 486
487 if (fVerboseLevel > 3) 487 if (fVerboseLevel > 3)
488 G4cout << "Calling SamplingSecondaries() o 488 G4cout << "Calling SamplingSecondaries() of G4PenelopeIonisationModel" << G4endl;
489 489
490 G4double kineticEnergy0 = aDynamicParticle-> 490 G4double kineticEnergy0 = aDynamicParticle->GetKineticEnergy();
491 const G4ParticleDefinition* theParticle = aD 491 const G4ParticleDefinition* theParticle = aDynamicParticle->GetDefinition();
492 492
493 if (kineticEnergy0 <= fIntrinsicLowEnergyLim 493 if (kineticEnergy0 <= fIntrinsicLowEnergyLimit)
494 { 494 {
495 fParticleChange->SetProposedKineticEnergy( 495 fParticleChange->SetProposedKineticEnergy(0.);
496 fParticleChange->ProposeLocalEnergyDeposit 496 fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy0);
497 return ; 497 return ;
498 } 498 }
499 499
500 const G4Material* material = couple->GetMate 500 const G4Material* material = couple->GetMaterial();
501 const G4PenelopeOscillatorTable* theTable = 501 const G4PenelopeOscillatorTable* theTable = fOscManager->GetOscillatorTableIonisation(material);
502 502
503 G4ParticleMomentum particleDirection0 = aDyn 503 G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection();
504 504
505 //Initialise final-state variables. The prop 505 //Initialise final-state variables. The proper values will be set by the methods
506 // SampleFinalStateElectron() and SampleFina 506 // SampleFinalStateElectron() and SampleFinalStatePositron()
507 fKineticEnergy1=kineticEnergy0; 507 fKineticEnergy1=kineticEnergy0;
508 fCosThetaPrimary=1.0; 508 fCosThetaPrimary=1.0;
509 fEnergySecondary=0.0; 509 fEnergySecondary=0.0;
510 fCosThetaSecondary=1.0; 510 fCosThetaSecondary=1.0;
511 fTargetOscillator = -1; 511 fTargetOscillator = -1;
512 512
513 if (theParticle == G4Electron::Electron()) 513 if (theParticle == G4Electron::Electron())
514 SampleFinalStateElectron(material,cutE,kin 514 SampleFinalStateElectron(material,cutE,kineticEnergy0);
515 else if (theParticle == G4Positron::Positron 515 else if (theParticle == G4Positron::Positron())
516 SampleFinalStatePositron(material,cutE,kin 516 SampleFinalStatePositron(material,cutE,kineticEnergy0);
517 else 517 else
518 { 518 {
519 G4ExceptionDescription ed; 519 G4ExceptionDescription ed;
520 ed << "Invalid particle " << theParticle 520 ed << "Invalid particle " << theParticle->GetParticleName() << G4endl;
521 G4Exception("G4PenelopeIonisationModel:: 521 G4Exception("G4PenelopeIonisationModel::SamplingSecondaries()",
522 "em0001",FatalException,ed); 522 "em0001",FatalException,ed);
523 523
524 } 524 }
525 if (fEnergySecondary == 0) return; 525 if (fEnergySecondary == 0) return;
526 526
527 if (fVerboseLevel > 3) 527 if (fVerboseLevel > 3)
528 { 528 {
529 G4cout << "G4PenelopeIonisationModel::Sam 529 G4cout << "G4PenelopeIonisationModel::SamplingSecondaries() for " <<
530 theParticle->GetParticleName() << G4endl; 530 theParticle->GetParticleName() << G4endl;
531 G4cout << "Final eKin = " << fKineticEne 531 G4cout << "Final eKin = " << fKineticEnergy1 << " keV" << G4endl;
532 G4cout << "Final cosTheta = " << fCosThe 532 G4cout << "Final cosTheta = " << fCosThetaPrimary << G4endl;
533 G4cout << "Delta-ray eKin = " << fEnergy 533 G4cout << "Delta-ray eKin = " << fEnergySecondary << " keV" << G4endl;
534 G4cout << "Delta-ray cosTheta = " << fCo 534 G4cout << "Delta-ray cosTheta = " << fCosThetaSecondary << G4endl;
535 G4cout << "Oscillator: " << fTargetOscil 535 G4cout << "Oscillator: " << fTargetOscillator << G4endl;
536 } 536 }
537 537
538 //Update the primary particle 538 //Update the primary particle
539 G4double sint = std::sqrt(1. - fCosThetaPrim 539 G4double sint = std::sqrt(1. - fCosThetaPrimary*fCosThetaPrimary);
540 G4double phiPrimary = twopi * G4UniformRand 540 G4double phiPrimary = twopi * G4UniformRand();
541 G4double dirx = sint * std::cos(phiPrimary); 541 G4double dirx = sint * std::cos(phiPrimary);
542 G4double diry = sint * std::sin(phiPrimary); 542 G4double diry = sint * std::sin(phiPrimary);
543 G4double dirz = fCosThetaPrimary; 543 G4double dirz = fCosThetaPrimary;
544 544
545 G4ThreeVector electronDirection1(dirx,diry,d 545 G4ThreeVector electronDirection1(dirx,diry,dirz);
546 electronDirection1.rotateUz(particleDirectio 546 electronDirection1.rotateUz(particleDirection0);
547 547
548 if (fKineticEnergy1 > 0) 548 if (fKineticEnergy1 > 0)
549 { 549 {
550 fParticleChange->ProposeMomentumDirectio 550 fParticleChange->ProposeMomentumDirection(electronDirection1);
551 fParticleChange->SetProposedKineticEnerg 551 fParticleChange->SetProposedKineticEnergy(fKineticEnergy1);
552 } 552 }
553 else 553 else
554 fParticleChange->SetProposedKineticEnergy( 554 fParticleChange->SetProposedKineticEnergy(0.);
555 555
556 //Generate the delta ray 556 //Generate the delta ray
557 G4double ionEnergyInPenelopeDatabase = 557 G4double ionEnergyInPenelopeDatabase =
558 (*theTable)[fTargetOscillator]->GetIonisat 558 (*theTable)[fTargetOscillator]->GetIonisationEnergy();
559 559
560 //Now, try to handle fluorescence 560 //Now, try to handle fluorescence
561 //Notice: merged levels are indicated with Z 561 //Notice: merged levels are indicated with Z=0 and flag=30
562 G4int shFlag = (*theTable)[fTargetOscillator 562 G4int shFlag = (*theTable)[fTargetOscillator]->GetShellFlag();
563 G4int Z = (G4int) (*theTable)[fTargetOscilla 563 G4int Z = (G4int) (*theTable)[fTargetOscillator]->GetParentZ();
564 564
565 //initialize here, then check photons create 565 //initialize here, then check photons created by Atomic-Deexcitation, and the final state e-
566 const G4AtomicTransitionManager* transitionM 566 const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
567 G4double bindingEnergy = 0.*eV; 567 G4double bindingEnergy = 0.*eV;
568 568
569 const G4AtomicShell* shell = nullptr; 569 const G4AtomicShell* shell = nullptr;
570 //Real level 570 //Real level
571 if (Z > 0 && shFlag<30) 571 if (Z > 0 && shFlag<30)
572 { 572 {
573 shell = transitionManager->Shell(Z,shFla 573 shell = transitionManager->Shell(Z,shFlag-1);
574 bindingEnergy = shell->BindingEnergy(); 574 bindingEnergy = shell->BindingEnergy();
575 //shellId = shell->ShellId(); 575 //shellId = shell->ShellId();
576 } 576 }
577 577
578 //correct the fEnergySecondary to account fo 578 //correct the fEnergySecondary to account for the fact that the Penelope
579 //database of ionisation energies is in gene 579 //database of ionisation energies is in general (slightly) different
580 //from the fluorescence database used in Gea 580 //from the fluorescence database used in Geant4.
581 fEnergySecondary += ionEnergyInPenelopeDatab 581 fEnergySecondary += ionEnergyInPenelopeDatabase-bindingEnergy;
582 582
583 G4double localEnergyDeposit = bindingEnergy; 583 G4double localEnergyDeposit = bindingEnergy;
584 //testing purposes only 584 //testing purposes only
585 G4double energyInFluorescence = 0; 585 G4double energyInFluorescence = 0;
586 G4double energyInAuger = 0; 586 G4double energyInAuger = 0;
587 587
588 if (fEnergySecondary < 0) 588 if (fEnergySecondary < 0)
589 { 589 {
590 //It means that there was some problem/m 590 //It means that there was some problem/mismatch between the two databases.
591 //In this case, the available energy is 591 //In this case, the available energy is ok to excite the level according
592 //to the Penelope database, but not acco 592 //to the Penelope database, but not according to the Geant4 database
593 //Full residual energy is deposited loca 593 //Full residual energy is deposited locally
594 localEnergyDeposit += fEnergySecondary; 594 localEnergyDeposit += fEnergySecondary;
595 fEnergySecondary = 0.0; 595 fEnergySecondary = 0.0;
596 } 596 }
597 597
598 //Notice: shell might be nullptr (invalid!) 598 //Notice: shell might be nullptr (invalid!) if shFlag=30. Must be protected
599 //Disable the built-in de-excitation of the 599 //Disable the built-in de-excitation of the PIXE flag is active. In this
600 //case, the PIXE interface takes care (stati 600 //case, the PIXE interface takes care (statistically) of producing the
601 //de-excitation. 601 //de-excitation.
602 //Now, take care of fluorescence, if require 602 //Now, take care of fluorescence, if required
603 if (fAtomDeexcitation && !fPIXEflag && shell 603 if (fAtomDeexcitation && !fPIXEflag && shell)
604 { 604 {
605 G4int index = couple->GetIndex(); 605 G4int index = couple->GetIndex();
606 if (fAtomDeexcitation->CheckDeexcitation 606 if (fAtomDeexcitation->CheckDeexcitationActiveRegion(index))
607 { 607 {
608 std::size_t nBefore = fvect->size(); 608 std::size_t nBefore = fvect->size();
609 fAtomDeexcitation->GenerateParticles(fvect 609 fAtomDeexcitation->GenerateParticles(fvect,shell,Z,index);
610 std::size_t nAfter = fvect->size(); 610 std::size_t nAfter = fvect->size();
611 611
612 if (nAfter>nBefore) //actual production of 612 if (nAfter>nBefore) //actual production of fluorescence
613 { 613 {
614 for (std::size_t j=nBefore;j<nAfter;++ 614 for (std::size_t j=nBefore;j<nAfter;++j) //loop on products
615 { 615 {
616 G4double itsEnergy = ((*fvect)[j])->GetK 616 G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy();
617 if (itsEnergy < localEnergyD 617 if (itsEnergy < localEnergyDeposit) // valid secondary, generate it
618 { 618 {
619 localEnergyDeposit -= itsEnergy; 619 localEnergyDeposit -= itsEnergy;
620 if (((*fvect)[j])->GetParticleDefini 620 if (((*fvect)[j])->GetParticleDefinition() == G4Gamma::Definition())
621 energyInFluorescence += itsEnergy; 621 energyInFluorescence += itsEnergy;
622 else if (((*fvect)[j])->GetParticleD 622 else if (((*fvect)[j])->GetParticleDefinition() == G4Electron::Definition())
623 energyInAuger += itsEnergy; 623 energyInAuger += itsEnergy;
624 } 624 }
625 else //invalid secondary: ta 625 else //invalid secondary: takes more than the available energy: delete it
626 { 626 {
627 delete (*fvect)[j]; 627 delete (*fvect)[j];
628 (*fvect)[j] = nullptr; 628 (*fvect)[j] = nullptr;
629 } 629 }
630 } 630 }
631 } 631 }
632 } 632 }
633 } 633 }
634 634
635 // Generate the delta ray --> to be done onl 635 // Generate the delta ray --> to be done only if above cut
636 if (fEnergySecondary > cutE) 636 if (fEnergySecondary > cutE)
637 { 637 {
638 G4DynamicParticle* electron = nullptr; 638 G4DynamicParticle* electron = nullptr;
639 G4double sinThetaE = std::sqrt(1.-fCosTh 639 G4double sinThetaE = std::sqrt(1.-fCosThetaSecondary*fCosThetaSecondary);
640 G4double phiEl = phiPrimary+pi; //pi wit 640 G4double phiEl = phiPrimary+pi; //pi with respect to the primary electron/positron
641 G4double xEl = sinThetaE * std::cos(phiE 641 G4double xEl = sinThetaE * std::cos(phiEl);
642 G4double yEl = sinThetaE * std::sin(phiE 642 G4double yEl = sinThetaE * std::sin(phiEl);
643 G4double zEl = fCosThetaSecondary; 643 G4double zEl = fCosThetaSecondary;
644 G4ThreeVector eDirection(xEl,yEl,zEl); / 644 G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction
645 eDirection.rotateUz(particleDirection0); 645 eDirection.rotateUz(particleDirection0);
646 electron = new G4DynamicParticle (G4Elec 646 electron = new G4DynamicParticle (G4Electron::Electron(),
647 eDirection,fEnergySecondary) ; 647 eDirection,fEnergySecondary) ;
648 fvect->push_back(electron); 648 fvect->push_back(electron);
649 } 649 }
650 else 650 else
651 { 651 {
652 localEnergyDeposit += fEnergySecondary; 652 localEnergyDeposit += fEnergySecondary;
653 fEnergySecondary = 0; 653 fEnergySecondary = 0;
654 } 654 }
655 655
656 if (localEnergyDeposit < 0) //Should not be: 656 if (localEnergyDeposit < 0) //Should not be: issue a G4Exception (warning)
657 { 657 {
658 G4Exception("G4PenelopeIonisationModel:: 658 G4Exception("G4PenelopeIonisationModel::SampleSecondaries()",
659 "em2099",JustWarning,"WARNING: Negative 659 "em2099",JustWarning,"WARNING: Negative local energy deposit");
660 localEnergyDeposit=0.; 660 localEnergyDeposit=0.;
661 } 661 }
662 fParticleChange->ProposeLocalEnergyDeposit(l 662 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
663 663
664 if (fVerboseLevel > 1) 664 if (fVerboseLevel > 1)
665 { 665 {
666 G4cout << "----------------------------- 666 G4cout << "-----------------------------------------------------------" << G4endl;
667 G4cout << "Energy balance from G4Penelop 667 G4cout << "Energy balance from G4PenelopeIonisation" << G4endl;
668 G4cout << "Incoming primary energy: " << 668 G4cout << "Incoming primary energy: " << kineticEnergy0/keV << " keV" << G4endl;
669 G4cout << "----------------------------- 669 G4cout << "-----------------------------------------------------------" << G4endl;
670 G4cout << "Outgoing primary energy: " << 670 G4cout << "Outgoing primary energy: " << fKineticEnergy1/keV << " keV" << G4endl;
671 G4cout << "Delta ray " << fEnergySeconda 671 G4cout << "Delta ray " << fEnergySecondary/keV << " keV" << G4endl;
672 if (energyInFluorescence) 672 if (energyInFluorescence)
673 G4cout << "Fluorescence x-rays: " << energyI 673 G4cout << "Fluorescence x-rays: " << energyInFluorescence/keV << " keV" << G4endl;
674 if (energyInAuger) 674 if (energyInAuger)
675 G4cout << "Auger electrons: " << energyInAug 675 G4cout << "Auger electrons: " << energyInAuger/keV << " keV" << G4endl;
676 G4cout << "Local energy deposit " << loc 676 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
677 G4cout << "Total final state: " << (fEne 677 G4cout << "Total final state: " << (fEnergySecondary+energyInFluorescence+fKineticEnergy1+
678 localEnergyDeposit+energyInAuger)/ 678 localEnergyDeposit+energyInAuger)/keV <<
679 " keV" << G4endl; 679 " keV" << G4endl;
680 G4cout << "----------------------------- 680 G4cout << "-----------------------------------------------------------" << G4endl;
681 } 681 }
682 682
683 if (fVerboseLevel > 0) 683 if (fVerboseLevel > 0)
684 { 684 {
685 G4double energyDiff = std::fabs(fEnergyS 685 G4double energyDiff = std::fabs(fEnergySecondary+energyInFluorescence+fKineticEnergy1+
686 localEnergyDeposit+energyInAuger 686 localEnergyDeposit+energyInAuger-kineticEnergy0);
687 if (energyDiff > 0.05*keV) 687 if (energyDiff > 0.05*keV)
688 G4cout << "Warning from G4PenelopeIonisation 688 G4cout << "Warning from G4PenelopeIonisation: problem with energy conservation: " <<
689 (fEnergySecondary+energyInFluorescence+fKi 689 (fEnergySecondary+energyInFluorescence+fKineticEnergy1+localEnergyDeposit+energyInAuger)/keV <<
690 " keV (final) vs. " << 690 " keV (final) vs. " <<
691 kineticEnergy0/keV << " keV (initial)" << 691 kineticEnergy0/keV << " keV (initial)" << G4endl;
692 } 692 }
693 693
694 } 694 }
695 695
696 //....oooOO0OOooo........oooOO0OOooo........oo 696 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
697 void G4PenelopeIonisationModel::SampleFinalSta 697 void G4PenelopeIonisationModel::SampleFinalStateElectron(const G4Material* mat,
698 G4double cutEnergy, 698 G4double cutEnergy,
699 G4double kineticEnergy) 699 G4double kineticEnergy)
700 { 700 {
701 // This method sets the final ionisation par 701 // This method sets the final ionisation parameters
702 // fKineticEnergy1, fCosThetaPrimary (= upda 702 // fKineticEnergy1, fCosThetaPrimary (= updates of the primary e-)
703 // fEnergySecondary, fCosThetaSecondary (= i 703 // fEnergySecondary, fCosThetaSecondary (= info of delta-ray)
704 // fTargetOscillator (= ionised oscillator) 704 // fTargetOscillator (= ionised oscillator)
705 // 705 //
706 // The method implements SUBROUTINE EINa of 706 // The method implements SUBROUTINE EINa of Penelope
707 // 707 //
708 708
709 G4PenelopeOscillatorTable* theTable = fOscMa 709 G4PenelopeOscillatorTable* theTable = fOscManager->GetOscillatorTableIonisation(mat);
710 std::size_t numberOfOscillators = theTable-> 710 std::size_t numberOfOscillators = theTable->size();
711 const G4PenelopeCrossSection* theXS = 711 const G4PenelopeCrossSection* theXS =
712 fCrossSectionHandler->GetCrossSectionTable 712 fCrossSectionHandler->GetCrossSectionTableForCouple(G4Electron::Electron(),mat,
713 cutEnergy); 713 cutEnergy);
714 G4double delta = fCrossSectionHandler->GetDe 714 G4double delta = fCrossSectionHandler->GetDensityCorrection(mat,kineticEnergy);
715 715
716 // Selection of the active oscillator 716 // Selection of the active oscillator
717 G4double TST = G4UniformRand(); 717 G4double TST = G4UniformRand();
718 fTargetOscillator = G4int(numberOfOscillator 718 fTargetOscillator = G4int(numberOfOscillators-1); //initialization, last oscillator
719 G4double XSsum = 0.; 719 G4double XSsum = 0.;
720 720
721 for (std::size_t i=0;i<numberOfOscillators-1 721 for (std::size_t i=0;i<numberOfOscillators-1;++i)
722 { 722 {
723 XSsum += theXS->GetNormalizedShellCrossS 723 XSsum += theXS->GetNormalizedShellCrossSection(i,kineticEnergy);
724 724
725 if (XSsum > TST) 725 if (XSsum > TST)
726 { 726 {
727 fTargetOscillator = (G4int) i; 727 fTargetOscillator = (G4int) i;
728 break; 728 break;
729 } 729 }
730 } 730 }
731 731
732 if (fVerboseLevel > 3) 732 if (fVerboseLevel > 3)
733 { 733 {
734 G4cout << "SampleFinalStateElectron: sam 734 G4cout << "SampleFinalStateElectron: sampled oscillator #" <<
735 fTargetOscillator << "." << G4endl; 735 fTargetOscillator << "." << G4endl;
736 G4cout << "Ionisation energy: " << 736 G4cout << "Ionisation energy: " <<
737 (*theTable)[fTargetOscillator]->GetIonisatio 737 (*theTable)[fTargetOscillator]->GetIonisationEnergy()/eV <<
738 " eV " << G4endl; 738 " eV " << G4endl;
739 G4cout << "Resonance energy: : " << 739 G4cout << "Resonance energy: : " <<
740 (*theTable)[fTargetOscillator]->GetResonance 740 (*theTable)[fTargetOscillator]->GetResonanceEnergy()/eV << " eV "
741 << G4endl; 741 << G4endl;
742 } 742 }
743 //Constants 743 //Constants
744 G4double rb = kineticEnergy + 2.0*electron_m 744 G4double rb = kineticEnergy + 2.0*electron_mass_c2;
745 G4double gam = 1.0+kineticEnergy/electron_ma 745 G4double gam = 1.0+kineticEnergy/electron_mass_c2;
746 G4double gam2 = gam*gam; 746 G4double gam2 = gam*gam;
747 G4double beta2 = (gam2-1.0)/gam2; 747 G4double beta2 = (gam2-1.0)/gam2;
748 G4double amol = ((gam-1.0)/gam)*((gam-1.0)/g 748 G4double amol = ((gam-1.0)/gam)*((gam-1.0)/gam);
749 749
750 //Partial cross section of the active oscill 750 //Partial cross section of the active oscillator
751 G4double resEne = (*theTable)[fTargetOscilla 751 G4double resEne = (*theTable)[fTargetOscillator]->GetResonanceEnergy();
752 G4double invResEne = 1.0/resEne; 752 G4double invResEne = 1.0/resEne;
753 G4double ionEne = (*theTable)[fTargetOscilla 753 G4double ionEne = (*theTable)[fTargetOscillator]->GetIonisationEnergy();
754 G4double cutoffEne = (*theTable)[fTargetOsci 754 G4double cutoffEne = (*theTable)[fTargetOscillator]->GetCutoffRecoilResonantEnergy();
755 G4double XHDL = 0.; 755 G4double XHDL = 0.;
756 G4double XHDT = 0.; 756 G4double XHDT = 0.;
757 G4double QM = 0.; 757 G4double QM = 0.;
758 G4double cps = 0.; 758 G4double cps = 0.;
759 G4double cp = 0.; 759 G4double cp = 0.;
760 760
761 //Distant excitations 761 //Distant excitations
762 if (resEne > cutEnergy && resEne < kineticEn 762 if (resEne > cutEnergy && resEne < kineticEnergy)
763 { 763 {
764 cps = kineticEnergy*rb; 764 cps = kineticEnergy*rb;
765 cp = std::sqrt(cps); 765 cp = std::sqrt(cps);
766 G4double XHDT0 = std::max(G4Log(gam2)-be 766 G4double XHDT0 = std::max(G4Log(gam2)-beta2-delta,0.);
767 if (resEne > 1.0e-6*kineticEnergy) 767 if (resEne > 1.0e-6*kineticEnergy)
768 { 768 {
769 G4double cpp = std::sqrt((kineticEnergy-re 769 G4double cpp = std::sqrt((kineticEnergy-resEne)*(kineticEnergy-resEne+2.0*electron_mass_c2));
770 QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_ 770 QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
771 } 771 }
772 else 772 else
773 { 773 {
774 QM = resEne*resEne/(beta2*2.0*electron_mas 774 QM = resEne*resEne/(beta2*2.0*electron_mass_c2);
775 QM *= (1.0-QM*0.5/electron_mass_c2); 775 QM *= (1.0-QM*0.5/electron_mass_c2);
776 } 776 }
777 if (QM < cutoffEne) 777 if (QM < cutoffEne)
778 { 778 {
779 XHDL = G4Log(cutoffEne*(QM+2.0*electron_ma 779 XHDL = G4Log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)))
780 *invResEne; 780 *invResEne;
781 XHDT = XHDT0*invResEne; 781 XHDT = XHDT0*invResEne;
782 } 782 }
783 else 783 else
784 { 784 {
785 QM = cutoffEne; 785 QM = cutoffEne;
786 XHDL = 0.; 786 XHDL = 0.;
787 XHDT = 0.; 787 XHDT = 0.;
788 } 788 }
789 } 789 }
790 else 790 else
791 { 791 {
792 QM = cutoffEne; 792 QM = cutoffEne;
793 cps = 0.; 793 cps = 0.;
794 cp = 0.; 794 cp = 0.;
795 XHDL = 0.; 795 XHDL = 0.;
796 XHDT = 0.; 796 XHDT = 0.;
797 } 797 }
798 798
799 //Close collisions 799 //Close collisions
800 G4double EE = kineticEnergy + ionEne; 800 G4double EE = kineticEnergy + ionEne;
801 G4double wmaxc = 0.5*EE; 801 G4double wmaxc = 0.5*EE;
802 G4double wcl = std::max(cutEnergy,cutoffEne) 802 G4double wcl = std::max(cutEnergy,cutoffEne);
803 G4double rcl = wcl/EE; 803 G4double rcl = wcl/EE;
804 G4double XHC = 0.; 804 G4double XHC = 0.;
805 if (wcl < wmaxc) 805 if (wcl < wmaxc)
806 { 806 {
807 G4double rl1 = 1.0-rcl; 807 G4double rl1 = 1.0-rcl;
808 G4double rrl1 = 1.0/rl1; 808 G4double rrl1 = 1.0/rl1;
809 XHC = (amol*(0.5-rcl)+1.0/rcl-rrl1+ 809 XHC = (amol*(0.5-rcl)+1.0/rcl-rrl1+
810 (1.0-amol)*G4Log(rcl*rrl1))/EE; 810 (1.0-amol)*G4Log(rcl*rrl1))/EE;
811 } 811 }
812 812
813 //Total cross section per molecule for the a 813 //Total cross section per molecule for the active shell, in cm2
814 G4double XHTOT = XHC + XHDL + XHDT; 814 G4double XHTOT = XHC + XHDL + XHDT;
815 815
816 //very small cross section, do nothing 816 //very small cross section, do nothing
817 if (XHTOT < 1.e-14*barn) 817 if (XHTOT < 1.e-14*barn)
818 { 818 {
819 fKineticEnergy1=kineticEnergy; 819 fKineticEnergy1=kineticEnergy;
820 fCosThetaPrimary=1.0; 820 fCosThetaPrimary=1.0;
821 fEnergySecondary=0.0; 821 fEnergySecondary=0.0;
822 fCosThetaSecondary=1.0; 822 fCosThetaSecondary=1.0;
823 fTargetOscillator = G4int(numberOfOscill 823 fTargetOscillator = G4int(numberOfOscillators-1);
824 return; 824 return;
825 } 825 }
826 826
827 //decide which kind of interaction we'll hav 827 //decide which kind of interaction we'll have
828 TST = XHTOT*G4UniformRand(); 828 TST = XHTOT*G4UniformRand();
829 829
830 // Hard close collision 830 // Hard close collision
831 G4double TS1 = XHC; 831 G4double TS1 = XHC;
832 832
833 if (TST < TS1) 833 if (TST < TS1)
834 { 834 {
835 G4double A = 5.0*amol; 835 G4double A = 5.0*amol;
836 G4double ARCL = A*0.5*rcl; 836 G4double ARCL = A*0.5*rcl;
837 G4double rk=0.; 837 G4double rk=0.;
838 G4bool loopAgain = false; 838 G4bool loopAgain = false;
839 do 839 do
840 { 840 {
841 loopAgain = false; 841 loopAgain = false;
842 G4double fb = (1.0+ARCL)*G4UniformRand(); 842 G4double fb = (1.0+ARCL)*G4UniformRand();
843 if (fb < 1) 843 if (fb < 1)
844 rk = rcl/(1.0-fb*(1.0-(rcl+rcl))); 844 rk = rcl/(1.0-fb*(1.0-(rcl+rcl)));
845 else 845 else
846 rk = rcl + (fb-1.0)*(0.5-rcl)/ARCL; 846 rk = rcl + (fb-1.0)*(0.5-rcl)/ARCL;
847 G4double rk2 = rk*rk; 847 G4double rk2 = rk*rk;
848 G4double rkf = rk/(1.0-rk); 848 G4double rkf = rk/(1.0-rk);
849 G4double phi = 1.0+rkf*rkf-rkf+amol*(rk2+r 849 G4double phi = 1.0+rkf*rkf-rkf+amol*(rk2+rkf);
850 if (G4UniformRand()*(1.0+A*rk2) > phi) 850 if (G4UniformRand()*(1.0+A*rk2) > phi)
851 loopAgain = true; 851 loopAgain = true;
852 }while(loopAgain); 852 }while(loopAgain);
853 //energy and scattering angle (primary e 853 //energy and scattering angle (primary electron)
854 G4double deltaE = rk*EE; 854 G4double deltaE = rk*EE;
855 855
856 fKineticEnergy1 = kineticEnergy - deltaE 856 fKineticEnergy1 = kineticEnergy - deltaE;
857 fCosThetaPrimary = std::sqrt(fKineticEne 857 fCosThetaPrimary = std::sqrt(fKineticEnergy1*rb/(kineticEnergy*(rb-deltaE)));
858 //energy and scattering angle of the del 858 //energy and scattering angle of the delta ray
859 fEnergySecondary = deltaE - ionEne; //su 859 fEnergySecondary = deltaE - ionEne; //subtract ionisation energy
860 fCosThetaSecondary= std::sqrt(deltaE*rb/ 860 fCosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*electron_mass_c2)));
861 if (fVerboseLevel > 3) 861 if (fVerboseLevel > 3)
862 G4cout << "SampleFinalStateElectron: sampled 862 G4cout << "SampleFinalStateElectron: sampled close collision " << G4endl;
863 return; 863 return;
864 } 864 }
865 865
866 //Hard distant longitudinal collisions 866 //Hard distant longitudinal collisions
867 TS1 += XHDL; 867 TS1 += XHDL;
868 G4double deltaE = resEne; 868 G4double deltaE = resEne;
869 fKineticEnergy1 = kineticEnergy - deltaE; 869 fKineticEnergy1 = kineticEnergy - deltaE;
870 870
871 if (TST < TS1) 871 if (TST < TS1)
872 { 872 {
873 G4double QS = QM/(1.0+QM*0.5/electron_ma 873 G4double QS = QM/(1.0+QM*0.5/electron_mass_c2);
874 G4double Q = QS/(std::pow((QS/cutoffEne) 874 G4double Q = QS/(std::pow((QS/cutoffEne)*(1.0+cutoffEne*0.5/electron_mass_c2),G4UniformRand())
875 - (QS*0.5/electron_mass_c2)); 875 - (QS*0.5/electron_mass_c2));
876 G4double QTREV = Q*(Q+2.0*electron_mass_ 876 G4double QTREV = Q*(Q+2.0*electron_mass_c2);
877 G4double cpps = fKineticEnergy1*(fKineti 877 G4double cpps = fKineticEnergy1*(fKineticEnergy1+2.0*electron_mass_c2);
878 fCosThetaPrimary = (cpps+cps-QTREV)/(2.0 878 fCosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps));
879 if (fCosThetaPrimary > 1.) 879 if (fCosThetaPrimary > 1.)
880 fCosThetaPrimary = 1.0; 880 fCosThetaPrimary = 1.0;
881 //energy and emission angle of the delta 881 //energy and emission angle of the delta ray
882 fEnergySecondary = deltaE - ionEne; 882 fEnergySecondary = deltaE - ionEne;
883 fCosThetaSecondary = 0.5*(deltaE*(kineti 883 fCosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
884 if (fCosThetaSecondary > 1.0) 884 if (fCosThetaSecondary > 1.0)
885 fCosThetaSecondary = 1.0; 885 fCosThetaSecondary = 1.0;
886 if (fVerboseLevel > 3) 886 if (fVerboseLevel > 3)
887 G4cout << "SampleFinalStateElectron: sampled 887 G4cout << "SampleFinalStateElectron: sampled distant longitudinal collision " << G4endl;
888 return; 888 return;
889 } 889 }
890 890
891 //Hard distant transverse collisions 891 //Hard distant transverse collisions
892 fCosThetaPrimary = 1.0; 892 fCosThetaPrimary = 1.0;
893 //energy and emission angle of the delta ray 893 //energy and emission angle of the delta ray
894 fEnergySecondary = deltaE - ionEne; 894 fEnergySecondary = deltaE - ionEne;
895 fCosThetaSecondary = 0.5; 895 fCosThetaSecondary = 0.5;
896 if (fVerboseLevel > 3) 896 if (fVerboseLevel > 3)
897 G4cout << "SampleFinalStateElectron: sampl 897 G4cout << "SampleFinalStateElectron: sampled distant transverse collision " << G4endl;
898 898
899 return; 899 return;
900 } 900 }
901 901
902 //....oooOO0OOooo........oooOO0OOooo........oo 902 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
903 903
904 void G4PenelopeIonisationModel::SampleFinalSta 904 void G4PenelopeIonisationModel::SampleFinalStatePositron(const G4Material* mat,
905 G4double cutEnergy, 905 G4double cutEnergy,
906 G4double kineticEnergy) 906 G4double kineticEnergy)
907 { 907 {
908 // This method sets the final ionisation par 908 // This method sets the final ionisation parameters
909 // fKineticEnergy1, fCosThetaPrimary (= upda 909 // fKineticEnergy1, fCosThetaPrimary (= updates of the primary e-)
910 // fEnergySecondary, fCosThetaSecondary (= i 910 // fEnergySecondary, fCosThetaSecondary (= info of delta-ray)
911 // fTargetOscillator (= ionised oscillator) 911 // fTargetOscillator (= ionised oscillator)
912 // 912 //
913 // The method implements SUBROUTINE PINa of 913 // The method implements SUBROUTINE PINa of Penelope
914 // 914 //
915 915
916 G4PenelopeOscillatorTable* theTable = fOscMa 916 G4PenelopeOscillatorTable* theTable = fOscManager->GetOscillatorTableIonisation(mat);
917 std::size_t numberOfOscillators = theTable-> 917 std::size_t numberOfOscillators = theTable->size();
918 const G4PenelopeCrossSection* theXS = 918 const G4PenelopeCrossSection* theXS =
919 fCrossSectionHandler->GetCrossSectionTable 919 fCrossSectionHandler->GetCrossSectionTableForCouple(G4Positron::Positron(),mat,
920 cutEnergy); 920 cutEnergy);
921 G4double delta = fCrossSectionHandler->GetDe 921 G4double delta = fCrossSectionHandler->GetDensityCorrection(mat,kineticEnergy);
922 922
923 // Selection of the active oscillator 923 // Selection of the active oscillator
924 G4double TST = G4UniformRand(); 924 G4double TST = G4UniformRand();
925 fTargetOscillator = G4int(numberOfOscillator 925 fTargetOscillator = G4int(numberOfOscillators-1); //initialization, last oscillator
926 G4double XSsum = 0.; 926 G4double XSsum = 0.;
927 for (std::size_t i=0;i<numberOfOscillators-1 927 for (std::size_t i=0;i<numberOfOscillators-1;++i)
928 { 928 {
929 XSsum += theXS->GetNormalizedShellCrossS 929 XSsum += theXS->GetNormalizedShellCrossSection(i,kineticEnergy);
930 if (XSsum > TST) 930 if (XSsum > TST)
931 { 931 {
932 fTargetOscillator = (G4int) i; 932 fTargetOscillator = (G4int) i;
933 break; 933 break;
934 } 934 }
935 } 935 }
936 936
937 if (fVerboseLevel > 3) 937 if (fVerboseLevel > 3)
938 { 938 {
939 G4cout << "SampleFinalStatePositron: sam 939 G4cout << "SampleFinalStatePositron: sampled oscillator #" <<
940 fTargetOscillator << "." << G4endl; 940 fTargetOscillator << "." << G4endl;
941 G4cout << "Ionisation energy: " << (*the 941 G4cout << "Ionisation energy: " << (*theTable)[fTargetOscillator]->GetIonisationEnergy()/eV
942 << " eV " << G4endl; 942 << " eV " << G4endl;
943 G4cout << "Resonance energy: : " << (*th 943 G4cout << "Resonance energy: : " << (*theTable)[fTargetOscillator]->GetResonanceEnergy()/eV
944 << " eV " << G4endl; 944 << " eV " << G4endl;
945 } 945 }
946 946
947 //Constants 947 //Constants
948 G4double rb = kineticEnergy + 2.0*electron_m 948 G4double rb = kineticEnergy + 2.0*electron_mass_c2;
949 G4double gam = 1.0+kineticEnergy/electron_ma 949 G4double gam = 1.0+kineticEnergy/electron_mass_c2;
950 G4double gam2 = gam*gam; 950 G4double gam2 = gam*gam;
951 G4double beta2 = (gam2-1.0)/gam2; 951 G4double beta2 = (gam2-1.0)/gam2;
952 G4double g12 = (gam+1.0)*(gam+1.0); 952 G4double g12 = (gam+1.0)*(gam+1.0);
953 G4double amol = ((gam-1.0)/gam)*((gam-1.0)/g 953 G4double amol = ((gam-1.0)/gam)*((gam-1.0)/gam);
954 //Bhabha coefficients 954 //Bhabha coefficients
955 G4double bha1 = amol*(2.0*g12-1.0)/(gam2-1.0 955 G4double bha1 = amol*(2.0*g12-1.0)/(gam2-1.0);
956 G4double bha2 = amol*(3.0+1.0/g12); 956 G4double bha2 = amol*(3.0+1.0/g12);
957 G4double bha3 = amol*2.0*gam*(gam-1.0)/g12; 957 G4double bha3 = amol*2.0*gam*(gam-1.0)/g12;
958 G4double bha4 = amol*(gam-1.0)*(gam-1.0)/g12 958 G4double bha4 = amol*(gam-1.0)*(gam-1.0)/g12;
959 959
960 // 960 //
961 //Partial cross section of the active oscill 961 //Partial cross section of the active oscillator
962 // 962 //
963 G4double resEne = (*theTable)[fTargetOscilla 963 G4double resEne = (*theTable)[fTargetOscillator]->GetResonanceEnergy();
964 G4double invResEne = 1.0/resEne; 964 G4double invResEne = 1.0/resEne;
965 G4double ionEne = (*theTable)[fTargetOscilla 965 G4double ionEne = (*theTable)[fTargetOscillator]->GetIonisationEnergy();
966 G4double cutoffEne = (*theTable)[fTargetOsci 966 G4double cutoffEne = (*theTable)[fTargetOscillator]->GetCutoffRecoilResonantEnergy();
967 967
968 G4double XHDL = 0.; 968 G4double XHDL = 0.;
969 G4double XHDT = 0.; 969 G4double XHDT = 0.;
970 G4double QM = 0.; 970 G4double QM = 0.;
971 G4double cps = 0.; 971 G4double cps = 0.;
972 G4double cp = 0.; 972 G4double cp = 0.;
973 973
974 //Distant excitations XS (same as for electr 974 //Distant excitations XS (same as for electrons)
975 if (resEne > cutEnergy && resEne < kineticEn 975 if (resEne > cutEnergy && resEne < kineticEnergy)
976 { 976 {
977 cps = kineticEnergy*rb; 977 cps = kineticEnergy*rb;
978 cp = std::sqrt(cps); 978 cp = std::sqrt(cps);
979 G4double XHDT0 = std::max(G4Log(gam2)-be 979 G4double XHDT0 = std::max(G4Log(gam2)-beta2-delta,0.);
980 if (resEne > 1.0e-6*kineticEnergy) 980 if (resEne > 1.0e-6*kineticEnergy)
981 { 981 {
982 G4double cpp = std::sqrt((kineticEnergy-re 982 G4double cpp = std::sqrt((kineticEnergy-resEne)*(kineticEnergy-resEne+2.0*electron_mass_c2));
983 QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_ 983 QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
984 } 984 }
985 else 985 else
986 { 986 {
987 QM = resEne*resEne/(beta2*2.0*electron_mas 987 QM = resEne*resEne/(beta2*2.0*electron_mass_c2);
988 QM *= (1.0-QM*0.5/electron_mass_c2); 988 QM *= (1.0-QM*0.5/electron_mass_c2);
989 } 989 }
990 if (QM < cutoffEne) 990 if (QM < cutoffEne)
991 { 991 {
992 XHDL = G4Log(cutoffEne*(QM+2.0*electron_ma 992 XHDL = G4Log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)))
993 *invResEne; 993 *invResEne;
994 XHDT = XHDT0*invResEne; 994 XHDT = XHDT0*invResEne;
995 } 995 }
996 else 996 else
997 { 997 {
998 QM = cutoffEne; 998 QM = cutoffEne;
999 XHDL = 0.; 999 XHDL = 0.;
1000 XHDT = 0.; 1000 XHDT = 0.;
1001 } 1001 }
1002 } 1002 }
1003 else 1003 else
1004 { 1004 {
1005 QM = cutoffEne; 1005 QM = cutoffEne;
1006 cps = 0.; 1006 cps = 0.;
1007 cp = 0.; 1007 cp = 0.;
1008 XHDL = 0.; 1008 XHDL = 0.;
1009 XHDT = 0.; 1009 XHDT = 0.;
1010 } 1010 }
1011 //Close collisions (Bhabha) 1011 //Close collisions (Bhabha)
1012 G4double wmaxc = kineticEnergy; 1012 G4double wmaxc = kineticEnergy;
1013 G4double wcl = std::max(cutEnergy,cutoffEne 1013 G4double wcl = std::max(cutEnergy,cutoffEne);
1014 G4double rcl = wcl/kineticEnergy; 1014 G4double rcl = wcl/kineticEnergy;
1015 G4double XHC = 0.; 1015 G4double XHC = 0.;
1016 if (wcl < wmaxc) 1016 if (wcl < wmaxc)
1017 { 1017 {
1018 G4double rl1 = 1.0-rcl; 1018 G4double rl1 = 1.0-rcl;
1019 XHC = ((1.0/rcl-1.0)+bha1*G4Log(rcl)+bh 1019 XHC = ((1.0/rcl-1.0)+bha1*G4Log(rcl)+bha2*rl1
1020 + (bha3/2.0)*(rcl*rcl-1.0) 1020 + (bha3/2.0)*(rcl*rcl-1.0)
1021 + (bha4/3.0)*(1.0-rcl*rcl*rcl))/kineti 1021 + (bha4/3.0)*(1.0-rcl*rcl*rcl))/kineticEnergy;
1022 } 1022 }
1023 1023
1024 //Total cross section per molecule for the 1024 //Total cross section per molecule for the active shell, in cm2
1025 G4double XHTOT = XHC + XHDL + XHDT; 1025 G4double XHTOT = XHC + XHDL + XHDT;
1026 1026
1027 //very small cross section, do nothing 1027 //very small cross section, do nothing
1028 if (XHTOT < 1.e-14*barn) 1028 if (XHTOT < 1.e-14*barn)
1029 { 1029 {
1030 fKineticEnergy1=kineticEnergy; 1030 fKineticEnergy1=kineticEnergy;
1031 fCosThetaPrimary=1.0; 1031 fCosThetaPrimary=1.0;
1032 fEnergySecondary=0.0; 1032 fEnergySecondary=0.0;
1033 fCosThetaSecondary=1.0; 1033 fCosThetaSecondary=1.0;
1034 fTargetOscillator = G4int(numberOfOscil 1034 fTargetOscillator = G4int(numberOfOscillators-1);
1035 return; 1035 return;
1036 } 1036 }
1037 1037
1038 //decide which kind of interaction we'll ha 1038 //decide which kind of interaction we'll have
1039 TST = XHTOT*G4UniformRand(); 1039 TST = XHTOT*G4UniformRand();
1040 1040
1041 // Hard close collision 1041 // Hard close collision
1042 G4double TS1 = XHC; 1042 G4double TS1 = XHC;
1043 if (TST < TS1) 1043 if (TST < TS1)
1044 { 1044 {
1045 G4double rl1 = 1.0-rcl; 1045 G4double rl1 = 1.0-rcl;
1046 G4double rk=0.; 1046 G4double rk=0.;
1047 G4bool loopAgain = false; 1047 G4bool loopAgain = false;
1048 do 1048 do
1049 { 1049 {
1050 loopAgain = false; 1050 loopAgain = false;
1051 rk = rcl/(1.0-G4UniformRand()*rl1); 1051 rk = rcl/(1.0-G4UniformRand()*rl1);
1052 G4double phi = 1.0-rk*(bha1-rk*(bha2-rk*( 1052 G4double phi = 1.0-rk*(bha1-rk*(bha2-rk*(bha3-bha4*rk)));
1053 if (G4UniformRand() > phi) 1053 if (G4UniformRand() > phi)
1054 loopAgain = true; 1054 loopAgain = true;
1055 }while(loopAgain); 1055 }while(loopAgain);
1056 //energy and scattering angle (primary 1056 //energy and scattering angle (primary electron)
1057 G4double deltaE = rk*kineticEnergy; 1057 G4double deltaE = rk*kineticEnergy;
1058 fKineticEnergy1 = kineticEnergy - delta 1058 fKineticEnergy1 = kineticEnergy - deltaE;
1059 fCosThetaPrimary = std::sqrt(fKineticEn 1059 fCosThetaPrimary = std::sqrt(fKineticEnergy1*rb/(kineticEnergy*(rb-deltaE)));
1060 //energy and scattering angle of the de 1060 //energy and scattering angle of the delta ray
1061 fEnergySecondary = deltaE - ionEne; //s 1061 fEnergySecondary = deltaE - ionEne; //subtract ionisation energy
1062 fCosThetaSecondary= std::sqrt(deltaE*rb 1062 fCosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*electron_mass_c2)));
1063 if (fVerboseLevel > 3) 1063 if (fVerboseLevel > 3)
1064 G4cout << "SampleFinalStatePositron: sample 1064 G4cout << "SampleFinalStatePositron: sampled close collision " << G4endl;
1065 return; 1065 return;
1066 } 1066 }
1067 1067
1068 //Hard distant longitudinal collisions 1068 //Hard distant longitudinal collisions
1069 TS1 += XHDL; 1069 TS1 += XHDL;
1070 G4double deltaE = resEne; 1070 G4double deltaE = resEne;
1071 fKineticEnergy1 = kineticEnergy - deltaE; 1071 fKineticEnergy1 = kineticEnergy - deltaE;
1072 if (TST < TS1) 1072 if (TST < TS1)
1073 { 1073 {
1074 G4double QS = QM/(1.0+QM*0.5/electron_m 1074 G4double QS = QM/(1.0+QM*0.5/electron_mass_c2);
1075 G4double Q = QS/(std::pow((QS/cutoffEne 1075 G4double Q = QS/(std::pow((QS/cutoffEne)*(1.0+cutoffEne*0.5/electron_mass_c2),G4UniformRand())
1076 - (QS*0.5/electron_mass_c2)); 1076 - (QS*0.5/electron_mass_c2));
1077 G4double QTREV = Q*(Q+2.0*electron_mass 1077 G4double QTREV = Q*(Q+2.0*electron_mass_c2);
1078 G4double cpps = fKineticEnergy1*(fKinet 1078 G4double cpps = fKineticEnergy1*(fKineticEnergy1+2.0*electron_mass_c2);
1079 fCosThetaPrimary = (cpps+cps-QTREV)/(2. 1079 fCosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps));
1080 if (fCosThetaPrimary > 1.) 1080 if (fCosThetaPrimary > 1.)
1081 fCosThetaPrimary = 1.0; 1081 fCosThetaPrimary = 1.0;
1082 //energy and emission angle of the delt 1082 //energy and emission angle of the delta ray
1083 fEnergySecondary = deltaE - ionEne; 1083 fEnergySecondary = deltaE - ionEne;
1084 fCosThetaSecondary = 0.5*(deltaE*(kinet 1084 fCosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
1085 if (fCosThetaSecondary > 1.0) 1085 if (fCosThetaSecondary > 1.0)
1086 fCosThetaSecondary = 1.0; 1086 fCosThetaSecondary = 1.0;
1087 if (fVerboseLevel > 3) 1087 if (fVerboseLevel > 3)
1088 G4cout << "SampleFinalStatePositron: sample 1088 G4cout << "SampleFinalStatePositron: sampled distant longitudinal collision " << G4endl;
1089 return; 1089 return;
1090 } 1090 }
1091 1091
1092 //Hard distant transverse collisions 1092 //Hard distant transverse collisions
1093 fCosThetaPrimary = 1.0; 1093 fCosThetaPrimary = 1.0;
1094 //energy and emission angle of the delta ra 1094 //energy and emission angle of the delta ray
1095 fEnergySecondary = deltaE - ionEne; 1095 fEnergySecondary = deltaE - ionEne;
1096 fCosThetaSecondary = 0.5; 1096 fCosThetaSecondary = 0.5;
1097 1097
1098 if (fVerboseLevel > 3) 1098 if (fVerboseLevel > 3)
1099 G4cout << "SampleFinalStatePositron: samp 1099 G4cout << "SampleFinalStatePositron: sampled distant transverse collision " << G4endl;
1100 1100
1101 return; 1101 return;
1102 } 1102 }
1103 1103
1104 //....oooOO0OOooo........oooOO0OOooo........o 1104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
1105 1105
1106 void G4PenelopeIonisationModel::SetParticle(c 1106 void G4PenelopeIonisationModel::SetParticle(const G4ParticleDefinition* p)
1107 { 1107 {
1108 if(!fParticle) { 1108 if(!fParticle) {
1109 fParticle = p; 1109 fParticle = p;
1110 } 1110 }
1111 } 1111 }
1112 1112