| 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: G4PenelopeComptonModel.cc 99415 2016-09-21 09:05:43Z gcosmo $
26 // 27 //
27 // Author: Luciano Pandola 28 // Author: Luciano Pandola
28 // 29 //
29 // History: 30 // History:
30 // -------- 31 // --------
31 // 15 Feb 2010 L Pandola Implementation 32 // 15 Feb 2010 L Pandola Implementation
32 // 18 Mar 2010 L Pandola Removed GetAtomsPe 33 // 18 Mar 2010 L Pandola Removed GetAtomsPerMolecule(), now demanded
33 // to G4PenelopeOsc 34 // to G4PenelopeOscillatorManager
34 // 01 Feb 2011 L Pandola Suppress fake ener 35 // 01 Feb 2011 L Pandola Suppress fake energy-violation warning when Auger is
35 // active. 36 // active.
36 // Make sure that flu 37 // Make sure that fluorescence/Auger is generated only
37 // if above thresho 38 // if above threshold
38 // 24 May 2011 L Pandola Renamed (make v200 39 // 24 May 2011 L Pandola Renamed (make v2008 as default Penelope)
39 // 10 Jun 2011 L Pandola Migrate atomic dee 40 // 10 Jun 2011 L Pandola Migrate atomic deexcitation interface
40 // 09 Oct 2013 L Pandola Migration to MT 41 // 09 Oct 2013 L Pandola Migration to MT
41 // 25 Jul 2023 D Iuso Fix for possible i <<
42 // 42 //
43 #include "G4PenelopeComptonModel.hh" 43 #include "G4PenelopeComptonModel.hh"
44 #include "G4PhysicalConstants.hh" 44 #include "G4PhysicalConstants.hh"
45 #include "G4SystemOfUnits.hh" 45 #include "G4SystemOfUnits.hh"
46 #include "G4ParticleDefinition.hh" 46 #include "G4ParticleDefinition.hh"
47 #include "G4MaterialCutsCouple.hh" 47 #include "G4MaterialCutsCouple.hh"
48 #include "G4DynamicParticle.hh" 48 #include "G4DynamicParticle.hh"
49 #include "G4VEMDataSet.hh" 49 #include "G4VEMDataSet.hh"
50 #include "G4PhysicsTable.hh" 50 #include "G4PhysicsTable.hh"
51 #include "G4PhysicsLogVector.hh" 51 #include "G4PhysicsLogVector.hh"
52 #include "G4AtomicTransitionManager.hh" 52 #include "G4AtomicTransitionManager.hh"
53 #include "G4AtomicShell.hh" 53 #include "G4AtomicShell.hh"
54 #include "G4Gamma.hh" 54 #include "G4Gamma.hh"
55 #include "G4Electron.hh" 55 #include "G4Electron.hh"
56 #include "G4PenelopeOscillatorManager.hh" 56 #include "G4PenelopeOscillatorManager.hh"
57 #include "G4PenelopeOscillator.hh" 57 #include "G4PenelopeOscillator.hh"
58 #include "G4LossTableManager.hh" 58 #include "G4LossTableManager.hh"
59 #include "G4Exp.hh" 59 #include "G4Exp.hh"
60 60
61 //....oooOO0OOooo........oooOO0OOooo........oo 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
62 62
>> 63
63 G4PenelopeComptonModel::G4PenelopeComptonModel 64 G4PenelopeComptonModel::G4PenelopeComptonModel(const G4ParticleDefinition* part,
64 const G4String& nam) 65 const G4String& nam)
65 :G4VEmModel(nam),fParticleChange(nullptr),fP << 66 :G4VEmModel(nam),fParticleChange(0),fParticle(0),
66 fAtomDeexcitation(nullptr), << 67 isInitialised(false),fAtomDeexcitation(0),
67 fOscManager(nullptr),fIsInitialised(false) << 68 oscManager(0)
68 { 69 {
69 fIntrinsicLowEnergyLimit = 100.0*eV; 70 fIntrinsicLowEnergyLimit = 100.0*eV;
70 fIntrinsicHighEnergyLimit = 100.0*GeV; 71 fIntrinsicHighEnergyLimit = 100.0*GeV;
71 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 72 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
72 // 73 //
73 fOscManager = G4PenelopeOscillatorManager::G << 74 oscManager = G4PenelopeOscillatorManager::GetOscillatorManager();
74 75
75 if (part) 76 if (part)
76 SetParticle(part); 77 SetParticle(part);
77 78
78 fVerboseLevel= 0; << 79 verboseLevel= 0;
79 // Verbosity scale: 80 // Verbosity scale:
80 // 0 = nothing 81 // 0 = nothing
81 // 1 = warning for energy non-conservation 82 // 1 = warning for energy non-conservation
82 // 2 = details of energy budget 83 // 2 = details of energy budget
83 // 3 = calculation of cross sections, file o 84 // 3 = calculation of cross sections, file openings, sampling of atoms
84 // 4 = entering in methods 85 // 4 = entering in methods
85 86
86 //Mark this model as "applicable" for atomic 87 //Mark this model as "applicable" for atomic deexcitation
87 SetDeexcitationFlag(true); 88 SetDeexcitationFlag(true);
88 89
89 fTransitionManager = G4AtomicTransitionManag 90 fTransitionManager = G4AtomicTransitionManager::Instance();
90 } 91 }
91 92
92 //....oooOO0OOooo........oooOO0OOooo........oo 93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
93 94
94 G4PenelopeComptonModel::~G4PenelopeComptonMode 95 G4PenelopeComptonModel::~G4PenelopeComptonModel()
95 {;} 96 {;}
96 97
97 //....oooOO0OOooo........oooOO0OOooo........oo 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
98 99
99 void G4PenelopeComptonModel::Initialise(const 100 void G4PenelopeComptonModel::Initialise(const G4ParticleDefinition* part,
100 const G4DataVector&) 101 const G4DataVector&)
101 { 102 {
102 if (fVerboseLevel > 3) << 103 if (verboseLevel > 3)
103 G4cout << "Calling G4PenelopeComptonModel: 104 G4cout << "Calling G4PenelopeComptonModel::Initialise()" << G4endl;
104 105
105 fAtomDeexcitation = G4LossTableManager::Inst 106 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
106 //Issue warning if the AtomicDeexcitation ha 107 //Issue warning if the AtomicDeexcitation has not been declared
107 if (!fAtomDeexcitation) 108 if (!fAtomDeexcitation)
108 { 109 {
109 G4cout << G4endl; 110 G4cout << G4endl;
110 G4cout << "WARNING from G4PenelopeCompto 111 G4cout << "WARNING from G4PenelopeComptonModel " << G4endl;
111 G4cout << "Atomic de-excitation module i 112 G4cout << "Atomic de-excitation module is not instantiated, so there will not be ";
112 G4cout << "any fluorescence/Auger emissi 113 G4cout << "any fluorescence/Auger emission." << G4endl;
113 G4cout << "Please make sure this is inte 114 G4cout << "Please make sure this is intended" << G4endl;
114 } 115 }
115 116
116 SetParticle(part); 117 SetParticle(part);
117 118
118 if (IsMaster() && part == fParticle) 119 if (IsMaster() && part == fParticle)
119 { 120 {
120 121
121 if (fVerboseLevel > 0) << 122 if (verboseLevel > 0)
122 { 123 {
123 G4cout << "Penelope Compton model v2008 is 124 G4cout << "Penelope Compton model v2008 is initialized " << G4endl
124 << "Energy range: " 125 << "Energy range: "
125 << LowEnergyLimit() / keV << " keV - " 126 << LowEnergyLimit() / keV << " keV - "
126 << HighEnergyLimit() / GeV << " GeV"; 127 << HighEnergyLimit() / GeV << " GeV";
127 } 128 }
128 //Issue a warning, if the model is going 129 //Issue a warning, if the model is going to be used down to a
129 //energy which is outside the validity o 130 //energy which is outside the validity of the model itself
130 if (LowEnergyLimit() < fIntrinsicLowEner 131 if (LowEnergyLimit() < fIntrinsicLowEnergyLimit)
131 { 132 {
132 G4ExceptionDescription ed; 133 G4ExceptionDescription ed;
133 ed << "Using the Penelope Compton model ou 134 ed << "Using the Penelope Compton model outside its intrinsic validity range. "
134 << G4endl; 135 << G4endl;
135 ed << "-> LowEnergyLimit() in process = " 136 ed << "-> LowEnergyLimit() in process = " << LowEnergyLimit()/keV << "keV " << G4endl;
136 ed << "-> Instrinsic low-energy limit = " 137 ed << "-> Instrinsic low-energy limit = " << fIntrinsicLowEnergyLimit/keV << "keV "
137 << G4endl; 138 << G4endl;
138 ed << "Result of the simulation have to be 139 ed << "Result of the simulation have to be taken with care" << G4endl;
139 G4Exception("G4PenelopeComptonModel::Initi 140 G4Exception("G4PenelopeComptonModel::Initialise()",
140 "em2100",JustWarning,ed); 141 "em2100",JustWarning,ed);
141 } 142 }
142 } 143 }
143 144
144 if(fIsInitialised) return; << 145 if(isInitialised) return;
145 fParticleChange = GetParticleChangeForGamma( 146 fParticleChange = GetParticleChangeForGamma();
146 fIsInitialised = true; << 147 isInitialised = true;
147 148
148 } 149 }
149 150
150 //....oooOO0OOooo........oooOO0OOooo........oo 151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
151 152
152 void G4PenelopeComptonModel::InitialiseLocal(c 153 void G4PenelopeComptonModel::InitialiseLocal(const G4ParticleDefinition* part,
153 154 G4VEmModel *masterModel)
154 { 155 {
155 if (fVerboseLevel > 3) << 156 if (verboseLevel > 3)
156 G4cout << "Calling G4PenelopeComptonModel 157 G4cout << "Calling G4PenelopeComptonModel::InitialiseLocal()" << G4endl;
>> 158
157 // 159 //
158 //Check that particle matches: one might hav 160 //Check that particle matches: one might have multiple master models (e.g.
159 //for e+ and e-). 161 //for e+ and e-).
160 // 162 //
161 if (part == fParticle) 163 if (part == fParticle)
162 { 164 {
163 //Get the const table pointers from the 165 //Get the const table pointers from the master to the workers
164 const G4PenelopeComptonModel* theModel = 166 const G4PenelopeComptonModel* theModel =
165 static_cast<G4PenelopeComptonModel*> ( 167 static_cast<G4PenelopeComptonModel*> (masterModel);
166 168
167 //Same verbosity for all workers, as the 169 //Same verbosity for all workers, as the master
168 fVerboseLevel = theModel->fVerboseLevel; << 170 verboseLevel = theModel->verboseLevel;
169 } 171 }
>> 172
170 return; 173 return;
171 } 174 }
172 175
173 176
174 //....oooOO0OOooo........oooOO0OOooo........oo 177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
175 178
176 G4double G4PenelopeComptonModel::CrossSectionP 179 G4double G4PenelopeComptonModel::CrossSectionPerVolume(const G4Material* material,
177 180 const G4ParticleDefinition* p,
178 181 G4double energy,
179 182 G4double,
180 183 G4double)
181 { 184 {
182 // Penelope model v2008 to calculate the Com 185 // Penelope model v2008 to calculate the Compton scattering cross section:
183 // D. Brusa et al., Nucl. Instrum. Meth. A 3 186 // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
184 // 187 //
185 // The cross section for Compton scattering 188 // The cross section for Compton scattering is calculated according to the Klein-Nishina
186 // formula for energy > 5 MeV. 189 // formula for energy > 5 MeV.
187 // For E < 5 MeV it is used a parametrizatio 190 // For E < 5 MeV it is used a parametrization for the differential cross-section dSigma/dOmega,
188 // which is integrated numerically in cos(th 191 // which is integrated numerically in cos(theta), G4PenelopeComptonModel::DifferentialCrossSection().
189 // The parametrization includes the J(p) 192 // The parametrization includes the J(p)
190 // distribution profiles for the atomic shel 193 // distribution profiles for the atomic shells, that are tabulated from Hartree-Fock calculations
191 // from F. Biggs et al., At. Data Nucl. Data 194 // from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201
192 // 195 //
193 if (fVerboseLevel > 3) << 196 if (verboseLevel > 3)
194 G4cout << "Calling CrossSectionPerVolume() 197 G4cout << "Calling CrossSectionPerVolume() of G4PenelopeComptonModel" << G4endl;
195 SetupForMaterial(p, material, energy); 198 SetupForMaterial(p, material, energy);
196 199
>> 200
197 G4double cs = 0; 201 G4double cs = 0;
198 //Force null cross-section if below the low- 202 //Force null cross-section if below the low-energy edge of the table
199 if (energy < LowEnergyLimit()) 203 if (energy < LowEnergyLimit())
200 return cs; 204 return cs;
201 205
202 //Retrieve the oscillator table for this mat 206 //Retrieve the oscillator table for this material
203 G4PenelopeOscillatorTable* theTable = fOscMa << 207 G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableCompton(material);
204 208
205 if (energy < 5*MeV) //explicit calculation f 209 if (energy < 5*MeV) //explicit calculation for E < 5 MeV
206 { 210 {
207 size_t numberOfOscillators = theTable->s 211 size_t numberOfOscillators = theTable->size();
208 for (size_t i=0;i<numberOfOscillators;i+ 212 for (size_t i=0;i<numberOfOscillators;i++)
209 { 213 {
210 G4PenelopeOscillator* theOsc = (*theTable) 214 G4PenelopeOscillator* theOsc = (*theTable)[i];
211 //sum contributions coming from each oscil 215 //sum contributions coming from each oscillator
212 cs += OscillatorTotalCrossSection(energy,t 216 cs += OscillatorTotalCrossSection(energy,theOsc);
213 } 217 }
214 } 218 }
215 else //use Klein-Nishina for E>5 MeV 219 else //use Klein-Nishina for E>5 MeV
216 cs = KleinNishinaCrossSection(energy,mater 220 cs = KleinNishinaCrossSection(energy,material);
217 221
218 //cross sections are in units of pi*classic_ 222 //cross sections are in units of pi*classic_electr_radius^2
219 cs *= pi*classic_electr_radius*classic_elect 223 cs *= pi*classic_electr_radius*classic_electr_radius;
220 224
221 //Now, cs is the cross section *per molecule 225 //Now, cs is the cross section *per molecule*, let's calculate the
222 //cross section per volume 226 //cross section per volume
>> 227
223 G4double atomDensity = material->GetTotNbOfA 228 G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
224 G4double atPerMol = fOscManager->GetAtomsPe << 229 G4double atPerMol = oscManager->GetAtomsPerMolecule(material);
225 230
226 if (fVerboseLevel > 3) << 231 if (verboseLevel > 3)
227 G4cout << "Material " << material->GetName 232 G4cout << "Material " << material->GetName() << " has " << atPerMol <<
228 "atoms per molecule" << G4endl; 233 "atoms per molecule" << G4endl;
229 234
230 G4double moleculeDensity = 0.; 235 G4double moleculeDensity = 0.;
231 236
232 if (atPerMol) 237 if (atPerMol)
233 moleculeDensity = atomDensity/atPerMol; 238 moleculeDensity = atomDensity/atPerMol;
234 239
235 G4double csvolume = cs*moleculeDensity; 240 G4double csvolume = cs*moleculeDensity;
236 241
237 if (fVerboseLevel > 2) << 242 if (verboseLevel > 2)
238 G4cout << "Compton mean free path at " << 243 G4cout << "Compton mean free path at " << energy/keV << " keV for material " <<
239 material->GetName() << " = " << (1 244 material->GetName() << " = " << (1./csvolume)/mm << " mm" << G4endl;
240 return csvolume; 245 return csvolume;
241 } 246 }
242 247
243 //....oooOO0OOooo........oooOO0OOooo........oo 248 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
244 249
245 //This is a dummy method. Never inkoved by the 250 //This is a dummy method. Never inkoved by the tracking, it just issues
246 //a warning if one tries to get Cross Sections 251 //a warning if one tries to get Cross Sections per Atom via the
247 //G4EmCalculator. 252 //G4EmCalculator.
248 G4double G4PenelopeComptonModel::ComputeCrossS 253 G4double G4PenelopeComptonModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
249 254 G4double,
250 255 G4double,
251 256 G4double,
252 257 G4double,
253 258 G4double)
254 { 259 {
255 G4cout << "*** G4PenelopeComptonModel -- WAR 260 G4cout << "*** G4PenelopeComptonModel -- WARNING ***" << G4endl;
256 G4cout << "Penelope Compton model v2008 does 261 G4cout << "Penelope Compton model v2008 does not calculate cross section _per atom_ " << G4endl;
257 G4cout << "so the result is always zero. For 262 G4cout << "so the result is always zero. For physics values, please invoke " << G4endl;
258 G4cout << "GetCrossSectionPerVolume() or Get 263 G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl;
259 return 0; 264 return 0;
260 } 265 }
261 266
262 //....oooOO0OOooo........oooOO0OOooo........oo 267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
263 268
264 void G4PenelopeComptonModel::SampleSecondaries 269 void G4PenelopeComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
265 const G4MaterialCutsCouple* c 270 const G4MaterialCutsCouple* couple,
266 const G4DynamicParticle* aDyna 271 const G4DynamicParticle* aDynamicGamma,
267 G4double, 272 G4double,
268 G4double) 273 G4double)
269 { 274 {
>> 275
270 // Penelope model v2008 to sample the Compto 276 // Penelope model v2008 to sample the Compton scattering final state.
271 // D. Brusa et al., Nucl. Instrum. Meth. A 3 277 // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
272 // The model determines also the original sh 278 // The model determines also the original shell from which the electron is expelled,
273 // in order to produce fluorescence de-excit 279 // in order to produce fluorescence de-excitation (from G4DeexcitationManager)
274 // 280 //
275 // The final state for Compton scattering is 281 // The final state for Compton scattering is calculated according to the Klein-Nishina
276 // formula for energy > 5 MeV. In this case, 282 // formula for energy > 5 MeV. In this case, the Doppler broadening is negligible and
277 // one can assume that the target electron i 283 // one can assume that the target electron is at rest.
278 // For E < 5 MeV it is used the parametrizat 284 // For E < 5 MeV it is used the parametrization for the differential cross-section dSigma/dOmega,
279 // to sample the scattering angle and the en 285 // to sample the scattering angle and the energy of the emerging electron, which is
280 // G4PenelopeComptonModel::DifferentialCross 286 // G4PenelopeComptonModel::DifferentialCrossSection(). The rejection method is
281 // used to sample cos(theta). The efficiency 287 // used to sample cos(theta). The efficiency increases monotonically with photon energy and is
282 // nearly independent on the Z; typical valu 288 // nearly independent on the Z; typical values are 35%, 80% and 95% for 1 keV, 1 MeV and 10 MeV,
283 // respectively. 289 // respectively.
284 // The parametrization includes the J(p) dis 290 // The parametrization includes the J(p) distribution profiles for the atomic shells, that are
285 // tabulated 291 // tabulated
286 // from Hartree-Fock calculations from F. Bi 292 // from Hartree-Fock calculations from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201.
287 // Doppler broadening is included. 293 // Doppler broadening is included.
288 // 294 //
289 295
290 if (fVerboseLevel > 3) << 296 if (verboseLevel > 3)
291 G4cout << "Calling SampleSecondaries() of 297 G4cout << "Calling SampleSecondaries() of G4PenelopeComptonModel" << G4endl;
292 298
293 G4double photonEnergy0 = aDynamicGamma->GetK 299 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
294 300
295 // do nothing below the threshold 301 // do nothing below the threshold
296 // should never get here because the XS is z 302 // should never get here because the XS is zero below the limit
297 if(photonEnergy0 < LowEnergyLimit()) 303 if(photonEnergy0 < LowEnergyLimit())
298 return; 304 return;
299 305
300 G4ParticleMomentum photonDirection0 = aDynam 306 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
301 const G4Material* material = couple->GetMate 307 const G4Material* material = couple->GetMaterial();
302 308
303 G4PenelopeOscillatorTable* theTable = fOscMa << 309 G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableCompton(material);
304 310
305 const G4int nmax = 64; 311 const G4int nmax = 64;
306 G4double rn[nmax]={0.0}; 312 G4double rn[nmax]={0.0};
307 G4double pac[nmax]={0.0}; 313 G4double pac[nmax]={0.0};
308 314
309 G4double S=0.0; 315 G4double S=0.0;
310 G4double epsilon = 0.0; 316 G4double epsilon = 0.0;
311 G4double cosTheta = 1.0; 317 G4double cosTheta = 1.0;
312 G4double hartreeFunc = 0.0; 318 G4double hartreeFunc = 0.0;
313 G4double oscStren = 0.0; 319 G4double oscStren = 0.0;
314 size_t numberOfOscillators = theTable->size( 320 size_t numberOfOscillators = theTable->size();
315 size_t targetOscillator = 0; 321 size_t targetOscillator = 0;
316 G4double ionEnergy = 0.0*eV; 322 G4double ionEnergy = 0.0*eV;
317 323
318 G4double ek = photonEnergy0/electron_mass_c2 324 G4double ek = photonEnergy0/electron_mass_c2;
319 G4double ek2 = 2.*ek+1.0; 325 G4double ek2 = 2.*ek+1.0;
320 G4double eks = ek*ek; 326 G4double eks = ek*ek;
321 G4double ek1 = eks-ek2-1.0; 327 G4double ek1 = eks-ek2-1.0;
322 328
323 G4double taumin = 1.0/ek2; 329 G4double taumin = 1.0/ek2;
324 //This is meant to fix a possible (rare) flo << 330 G4double a1 = std::log(ek2);
325 //causing an infinite loop. The maximum of t <<
326 //be represented (i.e. ~ 1. - 1e-16). Fix by <<
327 static G4double taumax = std::nexttoward(1.0 <<
328 if (fVerboseLevel > 3) <<
329 G4cout << "G4PenelopeComptonModel: maximum <<
330 //To here. <<
331 G4double a1 = G4Log(ek2); <<
332 G4double a2 = a1+2.0*ek*(1.0+ek)/(ek2*ek2); 331 G4double a2 = a1+2.0*ek*(1.0+ek)/(ek2*ek2);
333 332
334 G4double TST = 0; 333 G4double TST = 0;
335 G4double tau = 0.; 334 G4double tau = 0.;
336 335
337 //If the incoming photon is above 5 MeV, the 336 //If the incoming photon is above 5 MeV, the quicker approach based on the
338 //pure Klein-Nishina formula is used 337 //pure Klein-Nishina formula is used
339 if (photonEnergy0 > 5*MeV) 338 if (photonEnergy0 > 5*MeV)
340 { 339 {
341 do{ 340 do{
342 do{ 341 do{
343 if ((a2*G4UniformRand()) < a1) 342 if ((a2*G4UniformRand()) < a1)
344 tau = std::pow(taumin,G4UniformRand()); 343 tau = std::pow(taumin,G4UniformRand());
345 else 344 else
346 tau = std::sqrt(1.0+G4UniformRand()*(tau 345 tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));
347 //rejection function 346 //rejection function
348 TST = (1.0+tau*(ek1+tau*(ek2+tau*eks)))/(e 347 TST = (1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau));
349 }while (G4UniformRand()> TST); 348 }while (G4UniformRand()> TST);
350 if (tau > taumax) tau = taumax; //prevent FP <<
351 epsilon=tau; 349 epsilon=tau;
352 cosTheta = 1.0 - (1.0-tau)/(ek*tau); 350 cosTheta = 1.0 - (1.0-tau)/(ek*tau);
353 351
354 //Target shell electrons 352 //Target shell electrons
355 TST = fOscManager->GetTotalZ(material)*G4Uni << 353 TST = oscManager->GetTotalZ(material)*G4UniformRand();
356 targetOscillator = numberOfOscillators-1; // 354 targetOscillator = numberOfOscillators-1; //last level
357 S=0.0; 355 S=0.0;
358 G4bool levelFound = false; 356 G4bool levelFound = false;
359 for (size_t j=0;j<numberOfOscillators && !le 357 for (size_t j=0;j<numberOfOscillators && !levelFound; j++)
360 { 358 {
361 S += (*theTable)[j]->GetOscillatorStreng 359 S += (*theTable)[j]->GetOscillatorStrength();
362 if (S > TST) 360 if (S > TST)
363 { 361 {
364 targetOscillator = j; 362 targetOscillator = j;
365 levelFound = true; 363 levelFound = true;
366 } 364 }
367 } 365 }
368 //check whether the level is valid 366 //check whether the level is valid
369 ionEnergy = (*theTable)[targetOscillator]->G 367 ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy();
370 }while((epsilon*photonEnergy0-photonEner 368 }while((epsilon*photonEnergy0-photonEnergy0+ionEnergy) >0);
371 } 369 }
372 else //photonEnergy0 < 5 MeV 370 else //photonEnergy0 < 5 MeV
373 { 371 {
374 //Incoherent scattering function for the 372 //Incoherent scattering function for theta=PI
375 G4double s0=0.0; 373 G4double s0=0.0;
376 G4double pzomc=0.0; 374 G4double pzomc=0.0;
377 G4double rni=0.0; 375 G4double rni=0.0;
378 G4double aux=0.0; 376 G4double aux=0.0;
379 for (size_t i=0;i<numberOfOscillators;i+ 377 for (size_t i=0;i<numberOfOscillators;i++)
380 { 378 {
381 ionEnergy = (*theTable)[i]->GetIonisationE 379 ionEnergy = (*theTable)[i]->GetIonisationEnergy();
382 if (photonEnergy0 > ionEnergy) 380 if (photonEnergy0 > ionEnergy)
383 { 381 {
384 G4double aux2 = photonEnergy0*(photonE 382 G4double aux2 = photonEnergy0*(photonEnergy0-ionEnergy)*2.0;
385 hartreeFunc = (*theTable)[i]->GetHartr 383 hartreeFunc = (*theTable)[i]->GetHartreeFactor();
386 oscStren = (*theTable)[i]->GetOscillat 384 oscStren = (*theTable)[i]->GetOscillatorStrength();
387 pzomc = hartreeFunc*(aux2-electron_mas 385 pzomc = hartreeFunc*(aux2-electron_mass_c2*ionEnergy)/
388 (electron_mass_c2*std::sqrt(2.0*aux2+ionEn 386 (electron_mass_c2*std::sqrt(2.0*aux2+ionEnergy*ionEnergy));
389 if (pzomc > 0) 387 if (pzomc > 0)
390 rni = 1.0-0.5*G4Exp(0.5-(std::sqrt(0.5)+st 388 rni = 1.0-0.5*G4Exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
391 (std::sqrt(0.5)+std::sqrt(2.0)* 389 (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
392 else 390 else
393 rni = 0.5*G4Exp(0.5-(std::sqrt(0.5)-std::s 391 rni = 0.5*G4Exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*
394 (std::sqrt(0.5)-std::sqrt(2.0)*pzom 392 (std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
395 s0 += oscStren*rni; 393 s0 += oscStren*rni;
396 } 394 }
397 } 395 }
398 //Sampling tau 396 //Sampling tau
399 G4double cdt1 = 0.; 397 G4double cdt1 = 0.;
400 do 398 do
401 { 399 {
402 if ((G4UniformRand()*a2) < a1) 400 if ((G4UniformRand()*a2) < a1)
403 tau = std::pow(taumin,G4UniformRand()); 401 tau = std::pow(taumin,G4UniformRand());
404 else 402 else
405 tau = std::sqrt(1.0+G4UniformRand()*(tau 403 tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));
406 if (tau > taumax) tau = taumax; //prevent <<
407 cdt1 = (1.0-tau)/(ek*tau); 404 cdt1 = (1.0-tau)/(ek*tau);
408 //Incoherent scattering function 405 //Incoherent scattering function
409 S = 0.; 406 S = 0.;
410 for (size_t i=0;i<numberOfOscillators;i++) 407 for (size_t i=0;i<numberOfOscillators;i++)
411 { 408 {
412 ionEnergy = (*theTable)[i]->GetIonisat 409 ionEnergy = (*theTable)[i]->GetIonisationEnergy();
413 if (photonEnergy0 > ionEnergy) //sum o 410 if (photonEnergy0 > ionEnergy) //sum only on excitable levels
414 { 411 {
415 aux = photonEnergy0*(photonEnergy0-ionEn 412 aux = photonEnergy0*(photonEnergy0-ionEnergy)*cdt1;
416 hartreeFunc = (*theTable)[i]->GetHartree 413 hartreeFunc = (*theTable)[i]->GetHartreeFactor();
417 oscStren = (*theTable)[i]->GetOscillator 414 oscStren = (*theTable)[i]->GetOscillatorStrength();
418 pzomc = hartreeFunc*(aux-electron_mass_c 415 pzomc = hartreeFunc*(aux-electron_mass_c2*ionEnergy)/
419 (electron_mass_c2*std::sqrt(2.0*aux+io 416 (electron_mass_c2*std::sqrt(2.0*aux+ionEnergy*ionEnergy));
420 if (pzomc > 0) 417 if (pzomc > 0)
421 rn[i] = 1.0-0.5*G4Exp(0.5-(std::sqrt(0 418 rn[i] = 1.0-0.5*G4Exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*
422 (std::sqrt(0.5)+std::sqrt(2.0)* 419 (std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
423 else 420 else
424 rn[i] = 0.5*G4Exp(0.5-(std::sqrt(0.5)- 421 rn[i] = 0.5*G4Exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*
425 (std::sqrt(0.5)-std::sqrt(2.0)*pzom 422 (std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
426 S += oscStren*rn[i]; 423 S += oscStren*rn[i];
427 pac[i] = S; 424 pac[i] = S;
428 } 425 }
429 else 426 else
430 pac[i] = S-1e-6; 427 pac[i] = S-1e-6;
431 } 428 }
432 //Rejection function 429 //Rejection function
433 TST = S*(1.0+tau*(ek1+tau*(ek2+tau*eks)))/ 430 TST = S*(1.0+tau*(ek1+tau*(ek2+tau*eks)))/(eks*tau*(1.0+tau*tau));
434 }while ((G4UniformRand()*s0) > TST); 431 }while ((G4UniformRand()*s0) > TST);
435 432
436 cosTheta = 1.0 - cdt1; 433 cosTheta = 1.0 - cdt1;
437 G4double fpzmax=0.0,fpz=0.0; 434 G4double fpzmax=0.0,fpz=0.0;
438 G4double A=0.0; 435 G4double A=0.0;
439 //Target electron shell 436 //Target electron shell
440 do 437 do
441 { 438 {
442 do 439 do
443 { 440 {
444 TST = S*G4UniformRand(); 441 TST = S*G4UniformRand();
445 targetOscillator = numberOfOscillators 442 targetOscillator = numberOfOscillators-1; //last level
446 G4bool levelFound = false; 443 G4bool levelFound = false;
447 for (size_t i=0;i<numberOfOscillators 444 for (size_t i=0;i<numberOfOscillators && !levelFound;i++)
448 { 445 {
449 if (pac[i]>TST) 446 if (pac[i]>TST)
450 { 447 {
451 targetOscillator = i; 448 targetOscillator = i;
452 levelFound = true; 449 levelFound = true;
453 } 450 }
454 } 451 }
455 A = G4UniformRand()*rn[targetOscillato 452 A = G4UniformRand()*rn[targetOscillator];
456 hartreeFunc = (*theTable)[targetOscill 453 hartreeFunc = (*theTable)[targetOscillator]->GetHartreeFactor();
457 oscStren = (*theTable)[targetOscillato 454 oscStren = (*theTable)[targetOscillator]->GetOscillatorStrength();
458 if (A < 0.5) 455 if (A < 0.5)
459 pzomc = (std::sqrt(0.5)-std::sqrt(0.5-G4Lo << 456 pzomc = (std::sqrt(0.5)-std::sqrt(0.5-std::log(2.0*A)))/
460 (std::sqrt(2.0)*hartreeFunc); 457 (std::sqrt(2.0)*hartreeFunc);
461 else 458 else
462 pzomc = (std::sqrt(0.5-G4Log(2.0-2.0*A))-s << 459 pzomc = (std::sqrt(0.5-std::log(2.0-2.0*A))-std::sqrt(0.5))/
463 (std::sqrt(2.0)*hartreeFunc); 460 (std::sqrt(2.0)*hartreeFunc);
464 } while (pzomc < -1); 461 } while (pzomc < -1);
465 462
466 // F(EP) rejection 463 // F(EP) rejection
467 G4double XQC = 1.0+tau*(tau-2.0*cosTheta); 464 G4double XQC = 1.0+tau*(tau-2.0*cosTheta);
468 G4double AF = std::sqrt(XQC)*(1.0+tau*(tau 465 G4double AF = std::sqrt(XQC)*(1.0+tau*(tau-cosTheta)/XQC);
469 if (AF > 0) 466 if (AF > 0)
470 fpzmax = 1.0+AF*0.2; 467 fpzmax = 1.0+AF*0.2;
471 else 468 else
472 fpzmax = 1.0-AF*0.2; 469 fpzmax = 1.0-AF*0.2;
473 fpz = 1.0+AF*std::max(std::min(pzomc,0.2), 470 fpz = 1.0+AF*std::max(std::min(pzomc,0.2),-0.2);
474 }while ((fpzmax*G4UniformRand())>fpz); 471 }while ((fpzmax*G4UniformRand())>fpz);
475 472
476 //Energy of the scattered photon 473 //Energy of the scattered photon
477 G4double T = pzomc*pzomc; 474 G4double T = pzomc*pzomc;
478 G4double b1 = 1.0-T*tau*tau; 475 G4double b1 = 1.0-T*tau*tau;
479 G4double b2 = 1.0-T*tau*cosTheta; 476 G4double b2 = 1.0-T*tau*cosTheta;
480 if (pzomc > 0.0) 477 if (pzomc > 0.0)
481 epsilon = (tau/b1)*(b2+std::sqrt(std::abs(b2 478 epsilon = (tau/b1)*(b2+std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
482 else 479 else
483 epsilon = (tau/b1)*(b2-std::sqrt(std::abs(b2 480 epsilon = (tau/b1)*(b2-std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
484 } //energy < 5 MeV 481 } //energy < 5 MeV
485 482
486 //Ok, the kinematics has been calculated. 483 //Ok, the kinematics has been calculated.
487 G4double sinTheta = std::sqrt(1-cosTheta*cos 484 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
488 G4double phi = twopi * G4UniformRand() ; 485 G4double phi = twopi * G4UniformRand() ;
489 G4double dirx = sinTheta * std::cos(phi); 486 G4double dirx = sinTheta * std::cos(phi);
490 G4double diry = sinTheta * std::sin(phi); 487 G4double diry = sinTheta * std::sin(phi);
491 G4double dirz = cosTheta ; 488 G4double dirz = cosTheta ;
492 489
493 // Update G4VParticleChange for the scattere 490 // Update G4VParticleChange for the scattered photon
494 G4ThreeVector photonDirection1(dirx,diry,dir 491 G4ThreeVector photonDirection1(dirx,diry,dirz);
495 photonDirection1.rotateUz(photonDirection0); 492 photonDirection1.rotateUz(photonDirection0);
496 fParticleChange->ProposeMomentumDirection(ph 493 fParticleChange->ProposeMomentumDirection(photonDirection1) ;
497 494
498 G4double photonEnergy1 = epsilon * photonEne 495 G4double photonEnergy1 = epsilon * photonEnergy0;
499 496
500 if (photonEnergy1 > 0.) 497 if (photonEnergy1 > 0.)
501 fParticleChange->SetProposedKineticEnergy( 498 fParticleChange->SetProposedKineticEnergy(photonEnergy1) ;
502 else 499 else
503 { 500 {
504 fParticleChange->SetProposedKineticEnergy( 501 fParticleChange->SetProposedKineticEnergy(0.) ;
505 fParticleChange->ProposeTrackStatus(fStopA 502 fParticleChange->ProposeTrackStatus(fStopAndKill);
506 } 503 }
507 504
508 //Create scattered electron 505 //Create scattered electron
509 G4double diffEnergy = photonEnergy0*(1-epsil 506 G4double diffEnergy = photonEnergy0*(1-epsilon);
510 ionEnergy = (*theTable)[targetOscillator]->G 507 ionEnergy = (*theTable)[targetOscillator]->GetIonisationEnergy();
511 508
512 G4double Q2 = 509 G4double Q2 =
513 photonEnergy0*photonEnergy0+photonEnergy1* 510 photonEnergy0*photonEnergy0+photonEnergy1*(photonEnergy1-2.0*photonEnergy0*cosTheta);
514 G4double cosThetaE = 0.; //scattering angle 511 G4double cosThetaE = 0.; //scattering angle for the electron
515 512
516 if (Q2 > 1.0e-12) 513 if (Q2 > 1.0e-12)
517 cosThetaE = (photonEnergy0-photonEnergy1*c 514 cosThetaE = (photonEnergy0-photonEnergy1*cosTheta)/std::sqrt(Q2);
518 else 515 else
519 cosThetaE = 1.0; 516 cosThetaE = 1.0;
520 G4double sinThetaE = std::sqrt(1-cosThetaE*c 517 G4double sinThetaE = std::sqrt(1-cosThetaE*cosThetaE);
521 518
522 //Now, try to handle fluorescence 519 //Now, try to handle fluorescence
523 //Notice: merged levels are indicated with Z 520 //Notice: merged levels are indicated with Z=0 and flag=30
524 G4int shFlag = (*theTable)[targetOscillator] 521 G4int shFlag = (*theTable)[targetOscillator]->GetShellFlag();
525 G4int Z = (G4int) (*theTable)[targetOscillat 522 G4int Z = (G4int) (*theTable)[targetOscillator]->GetParentZ();
526 523
527 //initialize here, then check photons create 524 //initialize here, then check photons created by Atomic-Deexcitation, and the final state e-
528 G4double bindingEnergy = 0.*eV; 525 G4double bindingEnergy = 0.*eV;
529 const G4AtomicShell* shell = 0; 526 const G4AtomicShell* shell = 0;
530 527
531 //Real level 528 //Real level
532 if (Z > 0 && shFlag<30) 529 if (Z > 0 && shFlag<30)
533 { 530 {
534 shell = fTransitionManager->Shell(Z,shFl 531 shell = fTransitionManager->Shell(Z,shFlag-1);
535 bindingEnergy = shell->BindingEnergy(); 532 bindingEnergy = shell->BindingEnergy();
536 } 533 }
537 534
538 G4double ionEnergyInPenelopeDatabase = ionEn 535 G4double ionEnergyInPenelopeDatabase = ionEnergy;
539 //protection against energy non-conservation 536 //protection against energy non-conservation
540 ionEnergy = std::max(bindingEnergy,ionEnergy 537 ionEnergy = std::max(bindingEnergy,ionEnergyInPenelopeDatabase);
541 538
542 //subtract the excitation energy. If not emi 539 //subtract the excitation energy. If not emitted by fluorescence
543 //the ionization energy is deposited as loca 540 //the ionization energy is deposited as local energy deposition
544 G4double eKineticEnergy = diffEnergy - ionEn 541 G4double eKineticEnergy = diffEnergy - ionEnergy;
545 G4double localEnergyDeposit = ionEnergy; 542 G4double localEnergyDeposit = ionEnergy;
546 G4double energyInFluorescence = 0.; //testin 543 G4double energyInFluorescence = 0.; //testing purposes only
547 G4double energyInAuger = 0; //testing purpos 544 G4double energyInAuger = 0; //testing purposes
548 545
549 if (eKineticEnergy < 0) 546 if (eKineticEnergy < 0)
550 { 547 {
551 //It means that there was some problem/m 548 //It means that there was some problem/mismatch between the two databases.
552 //Try to make it work 549 //Try to make it work
553 //In this case available Energy (diffEne 550 //In this case available Energy (diffEnergy) < ionEnergy
554 //Full residual energy is deposited loca 551 //Full residual energy is deposited locally
555 localEnergyDeposit = diffEnergy; 552 localEnergyDeposit = diffEnergy;
556 eKineticEnergy = 0.0; 553 eKineticEnergy = 0.0;
557 } 554 }
558 555
559 //the local energy deposit is what remains: 556 //the local energy deposit is what remains: part of this may be spent for fluorescence.
560 //Notice: shell might be nullptr (invalid!) 557 //Notice: shell might be nullptr (invalid!) if shFlag=30. Must be protected
561 //Now, take care of fluorescence, if require 558 //Now, take care of fluorescence, if required
562 if (fAtomDeexcitation && shell) 559 if (fAtomDeexcitation && shell)
563 { 560 {
564 G4int index = couple->GetIndex(); 561 G4int index = couple->GetIndex();
565 if (fAtomDeexcitation->CheckDeexcitation 562 if (fAtomDeexcitation->CheckDeexcitationActiveRegion(index))
566 { 563 {
567 size_t nBefore = fvect->size(); 564 size_t nBefore = fvect->size();
568 fAtomDeexcitation->GenerateParticles(fvect 565 fAtomDeexcitation->GenerateParticles(fvect,shell,Z,index);
569 size_t nAfter = fvect->size(); 566 size_t nAfter = fvect->size();
570 567
571 if (nAfter > nBefore) //actual production 568 if (nAfter > nBefore) //actual production of fluorescence
572 { 569 {
573 for (size_t j=nBefore;j<nAfter;j++) // 570 for (size_t j=nBefore;j<nAfter;j++) //loop on products
574 { 571 {
575 G4double itsEnergy = ((*fvect)[j])->GetK 572 G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy();
576 if (itsEnergy < localEnergyDeposit) // v << 573 localEnergyDeposit -= itsEnergy;
577 { << 574 if (((*fvect)[j])->GetParticleDefinition() == G4Gamma::Definition())
578 localEnergyDeposit -= itsEnergy; << 575 energyInFluorescence += itsEnergy;
579 if (((*fvect)[j])->GetParticleDefini << 576 else if (((*fvect)[j])->GetParticleDefinition() == G4Electron::Definition())
580 energyInFluorescence += itsEnergy; << 577 energyInAuger += itsEnergy;
581 else if (((*fvect)[j])->GetParticleD << 578
582 G4Electron::Definition()) <<
583 energyInAuger += itsEnergy; <<
584 } <<
585 else //invalid secondary: takes more tha <<
586 { <<
587 delete (*fvect)[j]; <<
588 (*fvect)[j] = nullptr; <<
589 } <<
590 } 579 }
591 } 580 }
592 << 581
593 } 582 }
594 } 583 }
595 584
596 //Always produce explicitly the electron << 585
>> 586 /*
>> 587 if(DeexcitationFlag() && Z > 5 && shellId>0) {
>> 588
>> 589 const G4ProductionCutsTable* theCoupleTable=
>> 590 G4ProductionCutsTable::GetProductionCutsTable();
>> 591
>> 592 size_t index = couple->GetIndex();
>> 593 G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
>> 594 G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
>> 595
>> 596 // Generation of fluorescence
>> 597 // Data in EADL are available only for Z > 5
>> 598 // Protection to avoid generating photons in the unphysical case of
>> 599 // shell binding energy > photon energy
>> 600 if (localEnergyDeposit > cutg || localEnergyDeposit > cute)
>> 601 {
>> 602 G4DynamicParticle* aPhoton;
>> 603 deexcitationManager.SetCutForSecondaryPhotons(cutg);
>> 604 deexcitationManager.SetCutForAugerElectrons(cute);
>> 605
>> 606 photonVector = deexcitationManager.GenerateParticles(Z,shellId);
>> 607 if(photonVector)
>> 608 {
>> 609 size_t nPhotons = photonVector->size();
>> 610 for (size_t k=0; k<nPhotons; k++)
>> 611 {
>> 612 aPhoton = (*photonVector)[k];
>> 613 if (aPhoton)
>> 614 {
>> 615 G4double itsEnergy = aPhoton->GetKineticEnergy();
>> 616 G4bool keepIt = false;
>> 617 if (itsEnergy <= localEnergyDeposit)
>> 618 {
>> 619 //check if good!
>> 620 if(aPhoton->GetDefinition() == G4Gamma::Gamma()
>> 621 && itsEnergy >= cutg)
>> 622 {
>> 623 keepIt = true;
>> 624 energyInFluorescence += itsEnergy;
>> 625 }
>> 626 if (aPhoton->GetDefinition() == G4Electron::Electron() &&
>> 627 itsEnergy >= cute)
>> 628 {
>> 629 energyInAuger += itsEnergy;
>> 630 keepIt = true;
>> 631 }
>> 632 }
>> 633 //good secondary, register it
>> 634 if (keepIt)
>> 635 {
>> 636 localEnergyDeposit -= itsEnergy;
>> 637 fvect->push_back(aPhoton);
>> 638 }
>> 639 else
>> 640 {
>> 641 delete aPhoton;
>> 642 (*photonVector)[k] = 0;
>> 643 }
>> 644 }
>> 645 }
>> 646 delete photonVector;
>> 647 }
>> 648 }
>> 649 }
>> 650 */
>> 651
>> 652
>> 653 //Always produce explicitely the electron
597 G4DynamicParticle* electron = 0; 654 G4DynamicParticle* electron = 0;
598 655
599 G4double xEl = sinThetaE * std::cos(phi+pi); 656 G4double xEl = sinThetaE * std::cos(phi+pi);
600 G4double yEl = sinThetaE * std::sin(phi+pi); 657 G4double yEl = sinThetaE * std::sin(phi+pi);
601 G4double zEl = cosThetaE; 658 G4double zEl = cosThetaE;
602 G4ThreeVector eDirection(xEl,yEl,zEl); //ele 659 G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction
603 eDirection.rotateUz(photonDirection0); 660 eDirection.rotateUz(photonDirection0);
604 electron = new G4DynamicParticle (G4Electron 661 electron = new G4DynamicParticle (G4Electron::Electron(),
605 eDirection,eKineticEnergy) ; 662 eDirection,eKineticEnergy) ;
606 fvect->push_back(electron); 663 fvect->push_back(electron);
607 664
608 if (localEnergyDeposit < 0) //Should not be: << 665
>> 666 if (localEnergyDeposit < 0)
609 { 667 {
610 G4Exception("G4PenelopeComptonModel::Sam << 668 G4cout << "WARNING-"
611 "em2099",JustWarning,"WARNING: Negative << 669 << "G4PenelopeComptonModel::SampleSecondaries - Negative energy deposit"
>> 670 << G4endl;
612 localEnergyDeposit=0.; 671 localEnergyDeposit=0.;
613 } 672 }
614 fParticleChange->ProposeLocalEnergyDeposit(l 673 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
615 674
616 G4double electronEnergy = 0.; 675 G4double electronEnergy = 0.;
617 if (electron) 676 if (electron)
618 electronEnergy = eKineticEnergy; 677 electronEnergy = eKineticEnergy;
619 if (fVerboseLevel > 1) << 678 if (verboseLevel > 1)
620 { 679 {
621 G4cout << "----------------------------- 680 G4cout << "-----------------------------------------------------------" << G4endl;
622 G4cout << "Energy balance from G4Penelop 681 G4cout << "Energy balance from G4PenelopeCompton" << G4endl;
623 G4cout << "Incoming photon energy: " << 682 G4cout << "Incoming photon energy: " << photonEnergy0/keV << " keV" << G4endl;
624 G4cout << "----------------------------- 683 G4cout << "-----------------------------------------------------------" << G4endl;
625 G4cout << "Scattered photon: " << photon 684 G4cout << "Scattered photon: " << photonEnergy1/keV << " keV" << G4endl;
626 G4cout << "Scattered electron " << elect 685 G4cout << "Scattered electron " << electronEnergy/keV << " keV" << G4endl;
627 if (energyInFluorescence) 686 if (energyInFluorescence)
628 G4cout << "Fluorescence x-rays: " << energyI 687 G4cout << "Fluorescence x-rays: " << energyInFluorescence/keV << " keV" << G4endl;
629 if (energyInAuger) 688 if (energyInAuger)
630 G4cout << "Auger electrons: " << energyInAug 689 G4cout << "Auger electrons: " << energyInAuger/keV << " keV" << G4endl;
631 G4cout << "Local energy deposit " << loc 690 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
632 G4cout << "Total final state: " << (phot 691 G4cout << "Total final state: " << (photonEnergy1+electronEnergy+energyInFluorescence+
633 localEnergyDeposit+energyInAuger)/ 692 localEnergyDeposit+energyInAuger)/keV <<
634 " keV" << G4endl; 693 " keV" << G4endl;
635 G4cout << "----------------------------- 694 G4cout << "-----------------------------------------------------------" << G4endl;
636 } 695 }
637 if (fVerboseLevel > 0) << 696 if (verboseLevel > 0)
638 { 697 {
639 G4double energyDiff = std::fabs(photonEn 698 G4double energyDiff = std::fabs(photonEnergy1+
640 electronEnergy+energyInFluoresce 699 electronEnergy+energyInFluorescence+
641 localEnergyDeposit+energyInAuger 700 localEnergyDeposit+energyInAuger-photonEnergy0);
642 if (energyDiff > 0.05*keV) 701 if (energyDiff > 0.05*keV)
643 G4cout << "Warning from G4PenelopeCompton: p 702 G4cout << "Warning from G4PenelopeCompton: problem with energy conservation: " <<
644 (photonEnergy1+electronEnergy+energyInFluo 703 (photonEnergy1+electronEnergy+energyInFluorescence+energyInAuger+localEnergyDeposit)/keV <<
645 " keV (final) vs. " << 704 " keV (final) vs. " <<
646 photonEnergy0/keV << " keV (initial)" << G 705 photonEnergy0/keV << " keV (initial)" << G4endl;
647 } 706 }
648 } 707 }
649 708
650 //....oooOO0OOooo........oooOO0OOooo........oo 709 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
651 710
652 G4double G4PenelopeComptonModel::DifferentialC 711 G4double G4PenelopeComptonModel::DifferentialCrossSection(G4double cosTheta,G4double energy,
653 G4PenelopeOscillator* osc) 712 G4PenelopeOscillator* osc)
654 { 713 {
655 // 714 //
656 // Penelope model v2008. Single differential 715 // Penelope model v2008. Single differential cross section *per electron*
657 // for photon Compton scattering by 716 // for photon Compton scattering by
658 // electrons in the given atomic oscillator, 717 // electrons in the given atomic oscillator, differential in the direction of the
659 // scattering photon. This is in units of pi 718 // scattering photon. This is in units of pi*classic_electr_radius**2
660 // 719 //
661 // D. Brusa et al., Nucl. Instrum. Meth. A 3 720 // D. Brusa et al., Nucl. Instrum. Meth. A 379 (1996) 167
662 // The parametrization includes the J(p) dis 721 // The parametrization includes the J(p) distribution profiles for the atomic shells,
663 // that are tabulated from Hartree-Fock calc 722 // that are tabulated from Hartree-Fock calculations
664 // from F. Biggs et al., At. Data Nucl. Data 723 // from F. Biggs et al., At. Data Nucl. Data Tables 16 (1975) 201
665 // 724 //
666 G4double ionEnergy = osc->GetIonisationEnerg 725 G4double ionEnergy = osc->GetIonisationEnergy();
667 G4double harFunc = osc->GetHartreeFactor(); 726 G4double harFunc = osc->GetHartreeFactor();
668 727
669 static const G4double k2 = std::sqrt(2.); 728 static const G4double k2 = std::sqrt(2.);
670 static const G4double k1 = 1./k2; 729 static const G4double k1 = 1./k2;
671 730
672 if (energy < ionEnergy) 731 if (energy < ionEnergy)
673 return 0; 732 return 0;
674 733
675 //energy of the Compton line 734 //energy of the Compton line
676 G4double cdt1 = 1.0-cosTheta; 735 G4double cdt1 = 1.0-cosTheta;
677 G4double EOEC = 1.0+(energy/electron_mass_c2 736 G4double EOEC = 1.0+(energy/electron_mass_c2)*cdt1;
678 G4double ECOE = 1.0/EOEC; 737 G4double ECOE = 1.0/EOEC;
679 738
680 //Incoherent scattering function (analytical 739 //Incoherent scattering function (analytical profile)
681 G4double aux = energy*(energy-ionEnergy)*cdt 740 G4double aux = energy*(energy-ionEnergy)*cdt1;
682 G4double Pzimax = 741 G4double Pzimax =
683 (aux - electron_mass_c2*ionEnergy)/(electr 742 (aux - electron_mass_c2*ionEnergy)/(electron_mass_c2*std::sqrt(2*aux+ionEnergy*ionEnergy));
684 G4double sia = 0.0; 743 G4double sia = 0.0;
685 G4double x = harFunc*Pzimax; 744 G4double x = harFunc*Pzimax;
686 if (x > 0) 745 if (x > 0)
687 sia = 1.0-0.5*G4Exp(0.5-(k1+k2*x)*(k1+k2*x 746 sia = 1.0-0.5*G4Exp(0.5-(k1+k2*x)*(k1+k2*x));
688 else 747 else
689 sia = 0.5*G4Exp(0.5-(k1-k2*x)*(k1-k2*x)); 748 sia = 0.5*G4Exp(0.5-(k1-k2*x)*(k1-k2*x));
690 749
691 //1st order correction, integral of Pz times 750 //1st order correction, integral of Pz times the Compton profile.
692 //Calculated approximately using a free-elec 751 //Calculated approximately using a free-electron gas profile
693 G4double pf = 3.0/(4.0*harFunc); 752 G4double pf = 3.0/(4.0*harFunc);
694 if (std::fabs(Pzimax) < pf) 753 if (std::fabs(Pzimax) < pf)
695 { 754 {
696 G4double QCOE2 = 1.0+ECOE*ECOE-2.0*ECOE* 755 G4double QCOE2 = 1.0+ECOE*ECOE-2.0*ECOE*cosTheta;
697 G4double p2 = Pzimax*Pzimax; 756 G4double p2 = Pzimax*Pzimax;
698 G4double dspz = std::sqrt(QCOE2)* 757 G4double dspz = std::sqrt(QCOE2)*
699 (1.0+ECOE*(ECOE-cosTheta)/QCOE2)*harFunc 758 (1.0+ECOE*(ECOE-cosTheta)/QCOE2)*harFunc
700 *0.25*(2*p2-(p2*p2)/(pf*pf)-(pf*pf)); 759 *0.25*(2*p2-(p2*p2)/(pf*pf)-(pf*pf));
701 sia += std::max(dspz,-1.0*sia); 760 sia += std::max(dspz,-1.0*sia);
702 } 761 }
703 762
704 G4double XKN = EOEC+ECOE-1.0+cosTheta*cosThe 763 G4double XKN = EOEC+ECOE-1.0+cosTheta*cosTheta;
705 764
706 //Differential cross section (per electron, 765 //Differential cross section (per electron, in units of pi*classic_electr_radius**2)
707 G4double diffCS = ECOE*ECOE*XKN*sia; 766 G4double diffCS = ECOE*ECOE*XKN*sia;
708 767
709 return diffCS; 768 return diffCS;
710 } 769 }
711 770
712 //....oooOO0OOooo........oooOO0OOooo........oo 771 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
713 772
714 G4double G4PenelopeComptonModel::OscillatorTot 773 G4double G4PenelopeComptonModel::OscillatorTotalCrossSection(G4double energy,G4PenelopeOscillator* osc)
715 { 774 {
716 //Total cross section (integrated) for the g 775 //Total cross section (integrated) for the given oscillator in units of
717 //pi*classic_electr_radius^2 776 //pi*classic_electr_radius^2
718 777
719 //Integrate differential cross section for e 778 //Integrate differential cross section for each oscillator
720 G4double stre = osc->GetOscillatorStrength() 779 G4double stre = osc->GetOscillatorStrength();
721 780
722 // here one uses the using the 20-point 781 // here one uses the using the 20-point
723 // Gauss quadrature method with an adaptive 782 // Gauss quadrature method with an adaptive bipartition scheme
724 const G4int npoints=10; 783 const G4int npoints=10;
725 const G4int ncallsmax=20000; 784 const G4int ncallsmax=20000;
726 const G4int nst=256; 785 const G4int nst=256;
727 static const G4double Abscissas[10] = {7.652 786 static const G4double Abscissas[10] = {7.652651133497334e-02,2.2778585114164508e-01,3.7370608871541956e-01,
728 5.1086700195082710e-01,6.36053680726 787 5.1086700195082710e-01,6.3605368072651503e-01,7.4633190646015079e-01,
729 8.3911697182221882e-01,9.12234428251 788 8.3911697182221882e-01,9.1223442825132591e-01,9.6397192727791379e-01,
730 9.9312859918509492e-01}; 789 9.9312859918509492e-01};
731 static const G4double Weights[10] = {1.52753 790 static const G4double Weights[10] = {1.5275338713072585e-01,1.4917298647260375e-01,1.4209610931838205e-01,
732 1.3168863844917663e-01,1.1819453196151 791 1.3168863844917663e-01,1.1819453196151842e-01,1.0193011981724044e-01,
733 8.3276741576704749e-02,6.2672048334109 792 8.3276741576704749e-02,6.2672048334109064e-02,4.0601429800386941e-02,
734 1.7614007139152118e-02}; 793 1.7614007139152118e-02};
735 794
736 G4double MaxError = 1e-5; 795 G4double MaxError = 1e-5;
737 //Error control 796 //Error control
738 G4double Ctol = std::min(std::max(MaxError,1 797 G4double Ctol = std::min(std::max(MaxError,1e-13),1e-02);
739 G4double Ptol = 0.01*Ctol; 798 G4double Ptol = 0.01*Ctol;
740 G4double Err=1e35; 799 G4double Err=1e35;
741 800
742 //Gauss integration from -1 to 1 801 //Gauss integration from -1 to 1
743 G4double LowPoint = -1.0; 802 G4double LowPoint = -1.0;
744 G4double HighPoint = 1.0; 803 G4double HighPoint = 1.0;
745 804
746 G4double h=HighPoint-LowPoint; 805 G4double h=HighPoint-LowPoint;
747 G4double sumga=0.0; 806 G4double sumga=0.0;
748 G4double a=0.5*(HighPoint-LowPoint); 807 G4double a=0.5*(HighPoint-LowPoint);
749 G4double b=0.5*(HighPoint+LowPoint); 808 G4double b=0.5*(HighPoint+LowPoint);
750 G4double c=a*Abscissas[0]; 809 G4double c=a*Abscissas[0];
751 G4double d= Weights[0]* 810 G4double d= Weights[0]*
752 (DifferentialCrossSection(b+c,energy,osc)+ 811 (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
753 for (G4int i=2;i<=npoints;i++) 812 for (G4int i=2;i<=npoints;i++)
754 { 813 {
755 c=a*Abscissas[i-1]; 814 c=a*Abscissas[i-1];
756 d += Weights[i-1]* 815 d += Weights[i-1]*
757 (DifferentialCrossSection(b+c,energy,osc)+Di 816 (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
758 } 817 }
759 G4int icall = 2*npoints; 818 G4int icall = 2*npoints;
760 G4int LH=1; 819 G4int LH=1;
761 G4double S[nst],x[nst],sn[nst],xrn[nst]; 820 G4double S[nst],x[nst],sn[nst],xrn[nst];
762 S[0]=d*a; 821 S[0]=d*a;
763 x[0]=LowPoint; 822 x[0]=LowPoint;
764 823
765 G4bool loopAgain = true; 824 G4bool loopAgain = true;
766 825
767 //Adaptive bipartition scheme 826 //Adaptive bipartition scheme
768 do{ 827 do{
769 G4double h0=h; 828 G4double h0=h;
770 h=0.5*h; //bipartition 829 h=0.5*h; //bipartition
771 G4double sumr=0; 830 G4double sumr=0;
772 G4int LHN=0; 831 G4int LHN=0;
773 G4double si,xa,xb,xc; 832 G4double si,xa,xb,xc;
774 for (G4int i=1;i<=LH;i++){ 833 for (G4int i=1;i<=LH;i++){
775 si=S[i-1]; 834 si=S[i-1];
776 xa=x[i-1]; 835 xa=x[i-1];
777 xb=xa+h; 836 xb=xa+h;
778 xc=xa+h0; 837 xc=xa+h0;
779 a=0.5*(xb-xa); 838 a=0.5*(xb-xa);
780 b=0.5*(xb+xa); 839 b=0.5*(xb+xa);
781 c=a*Abscissas[0]; 840 c=a*Abscissas[0];
782 G4double dLocal = Weights[0]* 841 G4double dLocal = Weights[0]*
783 (DifferentialCrossSection(b+c,energy,osc)+Di 842 (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
784 843
785 for (G4int j=1;j<npoints;j++) 844 for (G4int j=1;j<npoints;j++)
786 { 845 {
787 c=a*Abscissas[j]; 846 c=a*Abscissas[j];
788 dLocal += Weights[j]* 847 dLocal += Weights[j]*
789 (DifferentialCrossSection(b+c,energy,osc 848 (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
790 } 849 }
791 G4double s1=dLocal*a; 850 G4double s1=dLocal*a;
792 a=0.5*(xc-xb); 851 a=0.5*(xc-xb);
793 b=0.5*(xc+xb); 852 b=0.5*(xc+xb);
794 c=a*Abscissas[0]; 853 c=a*Abscissas[0];
795 dLocal=Weights[0]* 854 dLocal=Weights[0]*
796 (DifferentialCrossSection(b+c,energy,osc)+Di 855 (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
797 856
798 for (G4int j=1;j<npoints;j++) 857 for (G4int j=1;j<npoints;j++)
799 { 858 {
800 c=a*Abscissas[j]; 859 c=a*Abscissas[j];
801 dLocal += Weights[j]* 860 dLocal += Weights[j]*
802 (DifferentialCrossSection(b+c,energy,osc 861 (DifferentialCrossSection(b+c,energy,osc)+DifferentialCrossSection(b-c,energy,osc));
803 } 862 }
804 G4double s2=dLocal*a; 863 G4double s2=dLocal*a;
805 icall=icall+4*npoints; 864 icall=icall+4*npoints;
806 G4double s12=s1+s2; 865 G4double s12=s1+s2;
807 if (std::abs(s12-si)<std::max(Ptol*std:: 866 if (std::abs(s12-si)<std::max(Ptol*std::abs(s12),1e-35))
808 sumga += s12; 867 sumga += s12;
809 else 868 else
810 { 869 {
811 sumr += s12; 870 sumr += s12;
812 LHN += 2; 871 LHN += 2;
813 sn[LHN-1]=s2; 872 sn[LHN-1]=s2;
814 xrn[LHN-1]=xb; 873 xrn[LHN-1]=xb;
815 sn[LHN-2]=s1; 874 sn[LHN-2]=s1;
816 xrn[LHN-2]=xa; 875 xrn[LHN-2]=xa;
817 } 876 }
818 877
819 if (icall>ncallsmax || LHN>nst) 878 if (icall>ncallsmax || LHN>nst)
820 { 879 {
821 G4cout << "G4PenelopeComptonModel: " << G4 880 G4cout << "G4PenelopeComptonModel: " << G4endl;
822 G4cout << "LowPoint: " << LowPoint << ", H 881 G4cout << "LowPoint: " << LowPoint << ", High Point: " << HighPoint << G4endl;
823 G4cout << "Tolerance: " << MaxError << G4e 882 G4cout << "Tolerance: " << MaxError << G4endl;
824 G4cout << "Calls: " << icall << ", Integra 883 G4cout << "Calls: " << icall << ", Integral: " << sumga << ", Error: " << Err << G4endl;
825 G4cout << "Number of open subintervals: " 884 G4cout << "Number of open subintervals: " << LHN << G4endl;
826 G4cout << "WARNING: the required accuracy 885 G4cout << "WARNING: the required accuracy has not been attained" << G4endl;
827 loopAgain = false; 886 loopAgain = false;
828 } 887 }
829 } 888 }
830 Err=std::abs(sumr)/std::max(std::abs(sumr+ 889 Err=std::abs(sumr)/std::max(std::abs(sumr+sumga),1e-35);
831 if (Err < Ctol || LHN == 0) 890 if (Err < Ctol || LHN == 0)
832 loopAgain = false; //end of cycle 891 loopAgain = false; //end of cycle
833 LH=LHN; 892 LH=LHN;
834 for (G4int i=0;i<LH;i++) 893 for (G4int i=0;i<LH;i++)
835 { 894 {
836 S[i]=sn[i]; 895 S[i]=sn[i];
837 x[i]=xrn[i]; 896 x[i]=xrn[i];
838 } 897 }
839 }while(Ctol < 1.0 && loopAgain); 898 }while(Ctol < 1.0 && loopAgain);
840 899
>> 900
841 G4double xs = stre*sumga; 901 G4double xs = stre*sumga;
842 902
843 return xs; 903 return xs;
844 } 904 }
845 905
846 //....oooOO0OOooo........oooOO0OOooo........oo 906 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
847 907
848 G4double G4PenelopeComptonModel::KleinNishinaC 908 G4double G4PenelopeComptonModel::KleinNishinaCrossSection(G4double energy,
849 const G4Material* material) 909 const G4Material* material)
850 { 910 {
851 // use Klein-Nishina formula 911 // use Klein-Nishina formula
852 // total cross section in units of pi*classi 912 // total cross section in units of pi*classic_electr_radius^2
>> 913
853 G4double cs = 0; 914 G4double cs = 0;
854 915
855 G4double ek =energy/electron_mass_c2; 916 G4double ek =energy/electron_mass_c2;
856 G4double eks = ek*ek; 917 G4double eks = ek*ek;
857 G4double ek2 = 1.0+ek+ek; 918 G4double ek2 = 1.0+ek+ek;
858 G4double ek1 = eks-ek2-1.0; 919 G4double ek1 = eks-ek2-1.0;
859 920
860 G4double t0 = 1.0/ek2; 921 G4double t0 = 1.0/ek2;
861 G4double csl = 0.5*eks*t0*t0+ek2*t0+ek1*G4Lo << 922 G4double csl = 0.5*eks*t0*t0+ek2*t0+ek1*std::log(t0)-(1.0/t0);
862 923
863 G4PenelopeOscillatorTable* theTable = fOscMa << 924 G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableCompton(material);
864 925
865 for (size_t i=0;i<theTable->size();i++) 926 for (size_t i=0;i<theTable->size();i++)
866 { 927 {
867 G4PenelopeOscillator* theOsc = (*theTabl 928 G4PenelopeOscillator* theOsc = (*theTable)[i];
868 G4double ionEnergy = theOsc->GetIonisati 929 G4double ionEnergy = theOsc->GetIonisationEnergy();
869 G4double tau=(energy-ionEnergy)/energy; 930 G4double tau=(energy-ionEnergy)/energy;
870 if (tau > t0) 931 if (tau > t0)
871 { 932 {
872 G4double csu = 0.5*eks*tau*tau+ek2*tau+ek1 << 933 G4double csu = 0.5*eks*tau*tau+ek2*tau+ek1*std::log(tau)-(1.0/tau);
873 G4double stre = theOsc->GetOscillatorStren 934 G4double stre = theOsc->GetOscillatorStrength();
874 935
875 cs += stre*(csu-csl); 936 cs += stre*(csu-csl);
876 } 937 }
877 } 938 }
>> 939
878 cs /= (ek*eks); 940 cs /= (ek*eks);
879 941
880 return cs; 942 return cs;
881 943
882 } 944 }
883 945
884 //....oooOO0OOooo........oooOO0OOooo........oo 946 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
885 947
886 void G4PenelopeComptonModel::SetParticle(const 948 void G4PenelopeComptonModel::SetParticle(const G4ParticleDefinition* p)
887 { 949 {
888 if(!fParticle) { 950 if(!fParticle) {
889 fParticle = p; 951 fParticle = p;
890 } 952 }
891 } 953 }
892 954