| 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: G4KleinNishinaModel.cc,v 1.4 2010/11/21 16:08:37 vnivanch Exp $
>> 27 // GEANT4 tag $Name: geant4-09-04 $
26 // 28 //
27 // ------------------------------------------- 29 // -------------------------------------------------------------------
28 // 30 //
29 // GEANT4 Class file 31 // GEANT4 Class file
30 // 32 //
31 // 33 //
32 // File name: G4KleinNishinaModel 34 // File name: G4KleinNishinaModel
33 // 35 //
34 // Author: Vladimir Ivanchenko on base 36 // Author: Vladimir Ivanchenko on base of G4KleinNishinaCompton
35 // 37 //
36 // Creation date: 13.06.2010 38 // Creation date: 13.06.2010
37 // 39 //
38 // Modifications: 40 // Modifications:
39 // 41 //
40 // Class Description: 42 // Class Description:
41 // 43 //
42 // ------------------------------------------- 44 // -------------------------------------------------------------------
43 // 45 //
44 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 48
47 #include "G4KleinNishinaModel.hh" 49 #include "G4KleinNishinaModel.hh"
48 #include "G4PhysicalConstants.hh" <<
49 #include "G4SystemOfUnits.hh" <<
50 #include "G4Electron.hh" 50 #include "G4Electron.hh"
51 #include "G4Gamma.hh" 51 #include "G4Gamma.hh"
52 #include "Randomize.hh" 52 #include "Randomize.hh"
53 #include "G4RandomDirection.hh" 53 #include "G4RandomDirection.hh"
54 #include "G4DataVector.hh" 54 #include "G4DataVector.hh"
55 #include "G4ParticleChangeForGamma.hh" 55 #include "G4ParticleChangeForGamma.hh"
56 #include "G4VAtomDeexcitation.hh" 56 #include "G4VAtomDeexcitation.hh"
57 #include "G4AtomicShells.hh" <<
58 #include "G4LossTableManager.hh" 57 #include "G4LossTableManager.hh"
59 #include "G4Log.hh" <<
60 #include "G4Exp.hh" <<
61 58
62 //....oooOO0OOooo........oooOO0OOooo........oo 59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63 60
64 using namespace std; 61 using namespace std;
65 62
66 G4KleinNishinaModel::G4KleinNishinaModel(const 63 G4KleinNishinaModel::G4KleinNishinaModel(const G4String& nam)
67 : G4VEmModel(nam), << 64 : G4VEmModel(nam),isInitialized(false)
68 lv1(0.,0.,0.,0.), <<
69 lv2(0.,0.,0.,0.), <<
70 bst(0.,0.,0.) <<
71 { 65 {
72 theGamma = G4Gamma::Gamma(); 66 theGamma = G4Gamma::Gamma();
73 theElectron = G4Electron::Electron(); 67 theElectron = G4Electron::Electron();
74 lowestSecondaryEnergy = 10*eV; << 68 lowestGammaEnergy = 1.0*eV;
75 limitFactor = 4; <<
76 fProbabilities.resize(9,0.0); 69 fProbabilities.resize(9,0.0);
77 SetDeexcitationFlag(true); <<
78 fParticleChange = nullptr; <<
79 fAtomDeexcitation = nullptr; <<
80 } 70 }
81 71
82 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83 73
84 G4KleinNishinaModel::~G4KleinNishinaModel() = << 74 G4KleinNishinaModel::~G4KleinNishinaModel()
>> 75 {}
85 76
86 //....oooOO0OOooo........oooOO0OOooo........oo 77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
87 78
88 void G4KleinNishinaModel::Initialise(const G4P 79 void G4KleinNishinaModel::Initialise(const G4ParticleDefinition* p,
89 const G4D << 80 const G4DataVector& cuts)
90 { 81 {
91 fAtomDeexcitation = G4LossTableManager::Inst 82 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
92 if(IsMaster()) { InitialiseElementSelectors( << 83 InitialiseElementSelectors(p, cuts);
93 if(nullptr == fParticleChange) { <<
94 fParticleChange = GetParticleChangeForGamm <<
95 } <<
96 } <<
97 <<
98 //....oooOO0OOooo........oooOO0OOooo........oo <<
99 84
100 void G4KleinNishinaModel::InitialiseLocal(cons << 85 if (isInitialized) { return; }
101 G4VE << 86 fParticleChange = GetParticleChangeForGamma();
102 { << 87 isInitialized = true;
103 SetElementSelectors(masterModel->GetElementS <<
104 } 88 }
105 89
106 //....oooOO0OOooo........oooOO0OOooo........oo 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107 91
108 G4double 92 G4double
109 G4KleinNishinaModel::ComputeCrossSectionPerAto 93 G4KleinNishinaModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
110 << 94 G4double GammaEnergy,
111 << 95 G4double Z, G4double,
112 << 96 G4double, G4double)
113 { 97 {
114 G4double xSection = 0.0 ; << 98 G4double CrossSection = 0.0 ;
115 if (gammaEnergy <= LowEnergyLimit()) { retur << 99 if ( Z < 0.9999 || GammaEnergy < 0.1*keV) { return CrossSection; }
116 100
117 static const G4double a = 20.0 , b = 230.0 , 101 static const G4double a = 20.0 , b = 230.0 , c = 440.0;
118 <<
119 static const G4double <<
120 d1= 2.7965e-1*CLHEP::barn, d2=-1.8300e-1*CLH <<
121 d3= 6.7527 *CLHEP::barn, d4=-1.9798e+1*CLH <<
122 e1= 1.9756e-5*CLHEP::barn, e2=-1.0205e-2*CLH <<
123 e3=-7.3913e-2*CLHEP::barn, e4= 2.7079e-2*CLH <<
124 f1=-3.9178e-7*CLHEP::barn, f2= 6.8241e-5*CLH <<
125 f3= 6.0480e-5*CLHEP::barn, f4= 3.0274e-4*CLH <<
126 102
>> 103 static const G4double
>> 104 d1= 2.7965e-1*barn, d2=-1.8300e-1*barn, d3= 6.7527 *barn, d4=-1.9798e+1*barn,
>> 105 e1= 1.9756e-5*barn, e2=-1.0205e-2*barn, e3=-7.3913e-2*barn, e4= 2.7079e-2*barn,
>> 106 f1=-3.9178e-7*barn, f2= 6.8241e-5*barn, f3= 6.0480e-5*barn, f4= 3.0274e-4*barn;
>> 107
127 G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z = 108 G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z = Z*(d2 + e2*Z + f2*Z*Z),
128 p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z = 109 p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z = Z*(d4 + e4*Z + f4*Z*Z);
129 110
130 G4double T0 = 15.0*keV; 111 G4double T0 = 15.0*keV;
131 if (Z < 1.5) { T0 = 40.0*keV; } 112 if (Z < 1.5) { T0 = 40.0*keV; }
132 113
133 G4double X = max(gammaEnergy, T0) / electr << 114 G4double X = max(GammaEnergy, T0) / electron_mass_c2;
134 xSection = p1Z*G4Log(1.+2.*X)/X << 115 CrossSection = p1Z*std::log(1.+2.*X)/X
135 + (p2Z + p3Z*X + p4Z*X*X)/(1. + 116 + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
136 << 117
137 // modification for low energy. (special ca 118 // modification for low energy. (special case for Hydrogen)
138 static const G4double dT0 = keV; << 119 if (GammaEnergy < T0) {
139 if (gammaEnergy < T0) { << 120 G4double dT0 = keV;
140 X = (T0+dT0) / electron_mass_c2 ; 121 X = (T0+dT0) / electron_mass_c2 ;
141 G4double sigma = p1Z*G4Log(1.+2*X)/X << 122 G4double sigma = p1Z*log(1.+2*X)/X
142 + (p2Z + p3Z*X + p4Z*X*X)/ 123 + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
143 G4double c1 = -T0*(sigma-xSection)/(xSec << 124 G4double c1 = -T0*(sigma-CrossSection)/(CrossSection*dT0);
144 G4double c2 = 0.150; 125 G4double c2 = 0.150;
145 if (Z > 1.5) { c2 = 0.375-0.0556*G4Log(Z); << 126 if (Z > 1.5) { c2 = 0.375-0.0556*log(Z); }
146 G4double y = G4Log(gammaEnergy/T0); << 127 G4double y = log(GammaEnergy/T0);
147 xSection *= G4Exp(-y*(c1+c2*y)); << 128 CrossSection *= exp(-y*(c1+c2*y));
148 } 129 }
149 <<
150 if(xSection < 0.0) { xSection = 0.0; } <<
151 // G4cout << "e= " << GammaEnergy << " Z= " 130 // G4cout << "e= " << GammaEnergy << " Z= " << Z
152 // << " cross= " << xSection << G4endl; << 131 // << " cross= " << CrossSection << G4endl;
153 return xSection; << 132 return CrossSection;
154 } 133 }
155 134
156 //....oooOO0OOooo........oooOO0OOooo........oo 135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
157 136
158 void G4KleinNishinaModel::SampleSecondaries( 137 void G4KleinNishinaModel::SampleSecondaries(
159 std::vector<G4Dyn << 138 std::vector<G4DynamicParticle*>* fvect,
160 const G4MaterialC << 139 const G4MaterialCutsCouple* couple,
161 const G4DynamicPa << 140 const G4DynamicParticle* aDynamicGamma,
162 G4double, << 141 G4double,
163 G4double) << 142 G4double)
164 { 143 {
165 // primary gamma <<
166 G4double energy = aDynamicGamma->GetKineticE 144 G4double energy = aDynamicGamma->GetKineticEnergy();
167 <<
168 // do nothing below the threshold <<
169 if(energy <= LowEnergyLimit()) { return; } <<
170 <<
171 G4ThreeVector direction = aDynamicGamma->Get 145 G4ThreeVector direction = aDynamicGamma->GetMomentumDirection();
172 146
173 // select atom 147 // select atom
174 const G4Element* elm = SelectRandomAtom(coup 148 const G4Element* elm = SelectRandomAtom(couple, theGamma, energy);
175 149
176 // select shell first 150 // select shell first
>> 151 G4int Z = (G4int)elm->GetZ();
177 G4int nShells = elm->GetNbOfAtomicShells(); 152 G4int nShells = elm->GetNbOfAtomicShells();
178 if(nShells > (G4int)fProbabilities.size()) { 153 if(nShells > (G4int)fProbabilities.size()) { fProbabilities.resize(nShells); }
179 G4double totprob = 0.0; 154 G4double totprob = 0.0;
180 G4int i; << 155 G4int i = 0;
181 for(i=0; i<nShells; ++i) { << 156 for(; i<nShells; ++i) {
182 //G4double bindingEnergy = elm->GetAtomicS << 157 G4double prob = 0.0;
183 totprob += elm->GetNbOfShellElectrons(i); << 158 if(energy > elm->GetAtomicShell(i)) {
184 //totprob += elm->GetNbOfShellElectrons(i) << 159 prob = (G4double)elm->GetNbOfShellElectrons(i);
>> 160 }
>> 161 totprob += prob;
185 fProbabilities[i] = totprob; 162 fProbabilities[i] = totprob;
186 } 163 }
>> 164 if(totprob == 0.0) { return; }
187 165
188 // Loop on sampling << 166 G4LorentzVector lv1, lv2, lv3;
189 static const G4int nlooplim = 1000; << 167 G4LorentzVector lv0(energy*direction.x(),energy*direction.y(),
190 G4int nloop = 0; << 168 energy*direction.z(),energy);
191 << 169 G4double eKinEnergy = 0.0;
192 G4double bindingEnergy, ePotEnergy, eKinEner << 170 G4double gamEnergy1 = 0.0;
193 G4double gamEnergy0, gamEnergy1; <<
194 <<
195 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra <<
196 G4double rndm[4]; <<
197 171
>> 172 // Loop on sampling
>> 173 G4double bindingEnergy;
198 do { 174 do {
199 ++nloop; << 175 G4double xprob = totprob*G4UniformRand();
200 176
201 // 4 random numbers to select e- << 177 for(i=0; i<nShells; ++i) { if(xprob <= fProbabilities[i]) {break;} }
202 rndmEngineMod->flatArray(4, rndm); << 178 if( i == nShells ) { return; }
203 G4double xprob = totprob*rndm[0]; <<
204 <<
205 // select shell <<
206 for(i=0; i<nShells; ++i) { if(xprob <= fPr <<
207 179
208 bindingEnergy = elm->GetAtomicShell(i); << 180 bindingEnergy = elm->GetAtomicShell(i);
209 lv1.set(0.0,0.0,energy,energy); << 181 G4double tkin = bindingEnergy*0.5;
210 /* << 182 G4double eEnergy = tkin + electron_mass_c2;
211 G4cout << "nShells= " << nShells << " i= " << 183 G4double eTotMomentum = sqrt(tkin*(tkin + electron_mass_c2*2));
212 << " Egamma= " << energy << " Ebind= " << 184 G4ThreeVector eDir = G4RandomDirection();
213 << G4endl; << 185 lv1 = lv0;
214 */ << 186 lv2.set(eTotMomentum*eDir.x(),eTotMomentum*eDir.y(),
215 // for rest frame of the electron << 187 eTotMomentum*eDir.z(),eEnergy);
216 G4double x = -G4Log(rndm[1]); << 188 G4ThreeVector bst = lv2.boostVector();
217 eKinEnergy = bindingEnergy*x; <<
218 ePotEnergy = bindingEnergy*(1.0 + x); <<
219 <<
220 // for rest frame of the electron <<
221 G4double eTotMomentum = sqrt(eKinEnergy*(e <<
222 G4double phi = rndm[2]*twopi; <<
223 G4double costet = 2*rndm[3] - 1; <<
224 G4double sintet = sqrt((1 - costet)*(1 + c <<
225 lv2.set(eTotMomentum*sintet*cos(phi),eTotM <<
226 eTotMomentum*costet,eKinEnergy + e <<
227 bst = lv2.boostVector(); <<
228 lv1.boost(-bst); 189 lv1.boost(-bst);
229 190
230 gamEnergy0 = lv1.e(); << 191 // In the rest frame of an electron
231 << 192 // The scattered gamma energy is sampled according to Klein - Nishina formula.
232 // In the rest frame of the electron <<
233 // The scattered gamma energy is sampled a <<
234 // The random number techniques of Butcher 193 // The random number techniques of Butcher & Messel are used
235 // (Nuc Phys 20(1960),15). << 194 // (Nuc Phys 20(1960),15).
236 G4double E0_m = gamEnergy0/electron_mass_c << 195
>> 196 G4double gamEnergy0 = lv1.e();
>> 197 G4double E0_m = gamEnergy0 / electron_mass_c2 ;
>> 198
>> 199 G4ThreeVector gamDirection0 = (lv1.vect()).unit();
237 200
238 //G4cout << "Nloop= "<< nloop << " Ecm(keV <<
239 // 201 //
240 // sample the energy rate of the scattered 202 // sample the energy rate of the scattered gamma
241 // 203 //
242 204
243 G4double epsilon, epsilonsq, onecost, sint 205 G4double epsilon, epsilonsq, onecost, sint2, greject ;
244 206
245 G4double eps0 = 1./(1 + 2*E0_m); << 207 G4double epsilon0 = 1./(1. + 2.*E0_m);
246 G4double epsilon0sq = eps0*eps0; << 208 G4double epsilon0sq = epsilon0*epsilon0;
247 G4double alpha1 = - G4Log(eps0); << 209 G4double alpha1 = - log(epsilon0);
248 G4double alpha2 = alpha1 + 0.5*(1 - ep << 210 G4double alpha2 = 0.5*(1.- epsilon0sq);
249 211
250 do { 212 do {
251 ++nloop; << 213 if ( alpha1/(alpha1+alpha2) > G4UniformRand() ) {
252 // false interaction if too many iterati << 214 epsilon = exp(-alpha1*G4UniformRand()); // epsilon0**r
253 if(nloop > nlooplim) { return; } << 215 epsilonsq = epsilon*epsilon;
254 <<
255 // 3 random numbers to sample scattering <<
256 rndmEngineMod->flatArray(3, rndm); <<
257 <<
258 if ( alpha1 > alpha2*rndm[0] ) { <<
259 epsilon = G4Exp(-alpha1*rndm[1]); <<
260 epsilonsq = epsilon*epsilon; <<
261 216
262 } else { 217 } else {
263 epsilonsq = epsilon0sq + (1.- epsilon0 << 218 epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand();
264 epsilon = sqrt(epsilonsq); << 219 epsilon = sqrt(epsilonsq);
265 } << 220 };
266 221
267 onecost = (1.- epsilon)/(epsilon*E0_m); 222 onecost = (1.- epsilon)/(epsilon*E0_m);
268 sint2 = onecost*(2.-onecost); 223 sint2 = onecost*(2.-onecost);
269 greject = 1. - epsilon*sint2/(1.+ epsilo 224 greject = 1. - epsilon*sint2/(1.+ epsilonsq);
270 225
271 // Loop checking, 03-Aug-2015, Vladimir << 226 } while (greject < G4UniformRand());
272 } while (greject < rndm[2]); <<
273 gamEnergy1 = epsilon*gamEnergy0; <<
274 <<
275 // before scattering total 4-momentum in e <<
276 lv2.set(0.0,0.0,0.0,electron_mass_c2); <<
277 lv2 += lv1; <<
278 227
279 // 228 //
280 // scattered gamma angles. ( Z - axis alon 229 // scattered gamma angles. ( Z - axis along the parent gamma)
281 // 230 //
282 if(sint2 < 0.0) { sint2 = 0.0; } <<
283 costet = 1. - onecost; <<
284 sintet = sqrt(sint2); <<
285 phi = twopi * rndmEngineMod->flat(); <<
286 231
287 // e- recoil << 232 G4double cosTeta = 1. - onecost;
>> 233 G4double sinTeta = sqrt (sint2);
>> 234 G4double Phi = twopi * G4UniformRand();
>> 235 G4double dirx = sinTeta*cos(Phi), diry = sinTeta*sin(Phi), dirz = cosTeta;
>> 236
288 // 237 //
289 // in rest frame of the electron << 238 // update G4VParticleChange for the scattered gamma
290 G4ThreeVector gamDir = lv1.vect().unit(); << 239 //
291 G4ThreeVector v = G4ThreeVector(sintet*cos << 240
292 v.rotateUz(gamDir); << 241 G4ThreeVector gamDirection1 ( dirx,diry,dirz );
293 lv1.set(gamEnergy1*v.x(),gamEnergy1*v.y(), << 242 gamDirection1.rotateUz(gamDirection0);
>> 243 gamEnergy1 = epsilon*gamEnergy0;
>> 244
>> 245 // before scattering
>> 246 lv2.set(0.0,0.0,0.0,electron_mass_c2);
>> 247 lv2 += lv1;
>> 248
>> 249 // after scattering
>> 250 lv1.set(gamEnergy1*gamDirection1.x(),gamEnergy1*gamDirection1.y(),
>> 251 gamEnergy1*gamDirection1.z(),gamEnergy1);
294 lv2 -= lv1; 252 lv2 -= lv1;
295 //G4cout<<"Egam(keV)= " << lv1.e()/keV <<
296 // <<" Ee(keV)= " << (lv2.e()-ele <<
297 lv2.boost(bst); 253 lv2.boost(bst);
298 eKinEnergy = lv2.e() - electron_mass_c2 - << 254 lv1.boost(bst);
299 //G4cout << "Nloop= " << nloop << " eKinEn << 255 eKinEnergy = lv2.e() - electron_mass_c2 - bindingEnergy;
300 <<
301 // Loop checking, 03-Aug-2015, Vladimir Iv <<
302 } while ( eKinEnergy < 0.0 ); 256 } while ( eKinEnergy < 0.0 );
303 257
304 // << 258 // gamma kinematics
305 // update G4VParticleChange for the scattere <<
306 // <<
307 <<
308 lv1.boost(bst); <<
309 gamEnergy1 = lv1.e(); 259 gamEnergy1 = lv1.e();
310 if(gamEnergy1 > lowestSecondaryEnergy) { << 260 G4double edep = bindingEnergy;
311 G4ThreeVector gamDirection1 = lv1.vect().u << 261 if(gamEnergy1 > lowestGammaEnergy) {
312 gamDirection1.rotateUz(direction); << 262 fParticleChange->SetProposedKineticEnergy(gamEnergy1);
313 fParticleChange->ProposeMomentumDirection( << 263 fParticleChange->ProposeMomentumDirection((lv1.vect()).unit());
314 } else { 264 } else {
315 fParticleChange->ProposeTrackStatus(fStopA 265 fParticleChange->ProposeTrackStatus(fStopAndKill);
316 gamEnergy1 = 0.0; << 266 fParticleChange->SetProposedKineticEnergy(0.0);
>> 267 edep += gamEnergy1;
317 } 268 }
318 fParticleChange->SetProposedKineticEnergy(ga <<
319 <<
320 // 269 //
321 // kinematic of the scattered electron 270 // kinematic of the scattered electron
322 // 271 //
323 << 272 if(eKinEnergy > DBL_MIN) {
324 if(eKinEnergy > lowestSecondaryEnergy) { << 273 G4ThreeVector eDirection = (lv2.vect()).unit();
325 G4ThreeVector eDirection = lv2.vect().unit << 274 G4DynamicParticle* dp = new G4DynamicParticle(theElectron,eDirection,eKinEnergy);
326 eDirection.rotateUz(direction); <<
327 auto dp = new G4DynamicParticle(theElectro <<
328 fvect->push_back(dp); 275 fvect->push_back(dp);
329 } else { eKinEnergy = 0.0; } << 276 }
330 <<
331 G4double edep = energy - gamEnergy1 - eKinEn <<
332 G4double esec = 0.0; <<
333 <<
334 // sample deexcitation 277 // sample deexcitation
335 // 278 //
336 if(nullptr != fAtomDeexcitation) { << 279 if(fAtomDeexcitation) {
337 G4int index = couple->GetIndex(); 280 G4int index = couple->GetIndex();
338 if(fAtomDeexcitation->CheckDeexcitationAct 281 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
339 G4int Z = elm->GetZasInt(); << 282 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(i);
340 auto as = (G4AtomicShellEnumerator)(i); << 283 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
341 const G4AtomicShell* shell = fAtomDeexci << 284 size_t nbefore = fvect->size();
342 G4int nbefore = (G4int)fvect->size(); <<
343 fAtomDeexcitation->GenerateParticles(fve 285 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
344 G4int nafter = (G4int)fvect->size(); << 286 size_t nafter = fvect->size();
345 //G4cout << "N1= " << nbefore << " N2= << 287 if(nafter > nbefore) {
346 for (G4int j=nbefore; j<nafter; ++j) { << 288 for (size_t i=nbefore; i<nafter; ++i) {
347 G4double e = ((*fvect)[j])->GetKinetic << 289 edep -= ((*fvect)[i])->GetKineticEnergy();
348 if(esec + e > edep) { << 290 }
349 // correct energy in order to have e <<
350 e = edep - esec; <<
351 ((*fvect)[j])->SetKineticEnergy(e); <<
352 esec += e; <<
353 /* <<
354 G4cout << "### G4KleinNishinaModel <<
355 << " Esec(eV)= " << esec/eV <<
356 << " E["<< j << "](eV)= " < <<
357 << " N= " << nafter <<
358 << " Z= " << Z << " shell= <<
359 << " Ebind(keV)= " << bind <<
360 << " Eshell(keV)= " << she <<
361 << G4endl; <<
362 */ <<
363 // delete the rest of secondaries (s <<
364 for (G4int jj=nafter-1; jj>j; --jj) <<
365 delete (*fvect)[jj]; <<
366 fvect->pop_back(); <<
367 } <<
368 break; <<
369 } <<
370 esec += e; <<
371 } 291 }
372 edep -= esec; <<
373 } 292 }
374 } 293 }
375 if(std::abs(energy - gamEnergy1 - eKinEnergy <<
376 G4cout << "### G4KleinNishinaModel dE(eV)= <<
377 << (energy - gamEnergy1 - eKinEnerg <<
378 << " shell= " << i <<
379 << " E(keV)= " << energy/keV <<
380 << " Ebind(keV)= " << bindingEnerg <<
381 << " Eg(keV)= " << gamEnergy1/keV <<
382 << " Ee(keV)= " << eKinEnergy/keV <<
383 << " Esec(keV)= " << esec/keV <<
384 << " Edep(keV)= " << edep/keV <<
385 << G4endl; <<
386 } <<
387 // energy balance 294 // energy balance
388 if(edep > 0.0) { << 295 fParticleChange->ProposeLocalEnergyDeposit(edep);
389 fParticleChange->ProposeLocalEnergyDeposit <<
390 } <<
391 } 296 }
392 297
393 //....oooOO0OOooo........oooOO0OOooo........oo 298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
394 299
395 300