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