| 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: G4KleinNishinaCompton.cc,v 1.10 2009/05/15 17:12:33 vnivanch Exp $
>> 27 // GEANT4 tag $Name: geant4-09-03 $
26 // 28 //
27 // ------------------------------------------- 29 // -------------------------------------------------------------------
28 // 30 //
29 // GEANT4 Class file 31 // GEANT4 Class file
30 // 32 //
31 // 33 //
32 // File name: G4KleinNishinaCompton 34 // File name: G4KleinNishinaCompton
33 // 35 //
34 // Author: Vladimir Ivanchenko on base 36 // Author: Vladimir Ivanchenko on base of Michel Maire code
35 // 37 //
36 // Creation date: 15.03.2005 38 // Creation date: 15.03.2005
37 // 39 //
38 // Modifications: 40 // Modifications:
39 // 18-04-05 Use G4ParticleChangeForGamma (V.Iv 41 // 18-04-05 Use G4ParticleChangeForGamma (V.Ivantchenko)
40 // 27-03-06 Remove upper limit of cross sectio 42 // 27-03-06 Remove upper limit of cross section (V.Ivantchenko)
41 // 43 //
42 // Class Description: 44 // Class Description:
43 // 45 //
44 // ------------------------------------------- 46 // -------------------------------------------------------------------
45 // 47 //
46 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48 50
49 #include "G4KleinNishinaCompton.hh" 51 #include "G4KleinNishinaCompton.hh"
50 #include "G4PhysicalConstants.hh" <<
51 #include "G4SystemOfUnits.hh" <<
52 #include "G4Electron.hh" 52 #include "G4Electron.hh"
53 #include "G4Gamma.hh" 53 #include "G4Gamma.hh"
54 #include "Randomize.hh" 54 #include "Randomize.hh"
55 #include "G4DataVector.hh" 55 #include "G4DataVector.hh"
56 #include "G4ParticleChangeForGamma.hh" 56 #include "G4ParticleChangeForGamma.hh"
57 #include "G4Log.hh" <<
58 #include "G4Exp.hh" <<
59 57
60 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61 59
62 using namespace std; 60 using namespace std;
63 61
64 G4KleinNishinaCompton::G4KleinNishinaCompton(c 62 G4KleinNishinaCompton::G4KleinNishinaCompton(const G4ParticleDefinition*,
65 c 63 const G4String& nam)
66 : G4VEmModel(nam) 64 : G4VEmModel(nam)
67 { 65 {
68 theGamma = G4Gamma::Gamma(); 66 theGamma = G4Gamma::Gamma();
69 theElectron = G4Electron::Electron(); 67 theElectron = G4Electron::Electron();
70 lowestSecondaryEnergy = 100.0*eV; << 68 lowestGammaEnergy = 1.0*eV;
71 fParticleChange = nullptr; <<
72 } 69 }
73 70
74 //....oooOO0OOooo........oooOO0OOooo........oo 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 72
76 G4KleinNishinaCompton::~G4KleinNishinaCompton( << 73 G4KleinNishinaCompton::~G4KleinNishinaCompton()
>> 74 {}
77 75
78 //....oooOO0OOooo........oooOO0OOooo........oo 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
79 77
80 void G4KleinNishinaCompton::Initialise(const G << 78 void G4KleinNishinaCompton::Initialise(const G4ParticleDefinition*,
81 const G << 79 const G4DataVector&)
82 { 80 {
83 if(IsMaster()) { InitialiseElementSelectors( << 81 fParticleChange = GetParticleChangeForGamma();
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 } 82 }
96 83
97 //....oooOO0OOooo........oooOO0OOooo........oo 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
98 85
99 G4double G4KleinNishinaCompton::ComputeCrossSe 86 G4double G4KleinNishinaCompton::ComputeCrossSectionPerAtom(
100 const G 87 const G4ParticleDefinition*,
101 G 88 G4double GammaEnergy,
102 G 89 G4double Z, G4double,
103 G 90 G4double, G4double)
104 { 91 {
105 G4double xSection = 0.0 ; << 92 G4double CrossSection = 0.0 ;
106 if (GammaEnergy <= LowEnergyLimit()) { retur << 93 if ( Z < 0.9999 ) return CrossSection;
>> 94 if ( GammaEnergy < 0.1*keV ) return CrossSection;
>> 95 // if ( GammaEnergy > (100.*GeV/Z) ) return CrossSection;
107 96
108 static const G4double a = 20.0 , b = 230.0 , 97 static const G4double a = 20.0 , b = 230.0 , c = 440.0;
109 << 98
110 static const G4double 99 static const G4double
111 d1= 2.7965e-1*CLHEP::barn, d2=-1.8300e-1*CLH << 100 d1= 2.7965e-1*barn, d2=-1.8300e-1*barn, d3= 6.7527 *barn, d4=-1.9798e+1*barn,
112 d3= 6.7527 *CLHEP::barn, d4=-1.9798e+1*CLH << 101 e1= 1.9756e-5*barn, e2=-1.0205e-2*barn, e3=-7.3913e-2*barn, e4= 2.7079e-2*barn,
113 e1= 1.9756e-5*CLHEP::barn, e2=-1.0205e-2*CLH << 102 f1=-3.9178e-7*barn, f2= 6.8241e-5*barn, f3= 6.0480e-5*barn, f4= 3.0274e-4*barn;
114 e3=-7.3913e-2*CLHEP::barn, e4= 2.7079e-2*CLH << 103
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 = 104 G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z = Z*(d2 + e2*Z + f2*Z*Z),
119 p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z = 105 p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z = Z*(d4 + e4*Z + f4*Z*Z);
120 106
121 G4double T0 = 15.0*keV; 107 G4double T0 = 15.0*keV;
122 if (Z < 1.5) { T0 = 40.0*keV; } << 108 if (Z < 1.5) T0 = 40.0*keV;
123 109
124 G4double X = max(GammaEnergy, T0) / electr 110 G4double X = max(GammaEnergy, T0) / electron_mass_c2;
125 xSection = p1Z*G4Log(1.+2.*X)/X << 111 CrossSection = p1Z*std::log(1.+2.*X)/X
126 + (p2Z + p3Z*X + p4Z*X*X)/(1. + 112 + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
127 << 113
128 // modification for low energy. (special ca 114 // modification for low energy. (special case for Hydrogen)
129 if (GammaEnergy < T0) { 115 if (GammaEnergy < T0) {
130 static const G4double dT0 = keV; << 116 G4double dT0 = 1.*keV;
131 X = (T0+dT0) / electron_mass_c2 ; 117 X = (T0+dT0) / electron_mass_c2 ;
132 G4double sigma = p1Z*G4Log(1.+2*X)/X << 118 G4double sigma = p1Z*log(1.+2*X)/X
133 + (p2Z + p3Z*X + p4Z*X*X)/ 119 + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
134 G4double c1 = -T0*(sigma-xSection)/(xSec << 120 G4double c1 = -T0*(sigma-CrossSection)/(CrossSection*dT0);
135 G4double c2 = 0.150; 121 G4double c2 = 0.150;
136 if (Z > 1.5) { c2 = 0.375-0.0556*G4Log(Z); << 122 if (Z > 1.5) c2 = 0.375-0.0556*log(Z);
137 G4double y = G4Log(GammaEnergy/T0); << 123 G4double y = log(GammaEnergy/T0);
138 xSection *= G4Exp(-y*(c1+c2*y)); << 124 CrossSection *= exp(-y*(c1+c2*y));
139 } 125 }
140 // G4cout<<"e= "<< GammaEnergy<<" Z= "<<Z<<" << 126 // G4cout << "e= " << GammaEnergy << " Z= " << Z << " cross= " << CrossSection << G4endl;
141 return xSection; << 127 return CrossSection;
142 } 128 }
143 129
144 //....oooOO0OOooo........oooOO0OOooo........oo 130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
145 131
146 void G4KleinNishinaCompton::SampleSecondaries( << 132 void G4KleinNishinaCompton::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
147 std::vector<G4Dyna << 133 const G4MaterialCutsCouple*,
148 const G4MaterialCu << 134 const G4DynamicParticle* aDynamicGamma,
149 const G4DynamicPar << 135 G4double,
150 G4double, << 136 G4double)
151 G4double) <<
152 { 137 {
153 // The scattered gamma energy is sampled acc 138 // The scattered gamma energy is sampled according to Klein - Nishina formula.
154 // The random number techniques of Butcher & 139 // The random number techniques of Butcher & Messel are used
155 // (Nuc Phys 20(1960),15). 140 // (Nuc Phys 20(1960),15).
156 // Note : Effects due to binding of atomic e 141 // Note : Effects due to binding of atomic electrons are negliged.
157 142
158 G4double gamEnergy0 = aDynamicGamma->GetKine 143 G4double gamEnergy0 = aDynamicGamma->GetKineticEnergy();
159 <<
160 // do nothing below the threshold <<
161 if(gamEnergy0 <= LowEnergyLimit()) { return; <<
162 <<
163 G4double E0_m = gamEnergy0 / electron_mass_c 144 G4double E0_m = gamEnergy0 / electron_mass_c2 ;
164 145
165 G4ThreeVector gamDirection0 = aDynamicGamma- 146 G4ThreeVector gamDirection0 = aDynamicGamma->GetMomentumDirection();
166 147
167 // 148 //
168 // sample the energy rate of the scattered g 149 // sample the energy rate of the scattered gamma
169 // 150 //
170 151
171 G4double epsilon, epsilonsq, onecost, sint2, 152 G4double epsilon, epsilonsq, onecost, sint2, greject ;
172 153
173 G4double eps0 = 1./(1. + 2.*E0_m); << 154 G4double epsilon0 = 1./(1. + 2.*E0_m);
174 G4double epsilon0sq = eps0*eps0; << 155 G4double epsilon0sq = epsilon0*epsilon0;
175 G4double alpha1 = - G4Log(eps0); << 156 G4double alpha1 = - log(epsilon0);
176 G4double alpha2 = alpha1 + 0.5*(1.- epsi << 157 G4double alpha2 = 0.5*(1.- epsilon0sq);
177 <<
178 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra <<
179 G4double rndm[3]; <<
180 158
181 static const G4int nlooplim = 1000; <<
182 G4int nloop = 0; <<
183 do { 159 do {
184 ++nloop; << 160 if ( alpha1/(alpha1+alpha2) > G4UniformRand() ) {
185 // false interaction if too many iteration << 161 epsilon = exp(-alpha1*G4UniformRand()); // epsilon0**r
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; 162 epsilonsq = epsilon*epsilon;
194 163
195 } else { 164 } else {
196 epsilonsq = epsilon0sq + (1.- epsilon0sq << 165 epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand();
197 epsilon = sqrt(epsilonsq); 166 epsilon = sqrt(epsilonsq);
198 }; 167 };
199 168
200 onecost = (1.- epsilon)/(epsilon*E0_m); 169 onecost = (1.- epsilon)/(epsilon*E0_m);
201 sint2 = onecost*(2.-onecost); 170 sint2 = onecost*(2.-onecost);
202 greject = 1. - epsilon*sint2/(1.+ epsilons 171 greject = 1. - epsilon*sint2/(1.+ epsilonsq);
203 172
204 // Loop checking, 03-Aug-2015, Vladimir Iv << 173 } while (greject < G4UniformRand());
205 } while (greject < rndm[2]); <<
206 174
207 // 175 //
208 // scattered gamma angles. ( Z - axis along 176 // scattered gamma angles. ( Z - axis along the parent gamma)
209 // 177 //
210 178
211 if(sint2 < 0.0) { sint2 = 0.0; } <<
212 G4double cosTeta = 1. - onecost; 179 G4double cosTeta = 1. - onecost;
213 G4double sinTeta = sqrt (sint2); 180 G4double sinTeta = sqrt (sint2);
214 G4double Phi = twopi * rndmEngineMod->fl << 181 G4double Phi = twopi * G4UniformRand();
>> 182 G4double dirx = sinTeta*cos(Phi), diry = sinTeta*sin(Phi), dirz = cosTeta;
215 183
216 // 184 //
217 // update G4VParticleChange for the scattere 185 // update G4VParticleChange for the scattered gamma
218 // 186 //
219 187
220 G4ThreeVector gamDirection1(sinTeta*cos(Phi) << 188 G4ThreeVector gamDirection1 ( dirx,diry,dirz );
221 gamDirection1.rotateUz(gamDirection0); 189 gamDirection1.rotateUz(gamDirection0);
222 G4double gamEnergy1 = epsilon*gamEnergy0; 190 G4double gamEnergy1 = epsilon*gamEnergy0;
223 G4double edep = 0.0; << 191 fParticleChange->SetProposedKineticEnergy(gamEnergy1);
224 if(gamEnergy1 > lowestSecondaryEnergy) { << 192 if(gamEnergy1 > lowestGammaEnergy) {
225 fParticleChange->ProposeMomentumDirection( 193 fParticleChange->ProposeMomentumDirection(gamDirection1);
226 fParticleChange->SetProposedKineticEnergy( <<
227 } else { 194 } else {
228 fParticleChange->ProposeTrackStatus(fStopA 195 fParticleChange->ProposeTrackStatus(fStopAndKill);
229 fParticleChange->SetProposedKineticEnergy( << 196 gamEnergy1 += fParticleChange->GetLocalEnergyDeposit();
230 edep = gamEnergy1; << 197 fParticleChange->ProposeLocalEnergyDeposit(gamEnergy1);
231 } 198 }
232 199
233 // 200 //
234 // kinematic of the scattered electron 201 // kinematic of the scattered electron
235 // 202 //
236 203
237 G4double eKinEnergy = gamEnergy0 - gamEnergy 204 G4double eKinEnergy = gamEnergy0 - gamEnergy1;
238 205
239 if(eKinEnergy > lowestSecondaryEnergy) { << 206 if(eKinEnergy > DBL_MIN) {
240 G4ThreeVector eDirection = gamEnergy0*gamD 207 G4ThreeVector eDirection = gamEnergy0*gamDirection0 - gamEnergy1*gamDirection1;
241 eDirection = eDirection.unit(); 208 eDirection = eDirection.unit();
242 209
243 // create G4DynamicParticle object for the 210 // create G4DynamicParticle object for the electron.
244 auto dp = new G4DynamicParticle(theElectro << 211 G4DynamicParticle* dp = new G4DynamicParticle(theElectron,eDirection,eKinEnergy);
245 fvect->push_back(dp); 212 fvect->push_back(dp);
246 } else { <<
247 edep += eKinEnergy; <<
248 } <<
249 // energy balance <<
250 if(edep > 0.0) { <<
251 fParticleChange->ProposeLocalEnergyDeposit <<
252 } 213 }
253 } 214 }
254 215
255 //....oooOO0OOooo........oooOO0OOooo........oo 216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
256 217
257 218
258 219