| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4KleinNishinaCompton
33 //
34 // Author: Vladimir Ivanchenko on base
35 //
36 // Creation date: 15.03.2005
37 //
38 // Modifications:
39 // 18-04-05 Use G4ParticleChangeForGamma (V.Iv
40 // 27-03-06 Remove upper limit of cross sectio
41 //
42 // Class Description:
43 //
44 // -------------------------------------------
45 //
46 //....oooOO0OOooo........oooOO0OOooo........oo
47 //....oooOO0OOooo........oooOO0OOooo........oo
48
49 #include "G4KleinNishinaCompton.hh"
50 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4Electron.hh"
53 #include "G4Gamma.hh"
54 #include "Randomize.hh"
55 #include "G4DataVector.hh"
56 #include "G4ParticleChangeForGamma.hh"
57 #include "G4Log.hh"
58 #include "G4Exp.hh"
59
60 //....oooOO0OOooo........oooOO0OOooo........oo
61
62 using namespace std;
63
64 G4KleinNishinaCompton::G4KleinNishinaCompton(c
65 c
66 : G4VEmModel(nam)
67 {
68 theGamma = G4Gamma::Gamma();
69 theElectron = G4Electron::Electron();
70 lowestSecondaryEnergy = 100.0*eV;
71 fParticleChange = nullptr;
72 }
73
74 //....oooOO0OOooo........oooOO0OOooo........oo
75
76 G4KleinNishinaCompton::~G4KleinNishinaCompton(
77
78 //....oooOO0OOooo........oooOO0OOooo........oo
79
80 void G4KleinNishinaCompton::Initialise(const G
81 const G
82 {
83 if(IsMaster()) { InitialiseElementSelectors(
84 if(nullptr == fParticleChange) {
85 fParticleChange = GetParticleChangeForGamm
86 }
87 }
88
89 //....oooOO0OOooo........oooOO0OOooo........oo
90
91 void G4KleinNishinaCompton::InitialiseLocal(co
92 G4
93 {
94 SetElementSelectors(masterModel->GetElementS
95 }
96
97 //....oooOO0OOooo........oooOO0OOooo........oo
98
99 G4double G4KleinNishinaCompton::ComputeCrossSe
100 const G
101 G
102 G
103 G
104 {
105 G4double xSection = 0.0 ;
106 if (GammaEnergy <= LowEnergyLimit()) { retur
107
108 static const G4double a = 20.0 , b = 230.0 ,
109
110 static const G4double
111 d1= 2.7965e-1*CLHEP::barn, d2=-1.8300e-1*CLH
112 d3= 6.7527 *CLHEP::barn, d4=-1.9798e+1*CLH
113 e1= 1.9756e-5*CLHEP::barn, e2=-1.0205e-2*CLH
114 e3=-7.3913e-2*CLHEP::barn, e4= 2.7079e-2*CLH
115 f1=-3.9178e-7*CLHEP::barn, f2= 6.8241e-5*CLH
116 f3= 6.0480e-5*CLHEP::barn, f4= 3.0274e-4*CLH
117
118 G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z =
119 p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z =
120
121 G4double T0 = 15.0*keV;
122 if (Z < 1.5) { T0 = 40.0*keV; }
123
124 G4double X = max(GammaEnergy, T0) / electr
125 xSection = p1Z*G4Log(1.+2.*X)/X
126 + (p2Z + p3Z*X + p4Z*X*X)/(1. +
127
128 // modification for low energy. (special ca
129 if (GammaEnergy < T0) {
130 static const G4double dT0 = keV;
131 X = (T0+dT0) / electron_mass_c2 ;
132 G4double sigma = p1Z*G4Log(1.+2*X)/X
133 + (p2Z + p3Z*X + p4Z*X*X)/
134 G4double c1 = -T0*(sigma-xSection)/(xSec
135 G4double c2 = 0.150;
136 if (Z > 1.5) { c2 = 0.375-0.0556*G4Log(Z);
137 G4double y = G4Log(GammaEnergy/T0);
138 xSection *= G4Exp(-y*(c1+c2*y));
139 }
140 // G4cout<<"e= "<< GammaEnergy<<" Z= "<<Z<<"
141 return xSection;
142 }
143
144 //....oooOO0OOooo........oooOO0OOooo........oo
145
146 void G4KleinNishinaCompton::SampleSecondaries(
147 std::vector<G4Dyna
148 const G4MaterialCu
149 const G4DynamicPar
150 G4double,
151 G4double)
152 {
153 // The scattered gamma energy is sampled acc
154 // The random number techniques of Butcher &
155 // (Nuc Phys 20(1960),15).
156 // Note : Effects due to binding of atomic e
157
158 G4double gamEnergy0 = aDynamicGamma->GetKine
159
160 // do nothing below the threshold
161 if(gamEnergy0 <= LowEnergyLimit()) { return;
162
163 G4double E0_m = gamEnergy0 / electron_mass_c
164
165 G4ThreeVector gamDirection0 = aDynamicGamma-
166
167 //
168 // sample the energy rate of the scattered g
169 //
170
171 G4double epsilon, epsilonsq, onecost, sint2,
172
173 G4double eps0 = 1./(1. + 2.*E0_m);
174 G4double epsilon0sq = eps0*eps0;
175 G4double alpha1 = - G4Log(eps0);
176 G4double alpha2 = alpha1 + 0.5*(1.- epsi
177
178 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra
179 G4double rndm[3];
180
181 static const G4int nlooplim = 1000;
182 G4int nloop = 0;
183 do {
184 ++nloop;
185 // false interaction if too many iteration
186 if(nloop > nlooplim) { return; }
187
188 // 3 random numbers to sample scattering
189 rndmEngineMod->flatArray(3, rndm);
190
191 if ( alpha1 > alpha2*rndm[0] ) {
192 epsilon = G4Exp(-alpha1*rndm[1]); //
193 epsilonsq = epsilon*epsilon;
194
195 } else {
196 epsilonsq = epsilon0sq + (1.- epsilon0sq
197 epsilon = sqrt(epsilonsq);
198 };
199
200 onecost = (1.- epsilon)/(epsilon*E0_m);
201 sint2 = onecost*(2.-onecost);
202 greject = 1. - epsilon*sint2/(1.+ epsilons
203
204 // Loop checking, 03-Aug-2015, Vladimir Iv
205 } while (greject < rndm[2]);
206
207 //
208 // scattered gamma angles. ( Z - axis along
209 //
210
211 if(sint2 < 0.0) { sint2 = 0.0; }
212 G4double cosTeta = 1. - onecost;
213 G4double sinTeta = sqrt (sint2);
214 G4double Phi = twopi * rndmEngineMod->fl
215
216 //
217 // update G4VParticleChange for the scattere
218 //
219
220 G4ThreeVector gamDirection1(sinTeta*cos(Phi)
221 gamDirection1.rotateUz(gamDirection0);
222 G4double gamEnergy1 = epsilon*gamEnergy0;
223 G4double edep = 0.0;
224 if(gamEnergy1 > lowestSecondaryEnergy) {
225 fParticleChange->ProposeMomentumDirection(
226 fParticleChange->SetProposedKineticEnergy(
227 } else {
228 fParticleChange->ProposeTrackStatus(fStopA
229 fParticleChange->SetProposedKineticEnergy(
230 edep = gamEnergy1;
231 }
232
233 //
234 // kinematic of the scattered electron
235 //
236
237 G4double eKinEnergy = gamEnergy0 - gamEnergy
238
239 if(eKinEnergy > lowestSecondaryEnergy) {
240 G4ThreeVector eDirection = gamEnergy0*gamD
241 eDirection = eDirection.unit();
242
243 // create G4DynamicParticle object for the
244 auto dp = new G4DynamicParticle(theElectro
245 fvect->push_back(dp);
246 } else {
247 edep += eKinEnergy;
248 }
249 // energy balance
250 if(edep > 0.0) {
251 fParticleChange->ProposeLocalEnergyDeposit
252 }
253 }
254
255 //....oooOO0OOooo........oooOO0OOooo........oo
256
257
258