| 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 // ------------------------------------------- 27 // -------------------------------------------------------------------
28 // 28 //
29 // Geant4 Class file 29 // Geant4 Class file
30 // 30 //
31 // Authors: Alfonso Mantero (Alfonso.Mantero@g 31 // Authors: Alfonso Mantero (Alfonso.Mantero@ge.infn.it)
32 // 32 //
33 // Created 22 April 2010 from old G4UAtomicDee 33 // Created 22 April 2010 from old G4UAtomicDeexcitation class
34 // 34 //
35 // Modified: 35 // Modified:
36 // --------- 36 // ---------
37 // 20 Oct 2011 Alf modified to take into acc 37 // 20 Oct 2011 Alf modified to take into account ECPSSR form Form Factor
38 // 03 Nov 2011 Alf Extended Empirical and Fo 38 // 03 Nov 2011 Alf Extended Empirical and Form Factor ionisation XS models
39 // out thei ranges with Anal 39 // out thei ranges with Analytical one.
40 // 07 Nov 2011 Alf Restored original ioniati 40 // 07 Nov 2011 Alf Restored original ioniation XS for alphas,
41 // letting scaled ones for o 41 // letting scaled ones for other ions.
42 // 20 Mar 2012 LP Register G4PenelopeIonisa 42 // 20 Mar 2012 LP Register G4PenelopeIonisationCrossSection
43 // 43 //
44 // ------------------------------------------- 44 // -------------------------------------------------------------------
45 // 45 //
46 // Class description: 46 // Class description:
47 // Implementation of atomic deexcitation 47 // Implementation of atomic deexcitation
48 // 48 //
49 // ------------------------------------------- 49 // -------------------------------------------------------------------
50 50
51 #include "G4UAtomicDeexcitation.hh" 51 #include "G4UAtomicDeexcitation.hh"
52 #include "G4PhysicalConstants.hh" 52 #include "G4PhysicalConstants.hh"
53 #include "G4SystemOfUnits.hh" 53 #include "G4SystemOfUnits.hh"
54 #include "Randomize.hh" 54 #include "Randomize.hh"
55 #include "G4Gamma.hh" 55 #include "G4Gamma.hh"
56 #include "G4AtomicTransitionManager.hh" 56 #include "G4AtomicTransitionManager.hh"
57 #include "G4FluoTransition.hh" 57 #include "G4FluoTransition.hh"
58 #include "G4Electron.hh" 58 #include "G4Electron.hh"
59 #include "G4Positron.hh" 59 #include "G4Positron.hh"
60 #include "G4Proton.hh" 60 #include "G4Proton.hh"
61 #include "G4Alpha.hh" 61 #include "G4Alpha.hh"
62 62
63 #include "G4teoCrossSection.hh" 63 #include "G4teoCrossSection.hh"
64 #include "G4empCrossSection.hh" 64 #include "G4empCrossSection.hh"
65 #include "G4PenelopeIonisationCrossSection.hh" 65 #include "G4PenelopeIonisationCrossSection.hh"
66 #include "G4LivermoreIonisationCrossSection.hh 66 #include "G4LivermoreIonisationCrossSection.hh"
67 #include "G4EmCorrections.hh" 67 #include "G4EmCorrections.hh"
68 #include "G4LossTableManager.hh" 68 #include "G4LossTableManager.hh"
69 #include "G4EmParameters.hh" 69 #include "G4EmParameters.hh"
70 #include "G4Material.hh" 70 #include "G4Material.hh"
71 #include "G4AtomicShells.hh" 71 #include "G4AtomicShells.hh"
72 72
73 using namespace std; 73 using namespace std;
74 74
75 //....oooOO0OOooo........oooOO0OOooo........oo 75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76 76
77 G4UAtomicDeexcitation::G4UAtomicDeexcitation() 77 G4UAtomicDeexcitation::G4UAtomicDeexcitation():
78 G4VAtomDeexcitation("UAtomDeexcitation"), 78 G4VAtomDeexcitation("UAtomDeexcitation"),
79 minGammaEnergy(DBL_MAX), 79 minGammaEnergy(DBL_MAX),
80 minElectronEnergy(DBL_MAX), 80 minElectronEnergy(DBL_MAX),
81 newShellId(-1) 81 newShellId(-1)
82 { 82 {
83 anaPIXEshellCS = nullptr; 83 anaPIXEshellCS = nullptr;
84 PIXEshellCS = nullptr; 84 PIXEshellCS = nullptr;
85 ePIXEshellCS = nullptr; 85 ePIXEshellCS = nullptr;
86 emcorr = G4LossTableManager::Instance()->EmC 86 emcorr = G4LossTableManager::Instance()->EmCorrections();
87 theElectron = G4Electron::Electron(); 87 theElectron = G4Electron::Electron();
88 thePositron = G4Positron::Positron(); 88 thePositron = G4Positron::Positron();
89 transitionManager = G4AtomicTransitionManage 89 transitionManager = G4AtomicTransitionManager::Instance();
90 } 90 }
91 91
92 //....oooOO0OOooo........oooOO0OOooo........oo 92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
93 93
94 G4UAtomicDeexcitation::~G4UAtomicDeexcitation( 94 G4UAtomicDeexcitation::~G4UAtomicDeexcitation()
95 { 95 {
96 delete anaPIXEshellCS; 96 delete anaPIXEshellCS;
97 delete PIXEshellCS; 97 delete PIXEshellCS;
98 delete ePIXEshellCS; 98 delete ePIXEshellCS;
99 } 99 }
100 100
101 //....oooOO0OOooo........oooOO0OOooo........oo 101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102 102
103 void G4UAtomicDeexcitation::InitialiseForNewRu 103 void G4UAtomicDeexcitation::InitialiseForNewRun()
104 { 104 {
105 if(!IsFluoActive()) { return; } 105 if(!IsFluoActive()) { return; }
106 transitionManager->Initialise(); 106 transitionManager->Initialise();
107 if(!IsPIXEActive()) { return; } 107 if(!IsPIXEActive()) { return; }
108 108
109 if(!anaPIXEshellCS) { 109 if(!anaPIXEshellCS) {
110 anaPIXEshellCS = new G4teoCrossSection("EC 110 anaPIXEshellCS = new G4teoCrossSection("ECPSSR_Analytical");
111 } 111 }
112 G4cout << G4endl; 112 G4cout << G4endl;
113 G4cout << "### === G4UAtomicDeexcitation::In 113 G4cout << "### === G4UAtomicDeexcitation::InitialiseForNewRun()" << G4endl;
114 114
115 G4EmParameters* param = G4EmParameters::Inst 115 G4EmParameters* param = G4EmParameters::Instance();
116 G4String namePIXExsModel = param->PIXECrossS 116 G4String namePIXExsModel = param->PIXECrossSectionModel();
117 G4String namePIXExsElectronModel = param->PI 117 G4String namePIXExsElectronModel = param->PIXEElectronCrossSectionModel();
118 118
119 // Check if old cross section for p/ion shou 119 // Check if old cross section for p/ion should be deleted
120 if(PIXEshellCS && namePIXExsModel != PIXEshe 120 if(PIXEshellCS && namePIXExsModel != PIXEshellCS->GetName())
121 { 121 {
122 delete PIXEshellCS; 122 delete PIXEshellCS;
123 PIXEshellCS = nullptr; 123 PIXEshellCS = nullptr;
124 } 124 }
125 125
126 // Instantiate new proton/ion cross section 126 // Instantiate new proton/ion cross section
127 if(!PIXEshellCS) { 127 if(!PIXEshellCS) {
128 if (namePIXExsModel == "ECPSSR_FormFactor" 128 if (namePIXExsModel == "ECPSSR_FormFactor")
129 { 129 {
130 PIXEshellCS = new G4teoCrossSection(namePIXE 130 PIXEshellCS = new G4teoCrossSection(namePIXExsModel);
131 } 131 }
132 else if(namePIXExsModel == "ECPSSR_ANSTO") 132 else if(namePIXExsModel == "ECPSSR_ANSTO")
133 { 133 {
134 PIXEshellCS = new G4teoCrossSection(namePIXE 134 PIXEshellCS = new G4teoCrossSection(namePIXExsModel);
135 } 135 }
136 else if(namePIXExsModel == "Empirical") 136 else if(namePIXExsModel == "Empirical")
137 { 137 {
138 PIXEshellCS = new G4empCrossSection(namePIXE 138 PIXEshellCS = new G4empCrossSection(namePIXExsModel);
139 } 139 }
140 } 140 }
141 141
142 // Check if old cross section for e+- should 142 // Check if old cross section for e+- should be deleted
143 if(ePIXEshellCS && namePIXExsElectronModel ! 143 if(ePIXEshellCS && namePIXExsElectronModel != ePIXEshellCS->GetName())
144 { 144 {
145 delete ePIXEshellCS; 145 delete ePIXEshellCS;
146 ePIXEshellCS = nullptr; 146 ePIXEshellCS = nullptr;
147 } 147 }
148 148
149 // Instantiate new e+- cross section 149 // Instantiate new e+- cross section
150 if(nullptr == ePIXEshellCS) 150 if(nullptr == ePIXEshellCS)
151 { 151 {
152 if(namePIXExsElectronModel == "Empirical 152 if(namePIXExsElectronModel == "Empirical")
153 { 153 {
154 ePIXEshellCS = new G4empCrossSection("Empi 154 ePIXEshellCS = new G4empCrossSection("Empirical");
155 } 155 }
156 else if(namePIXExsElectronModel == "ECPS 156 else if(namePIXExsElectronModel == "ECPSSR_Analytical")
157 { 157 {
158 ePIXEshellCS = new G4teoCrossSection("ECPS 158 ePIXEshellCS = new G4teoCrossSection("ECPSSR_Analytical");
159 } 159 }
160 else if (namePIXExsElectronModel == "Pen 160 else if (namePIXExsElectronModel == "Penelope")
161 { 161 {
162 ePIXEshellCS = new G4PenelopeIonisationCro 162 ePIXEshellCS = new G4PenelopeIonisationCrossSection();
163 } 163 }
164 else 164 else
165 { 165 {
166 ePIXEshellCS = new G4LivermoreIonisationCr 166 ePIXEshellCS = new G4LivermoreIonisationCrossSection();
167 } 167 }
168 } 168 }
169 } 169 }
170 170
171 //....oooOO0OOooo........oooOO0OOooo........oo 171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
172 172
173 void G4UAtomicDeexcitation::InitialiseForExtra 173 void G4UAtomicDeexcitation::InitialiseForExtraAtom(G4int)
174 {} 174 {}
175 175
176 //....oooOO0OOooo........oooOO0OOooo........oo 176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
177 177
178 const G4AtomicShell* 178 const G4AtomicShell*
179 G4UAtomicDeexcitation::GetAtomicShell(G4int Z, 179 G4UAtomicDeexcitation::GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)
180 { 180 {
181 return transitionManager->Shell(Z, (std::siz 181 return transitionManager->Shell(Z, (std::size_t)shell);
182 } 182 }
183 183
184 //....oooOO0OOooo........oooOO0OOooo........oo 184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
185 185
186 void G4UAtomicDeexcitation::GenerateParticles( 186 void G4UAtomicDeexcitation::GenerateParticles(
187 std::vector<G4DynamicParticle* 187 std::vector<G4DynamicParticle*>* vectorOfParticles,
188 const G4AtomicShell* atomicShell, 188 const G4AtomicShell* atomicShell,
189 G4int Z, 189 G4int Z,
190 G4double gammaCut, 190 G4double gammaCut,
191 G4double eCut) 191 G4double eCut)
192 { 192 {
193 // Defined initial conditions 193 // Defined initial conditions
194 G4int givenShellId = atomicShell->ShellId(); 194 G4int givenShellId = atomicShell->ShellId();
195 minGammaEnergy = gammaCut; 195 minGammaEnergy = gammaCut;
196 minElectronEnergy = eCut; 196 minElectronEnergy = eCut;
197 vacancyArray.clear(); << 197
198 <<
199 // generation secondaries 198 // generation secondaries
200 G4DynamicParticle* aParticle=0; 199 G4DynamicParticle* aParticle=0;
201 G4int provShellId = 0; 200 G4int provShellId = 0;
202 201
203 //ORIGINAL METHOD BY ALFONSO MANTERO 202 //ORIGINAL METHOD BY ALFONSO MANTERO
204 if (!IsAugerCascadeActive()) 203 if (!IsAugerCascadeActive())
205 { 204 {
206 //---------------------------- 205 //----------------------------
207 G4int counter = 0; 206 G4int counter = 0;
208 207
209 // limits of the EPDL data 208 // limits of the EPDL data
210 if (Z>5 && Z<105) { 209 if (Z>5 && Z<105) {
211 210
212 // The aim of this loop is to generate more 211 // The aim of this loop is to generate more than one fluorecence photon
213 // from the same ionizing event 212 // from the same ionizing event
214 do 213 do
215 { 214 {
216 if (counter == 0) 215 if (counter == 0)
217 // First call to GenerateParticles(... 216 // First call to GenerateParticles(...):
218 // givenShellId is given by the proces 217 // givenShellId is given by the process
219 { 218 {
220 provShellId = SelectTypeOfTransition(Z, gi 219 provShellId = SelectTypeOfTransition(Z, givenShellId);
221 220
222 if (provShellId >0) 221 if (provShellId >0)
223 { 222 {
224 aParticle = 223 aParticle =
225 GenerateFluorescence(Z, givenShellId 224 GenerateFluorescence(Z, givenShellId, provShellId);
226 } 225 }
227 else if (provShellId == -1) 226 else if (provShellId == -1)
228 { 227 {
229 aParticle = GenerateAuger(Z, givenShel 228 aParticle = GenerateAuger(Z, givenShellId);
230 } 229 }
231 } 230 }
232 else 231 else
233 // Following calls to GenerateParticle 232 // Following calls to GenerateParticles(...):
234 // newShellId is given by GenerateFluo 233 // newShellId is given by GenerateFluorescence(...)
235 { 234 {
236 provShellId = SelectTypeOfTransition(Z,new 235 provShellId = SelectTypeOfTransition(Z,newShellId);
237 if (provShellId >0) 236 if (provShellId >0)
238 { 237 {
239 aParticle = GenerateFluorescence(Z,new 238 aParticle = GenerateFluorescence(Z,newShellId,provShellId);
240 } 239 }
241 else if ( provShellId == -1) 240 else if ( provShellId == -1)
242 { 241 {
243 aParticle = GenerateAuger(Z, newShellI 242 aParticle = GenerateAuger(Z, newShellId);
244 } 243 }
245 } 244 }
246 ++counter; 245 ++counter;
247 if (aParticle != 0) 246 if (aParticle != 0)
248 { 247 {
249 vectorOfParticles->push_back(aParticle); 248 vectorOfParticles->push_back(aParticle);
250 } 249 }
251 else {provShellId = -2;} 250 else {provShellId = -2;}
252 } 251 }
253 while (provShellId > -2); 252 while (provShellId > -2);
254 } 253 }
255 } // Auger cascade is not active 254 } // Auger cascade is not active
256 255
257 //END OF ORIGINAL METHOD BY ALFONSO MANTERO 256 //END OF ORIGINAL METHOD BY ALFONSO MANTERO
258 //---------------------- 257 //----------------------
259 258
260 // NEW METHOD 259 // NEW METHOD
261 // Auger cascade by Burkhant Suerfu on March 260 // Auger cascade by Burkhant Suerfu on March 24 2015 (Bugzilla 1727)
262 if (IsAugerCascadeActive()) 261 if (IsAugerCascadeActive())
263 { 262 {
264 //---------------------- 263 //----------------------
265 vacancyArray.push_back(givenShellId); 264 vacancyArray.push_back(givenShellId);
266 265
267 // let's check that 5<Z<100 266 // let's check that 5<Z<100
268 if (Z<6 || Z>104){ 267 if (Z<6 || Z>104){
269 return; 268 return;
270 } 269 }
271 270
272 // as long as there is vacancy to be fil 271 // as long as there is vacancy to be filled by either fluo or auger, stay in the loop.
273 while(!vacancyArray.empty()){ 272 while(!vacancyArray.empty()){
274 // prepare to process the last element, and 273 // prepare to process the last element, and then delete it from the vector.
275 givenShellId = vacancyArray[0]; 274 givenShellId = vacancyArray[0];
276 provShellId = SelectTypeOfTransition(Z,given 275 provShellId = SelectTypeOfTransition(Z,givenShellId);
277 276
278 //G4cout<<"\n------ Atom Transition with Z: 277 //G4cout<<"\n------ Atom Transition with Z: "<<Z<<"\tbetween current:"
279 // <<givenShellId<<" & target:"<<provShel 278 // <<givenShellId<<" & target:"<<provShellId<<G4endl;
280 if(provShellId>0){ 279 if(provShellId>0){
281 aParticle = GenerateFluorescence(Z,givenSh 280 aParticle = GenerateFluorescence(Z,givenShellId,provShellId);
282 } 281 }
283 else if(provShellId == -1){ 282 else if(provShellId == -1){
284 aParticle = GenerateAuger(Z, givenShellId) 283 aParticle = GenerateAuger(Z, givenShellId);
285 } 284 }
286 // if a particle is created, put it in the 285 // if a particle is created, put it in the vector of new particles
287 if(aParticle!=0) 286 if(aParticle!=0)
288 vectorOfParticles->push_back(aParticle); 287 vectorOfParticles->push_back(aParticle);
289 288
290 // one vacancy has been processed. Erase it 289 // one vacancy has been processed. Erase it.
291 vacancyArray.erase(vacancyArray.begin()); 290 vacancyArray.erase(vacancyArray.begin());
292 } 291 }
293 //---------------------- 292 //----------------------
294 //End of Auger cascade by Burkhant Suerf 293 //End of Auger cascade by Burkhant Suerfu on March 24 2015 (Bugzilla 1727)
295 294
296 } // Auger cascade is active 295 } // Auger cascade is active
297 } 296 }
298 297
299 //....oooOO0OOooo........oooOO0OOooo........oo 298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
300 299
301 G4double 300 G4double
302 G4UAtomicDeexcitation::GetShellIonisationCross 301 G4UAtomicDeexcitation::GetShellIonisationCrossSectionPerAtom(
303 const G4ParticleDefinition* pdef, 302 const G4ParticleDefinition* pdef,
304 G4int Z, 303 G4int Z,
305 G4AtomicShellEnumerator shellEnum, 304 G4AtomicShellEnumerator shellEnum,
306 G4double kineticEnergy, 305 G4double kineticEnergy,
307 const G4Material* mat) 306 const G4Material* mat)
308 { 307 {
309 // we must put a control on the shell that a 308 // we must put a control on the shell that are passed:
310 // some shells should not pass (line "0" or 309 // some shells should not pass (line "0" or "2")
311 310
312 // check atomic number 311 // check atomic number
313 G4double xsec = 0.0; 312 G4double xsec = 0.0;
314 if(Z > 93 || Z < 6 ) { return xsec; } //corr 313 if(Z > 93 || Z < 6 ) { return xsec; } //corrected by alf - Z<6 missing
315 G4int idx = G4int(shellEnum); 314 G4int idx = G4int(shellEnum);
316 if(idx >= G4AtomicShells::GetNumberOfShells( 315 if(idx >= G4AtomicShells::GetNumberOfShells(Z)) { return xsec; }
317 316
318 if(pdef == theElectron || pdef == thePositro 317 if(pdef == theElectron || pdef == thePositron) {
319 xsec = ePIXEshellCS->CrossSection(Z,shellE 318 xsec = ePIXEshellCS->CrossSection(Z,shellEnum,kineticEnergy,0.0,mat);
320 return xsec; 319 return xsec;
321 } 320 }
322 321
323 G4double mass = pdef->GetPDGMass(); 322 G4double mass = pdef->GetPDGMass();
324 G4double escaled = kineticEnergy; 323 G4double escaled = kineticEnergy;
325 G4double q2 = 0.0; 324 G4double q2 = 0.0;
326 325
327 // scaling to protons for all particles excl 326 // scaling to protons for all particles excluding protons and alpha
328 G4int pdg = pdef->GetPDGEncoding(); 327 G4int pdg = pdef->GetPDGEncoding();
329 if (pdg != 2212 && pdg != 1000020040) 328 if (pdg != 2212 && pdg != 1000020040)
330 { 329 {
331 mass = proton_mass_c2; 330 mass = proton_mass_c2;
332 escaled = kineticEnergy*mass/(pdef->GetP 331 escaled = kineticEnergy*mass/(pdef->GetPDGMass());
333 332
334 if(mat) { 333 if(mat) {
335 q2 = emcorr->EffectiveChargeSquareRatio(pdef 334 q2 = emcorr->EffectiveChargeSquareRatio(pdef,mat,kineticEnergy);
336 } else { 335 } else {
337 G4double q = pdef->GetPDGCharge()/eplus; 336 G4double q = pdef->GetPDGCharge()/eplus;
338 q2 = q*q; 337 q2 = q*q;
339 } 338 }
340 } 339 }
341 340
342 if(PIXEshellCS) { 341 if(PIXEshellCS) {
343 xsec = PIXEshellCS->CrossSection(Z,shellEn 342 xsec = PIXEshellCS->CrossSection(Z,shellEnum,escaled,mass,mat);
344 } 343 }
345 if(xsec < 1e-100) { 344 if(xsec < 1e-100) {
346 xsec = anaPIXEshellCS->CrossSection(Z,shel 345 xsec = anaPIXEshellCS->CrossSection(Z,shellEnum,escaled,mass,mat);
347 } 346 }
348 347
349 if (q2) {xsec *= q2;} 348 if (q2) {xsec *= q2;}
350 349
351 return xsec; 350 return xsec;
352 } 351 }
353 352
354 //....oooOO0OOooo........oooOO0OOooo........oo 353 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
355 354
356 void G4UAtomicDeexcitation::SetCutForSecondary 355 void G4UAtomicDeexcitation::SetCutForSecondaryPhotons(G4double cut)
357 { 356 {
358 minGammaEnergy = cut; 357 minGammaEnergy = cut;
359 } 358 }
360 359
361 //....oooOO0OOooo........oooOO0OOooo........oo 360 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
362 361
363 void G4UAtomicDeexcitation::SetCutForAugerElec 362 void G4UAtomicDeexcitation::SetCutForAugerElectrons(G4double cut)
364 { 363 {
365 minElectronEnergy = cut; 364 minElectronEnergy = cut;
366 } 365 }
367 366
368 //....oooOO0OOooo........oooOO0OOooo........oo 367 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
369 368
370 G4double G4UAtomicDeexcitation::ComputeShellIo 369 G4double G4UAtomicDeexcitation::ComputeShellIonisationCrossSectionPerAtom(
371 const G4ParticleDefinition* p, 370 const G4ParticleDefinition* p,
372 G4int Z, 371 G4int Z,
373 G4AtomicShellEnumerator shell, 372 G4AtomicShellEnumerator shell,
374 G4double kinE, 373 G4double kinE,
375 const G4Material* mat) 374 const G4Material* mat)
376 { 375 {
377 return GetShellIonisationCrossSectionPerAtom 376 return GetShellIonisationCrossSectionPerAtom(p,Z,shell,kinE,mat);
378 } 377 }
379 378
380 //....oooOO0OOooo........oooOO0OOooo........oo 379 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
381 380
382 G4int G4UAtomicDeexcitation::SelectTypeOfTrans 381 G4int G4UAtomicDeexcitation::SelectTypeOfTransition(G4int Z, G4int shellId)
383 { 382 {
384 if (shellId <=0 ) { 383 if (shellId <=0 ) {
385 return 0; 384 return 0;
386 } 385 }
387 386
388 G4int provShellId = -1; 387 G4int provShellId = -1;
389 G4int shellNum = 0; 388 G4int shellNum = 0;
390 G4int maxNumOfShells = transitionManager->Nu 389 G4int maxNumOfShells = transitionManager->NumberOfReachableShells(Z);
391 390
392 const G4FluoTransition* refShell = 391 const G4FluoTransition* refShell =
393 transitionManager->ReachableShell(Z,maxNum 392 transitionManager->ReachableShell(Z,maxNumOfShells-1);
394 393
395 // This loop gives shellNum the value of the 394 // This loop gives shellNum the value of the index of shellId
396 // in the vector storing the list of the she 395 // in the vector storing the list of the shells reachable through
397 // a radiative transition 396 // a radiative transition
398 if ( shellId <= refShell->FinalShellId()) 397 if ( shellId <= refShell->FinalShellId())
399 { 398 {
400 while (shellId != transitionManager->Rea 399 while (shellId != transitionManager->ReachableShell(Z,shellNum)->FinalShellId())
401 { 400 {
402 if(shellNum ==maxNumOfShells-1) 401 if(shellNum ==maxNumOfShells-1)
403 { 402 {
404 break; 403 break;
405 } 404 }
406 shellNum++; 405 shellNum++;
407 } 406 }
408 G4int transProb = 0; //AM change 29/6/07 407 G4int transProb = 0; //AM change 29/6/07 was 1
409 408
410 G4double partialProb = G4UniformRand(); 409 G4double partialProb = G4UniformRand();
411 G4double partSum = 0; 410 G4double partSum = 0;
412 const G4FluoTransition* aShell = transit 411 const G4FluoTransition* aShell = transitionManager->ReachableShell(Z,shellNum);
413 G4int trSize = (G4int)(aShell->Transiti 412 G4int trSize = (G4int)(aShell->TransitionProbabilities()).size();
414 413
415 // Loop over the shells wich can provide 414 // Loop over the shells wich can provide an electron for a
416 // radiative transition towards shellId: 415 // radiative transition towards shellId:
417 // in every loop the partial sum of the 416 // in every loop the partial sum of the first transProb shells
418 // is calculated and compared with a ran 417 // is calculated and compared with a random number [0,1].
419 // If the partial sum is greater, the sh 418 // If the partial sum is greater, the shell whose index is transProb
420 // is chosen as the starting shell for a 419 // is chosen as the starting shell for a radiative transition
421 // and its identity is returned 420 // and its identity is returned
422 // Else, terminateded the loop, -1 is re 421 // Else, terminateded the loop, -1 is returned
423 while(transProb < trSize){ 422 while(transProb < trSize){
424 partSum += aShell->TransitionProbability(tra 423 partSum += aShell->TransitionProbability(transProb);
425 424
426 if(partialProb <= partSum) 425 if(partialProb <= partSum)
427 { 426 {
428 provShellId = aShell->OriginatingShellId 427 provShellId = aShell->OriginatingShellId(transProb);
429 break; 428 break;
430 } 429 }
431 ++transProb; 430 ++transProb;
432 } 431 }
433 // here provShellId is the right one or 432 // here provShellId is the right one or is -1.
434 // if -1, the control is passed to the A 433 // if -1, the control is passed to the Auger generation part of the package
435 } 434 }
436 else 435 else
437 { 436 {
438 provShellId = -1; 437 provShellId = -1;
439 } 438 }
440 return provShellId; 439 return provShellId;
441 } 440 }
442 441
443 //....oooOO0OOooo........oooOO0OOooo........oo 442 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
444 443
445 G4DynamicParticle* 444 G4DynamicParticle*
446 G4UAtomicDeexcitation::GenerateFluorescence(G4 445 G4UAtomicDeexcitation::GenerateFluorescence(G4int Z, G4int shellId,
447 G4int provShellId ) 446 G4int provShellId )
448 { 447 {
449 if (shellId <=0 ) 448 if (shellId <=0 )
450 { 449 {
451 return nullptr; 450 return nullptr;
452 } 451 }
453 452
454 //isotropic angular distribution for the out 453 //isotropic angular distribution for the outcoming photon
455 G4double newcosTh = 1.-2.*G4UniformRand(); 454 G4double newcosTh = 1.-2.*G4UniformRand();
456 G4double newsinTh = std::sqrt((1.-newcosTh)* 455 G4double newsinTh = std::sqrt((1.-newcosTh)*(1. + newcosTh));
457 G4double newPhi = twopi*G4UniformRand(); 456 G4double newPhi = twopi*G4UniformRand();
458 457
459 G4double xDir = newsinTh*std::sin(newPhi); 458 G4double xDir = newsinTh*std::sin(newPhi);
460 G4double yDir = newsinTh*std::cos(newPhi); 459 G4double yDir = newsinTh*std::cos(newPhi);
461 G4double zDir = newcosTh; 460 G4double zDir = newcosTh;
462 461
463 G4ThreeVector newGammaDirection(xDir,yDir,zD 462 G4ThreeVector newGammaDirection(xDir,yDir,zDir);
464 463
465 G4int shellNum = 0; 464 G4int shellNum = 0;
466 G4int maxNumOfShells = transitionManager->Nu 465 G4int maxNumOfShells = transitionManager->NumberOfReachableShells(Z);
467 466
468 // find the index of the shell named shellId 467 // find the index of the shell named shellId
469 while (shellId != transitionManager-> 468 while (shellId != transitionManager->
470 ReachableShell(Z,shellNum)->FinalShellId()) 469 ReachableShell(Z,shellNum)->FinalShellId())
471 { 470 {
472 if(shellNum == maxNumOfShells-1) 471 if(shellNum == maxNumOfShells-1)
473 { 472 {
474 break; 473 break;
475 } 474 }
476 ++shellNum; 475 ++shellNum;
477 } 476 }
478 // number of shell from wich an electron can 477 // number of shell from wich an electron can reach shellId
479 G4int transitionSize = (G4int)transitionMana 478 G4int transitionSize = (G4int)transitionManager->
480 ReachableShell(Z,shellNum)->OriginatingShe 479 ReachableShell(Z,shellNum)->OriginatingShellIds().size();
481 480
482 G4int index = 0; 481 G4int index = 0;
483 482
484 // find the index of the shell named provShe 483 // find the index of the shell named provShellId in the vector
485 // storing the shells from which shellId can 484 // storing the shells from which shellId can be reached
486 while (provShellId != transitionManager-> 485 while (provShellId != transitionManager->
487 ReachableShell(Z,shellNum)->OriginatingShel 486 ReachableShell(Z,shellNum)->OriginatingShellId(index))
488 { 487 {
489 if(index == transitionSize-1) 488 if(index == transitionSize-1)
490 { 489 {
491 break; 490 break;
492 } 491 }
493 ++index; 492 ++index;
494 } 493 }
495 // energy of the gamma leaving provShellId f 494 // energy of the gamma leaving provShellId for shellId
496 G4double transitionEnergy = transitionManage 495 G4double transitionEnergy = transitionManager->
497 ReachableShell(Z,shellNum)->TransitionEner 496 ReachableShell(Z,shellNum)->TransitionEnergy(index);
498 497
499 if (transitionEnergy < minGammaEnergy) retur 498 if (transitionEnergy < minGammaEnergy) return nullptr;
500 499
501 // This is the shell where the new vacancy i 500 // This is the shell where the new vacancy is: it is the same
502 // shell where the electron came from 501 // shell where the electron came from
503 newShellId = transitionManager-> 502 newShellId = transitionManager->
504 ReachableShell(Z,shellNum)->OriginatingShe 503 ReachableShell(Z,shellNum)->OriginatingShellId(index);
505 504
506 G4DynamicParticle* newPart = new G4DynamicPa 505 G4DynamicParticle* newPart = new G4DynamicParticle(G4Gamma::Gamma(),
507 newGammaDirection, 506 newGammaDirection,
508 transitionEnergy); 507 transitionEnergy);
509 508
510 //Auger cascade by Burkhant Suerfu on March 509 //Auger cascade by Burkhant Suerfu on March 24 2015 (Bugzilla 1727)
511 if (IsAugerCascadeActive()) vacancyArray.pus 510 if (IsAugerCascadeActive()) vacancyArray.push_back(newShellId);
512 511
513 return newPart; 512 return newPart;
514 } 513 }
515 514
516 //....oooOO0OOooo........oooOO0OOooo........oo 515 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
517 516
518 G4DynamicParticle* G4UAtomicDeexcitation::Gene 517 G4DynamicParticle* G4UAtomicDeexcitation::GenerateAuger(G4int Z, G4int shellId)
519 { 518 {
520 if(!IsAugerActive()) { 519 if(!IsAugerActive()) {
521 // G4cout << "auger inactive!" << G4end 520 // G4cout << "auger inactive!" << G4endl; //debug
522 return nullptr; 521 return nullptr;
523 } 522 }
524 523
525 if (shellId <=0 ) { 524 if (shellId <=0 ) {
526 //G4Exception("G4UAtomicDeexcitation::Gene 525 //G4Exception("G4UAtomicDeexcitation::GenerateAuger()","de0002",
527 // JustWarning, "Energy deposited local 526 // JustWarning, "Energy deposited locally");
528 return nullptr; 527 return nullptr;
529 } 528 }
530 529
531 G4int maxNumOfShells = transitionManager->Nu 530 G4int maxNumOfShells = transitionManager->NumberOfReachableAugerShells(Z);
532 531
533 const G4AugerTransition* refAugerTransition 532 const G4AugerTransition* refAugerTransition =
534 transitionManager->ReachableAugerShell(Z,m 533 transitionManager->ReachableAugerShell(Z,maxNumOfShells-1);
535 534
536 // This loop gives to shellNum the value of 535 // This loop gives to shellNum the value of the index of shellId
537 // in the vector storing the list of the vac 536 // in the vector storing the list of the vacancies in the variuos shells
538 // that can originate a NON-radiative transi 537 // that can originate a NON-radiative transition
539 G4int shellNum = 0; 538 G4int shellNum = 0;
540 539
541 if ( shellId <= refAugerTransition->FinalShe 540 if ( shellId <= refAugerTransition->FinalShellId() )
542 // "FinalShellId" is final from the point 541 // "FinalShellId" is final from the point of view of the electron
543 // who makes the transition, 542 // who makes the transition,
544 // being the Id of the shell in which ther 543 // being the Id of the shell in which there is a vacancy
545 { 544 {
546 G4int pippo = transitionManager->Reachab 545 G4int pippo = transitionManager->ReachableAugerShell(Z,shellNum)->FinalShellId();
547 if (shellId != pippo ) { 546 if (shellId != pippo ) {
548 do { 547 do {
549 ++shellNum; 548 ++shellNum;
550 if(shellNum == maxNumOfShells) 549 if(shellNum == maxNumOfShells)
551 { 550 {
552 // G4cout << "No Auger transition foun 551 // G4cout << "No Auger transition found" << G4endl; //debug
553 return 0; 552 return 0;
554 } 553 }
555 } 554 }
556 while (shellId != (transitionManager->Reacha 555 while (shellId != (transitionManager->ReachableAugerShell(Z,shellNum)->FinalShellId()) );
557 } 556 }
558 557
559 // Now we have that shellnum is the shel 558 // Now we have that shellnum is the shellIndex of the shell named ShellId
560 // G4cout << " the index of the she 559 // G4cout << " the index of the shell is: "<<shellNum<<G4endl;
561 // But we have now to select two shells: 560 // But we have now to select two shells: one for the transition,
562 // and another for the auger emission. 561 // and another for the auger emission.
563 G4int transitionLoopShellIndex = 0; 562 G4int transitionLoopShellIndex = 0;
564 G4double partSum = 0; 563 G4double partSum = 0;
565 const G4AugerTransition* anAugerTransiti 564 const G4AugerTransition* anAugerTransition =
566 transitionManager->ReachableAugerShell(Z,she 565 transitionManager->ReachableAugerShell(Z,shellNum);
567 566
568 G4int transitionSize = (G4int) 567 G4int transitionSize = (G4int)
569 (anAugerTransition->TransitionOriginatingShe 568 (anAugerTransition->TransitionOriginatingShellIds())->size();
570 while (transitionLoopShellIndex < transi 569 while (transitionLoopShellIndex < transitionSize) {
571 570
572 std::vector<G4int>::const_iterator pos 571 std::vector<G4int>::const_iterator pos =
573 anAugerTransition->TransitionOriginatingSh 572 anAugerTransition->TransitionOriginatingShellIds()->cbegin();
574 573
575 G4int transitionLoopShellId = *(pos+tr 574 G4int transitionLoopShellId = *(pos+transitionLoopShellIndex);
576 G4int numberOfPossibleAuger = (G4int) 575 G4int numberOfPossibleAuger = (G4int)
577 (anAugerTransition->AugerTransitionProbabi 576 (anAugerTransition->AugerTransitionProbabilities(transitionLoopShellId))->size();
578 G4int augerIndex = 0; 577 G4int augerIndex = 0;
579 578
580 if (augerIndex < numberOfPossibleAuger) { 579 if (augerIndex < numberOfPossibleAuger) {
581 do 580 do
582 { 581 {
583 G4double thisProb = anAugerTransition- 582 G4double thisProb = anAugerTransition->AugerTransitionProbability(augerIndex,
584 transitionLoopShellId); 583 transitionLoopShellId);
585 partSum += thisProb; 584 partSum += thisProb;
586 augerIndex++; 585 augerIndex++;
587 586
588 } while (augerIndex < numberOfPossibleAu 587 } while (augerIndex < numberOfPossibleAuger);
589 } 588 }
590 ++transitionLoopShellIndex; 589 ++transitionLoopShellIndex;
591 } 590 }
592 591
593 G4double totalVacancyAugerProbability = 592 G4double totalVacancyAugerProbability = partSum;
594 593
595 //And now we start to select the right a 594 //And now we start to select the right auger transition and emission
596 G4int transitionRandomShellIndex = 0; 595 G4int transitionRandomShellIndex = 0;
597 G4int transitionRandomShellId = 1; 596 G4int transitionRandomShellId = 1;
598 G4int augerIndex = 0; 597 G4int augerIndex = 0;
599 partSum = 0; 598 partSum = 0;
600 G4double partialProb = G4UniformRand(); 599 G4double partialProb = G4UniformRand();
601 600
602 G4int numberOfPossibleAuger = 0; 601 G4int numberOfPossibleAuger = 0;
603 G4bool foundFlag = false; 602 G4bool foundFlag = false;
604 603
605 while (transitionRandomShellIndex < tran 604 while (transitionRandomShellIndex < transitionSize) {
606 605
607 std::vector<G4int>::const_iterator pos 606 std::vector<G4int>::const_iterator pos =
608 anAugerTransition->TransitionOriginatingSh 607 anAugerTransition->TransitionOriginatingShellIds()->begin();
609 608
610 transitionRandomShellId = *(pos+transi 609 transitionRandomShellId = *(pos+transitionRandomShellIndex);
611 610
612 augerIndex = 0; 611 augerIndex = 0;
613 numberOfPossibleAuger = (G4int)(anAugerTrans 612 numberOfPossibleAuger = (G4int)(anAugerTransition->
614 AugerTransitionProbabilities(transiti 613 AugerTransitionProbabilities(transitionRandomShellId))->size();
615 614
616 while (augerIndex < numberOfPossibleAu 615 while (augerIndex < numberOfPossibleAuger) {
617 G4double thisProb =anAugerTransition->Auge 616 G4double thisProb =anAugerTransition->AugerTransitionProbability(augerIndex,
618 transitionRandomShellId); 617 transitionRandomShellId);
619 618
620 partSum += thisProb; 619 partSum += thisProb;
621 620
622 if (partSum >= (partialProb*totalVac 621 if (partSum >= (partialProb*totalVacancyAugerProbability) ) { // was /
623 foundFlag = true; 622 foundFlag = true;
624 break; 623 break;
625 } 624 }
626 augerIndex++; 625 augerIndex++;
627 } 626 }
628 if (partSum >= (partialProb*totalVacan 627 if (partSum >= (partialProb*totalVacancyAugerProbability) ) {break;} // was /
629 ++transitionRandomShellIndex; 628 ++transitionRandomShellIndex;
630 } 629 }
631 630
632 // Now we have the index of the shell fr 631 // Now we have the index of the shell from wich comes the auger electron (augerIndex),
633 // and the id of the shell, from which t 632 // and the id of the shell, from which the transition e- come (transitionRandomShellid)
634 // If no Transition has been found, 0 is 633 // If no Transition has been found, 0 is returned.
635 if (!foundFlag) { 634 if (!foundFlag) {
636 return nullptr; 635 return nullptr;
637 } 636 }
638 637
639 // Isotropic angular distribution for th 638 // Isotropic angular distribution for the outcoming e-
640 G4double newcosTh = 1.-2.*G4UniformRand( 639 G4double newcosTh = 1.-2.*G4UniformRand();
641 G4double newsinTh = std::sqrt(1.-newcosT 640 G4double newsinTh = std::sqrt(1.-newcosTh*newcosTh);
642 G4double newPhi = twopi*G4UniformRand(); 641 G4double newPhi = twopi*G4UniformRand();
643 642
644 G4double xDir = newsinTh*std::sin(newPhi 643 G4double xDir = newsinTh*std::sin(newPhi);
645 G4double yDir = newsinTh*std::cos(newPhi 644 G4double yDir = newsinTh*std::cos(newPhi);
646 G4double zDir = newcosTh; 645 G4double zDir = newcosTh;
647 646
648 G4ThreeVector newElectronDirection(xDir, 647 G4ThreeVector newElectronDirection(xDir,yDir,zDir);
649 648
650 // energy of the auger electron emitted 649 // energy of the auger electron emitted
651 G4double transitionEnergy = 650 G4double transitionEnergy =
652 anAugerTransition->AugerTransitionEnergy(aug 651 anAugerTransition->AugerTransitionEnergy(augerIndex, transitionRandomShellId);
653 652
654 if (transitionEnergy < minElectronEnergy 653 if (transitionEnergy < minElectronEnergy) {
655 return nullptr; 654 return nullptr;
656 } 655 }
657 656
658 // This is the shell where the new vacan 657 // This is the shell where the new vacancy is: it is the same
659 // shell where the electron came from 658 // shell where the electron came from
660 newShellId = transitionRandomShellId; 659 newShellId = transitionRandomShellId;
661 660
662 //Auger cascade by Burkhant Suerfu on Ma 661 //Auger cascade by Burkhant Suerfu on March 24 2015 (Bugzilla 1727)
663 if (IsAugerCascadeActive()) 662 if (IsAugerCascadeActive())
664 { 663 {
665 vacancyArray.push_back(newShellId); 664 vacancyArray.push_back(newShellId);
666 vacancyArray.push_back(anAugerTransition-> 665 vacancyArray.push_back(anAugerTransition->AugerOriginatingShellId(augerIndex,transitionRandomShellId));
667 } 666 }
668 667
669 return new G4DynamicParticle(G4Electron: 668 return new G4DynamicParticle(G4Electron::Electron(),
670 newElectronDirection, 669 newElectronDirection,
671 transitionEnergy); 670 transitionEnergy);
672 } 671 }
673 else 672 else
674 { 673 {
675 return nullptr; 674 return nullptr;
676 } 675 }
677 } 676 }
678 677