source: trunk/source/processes/electromagnetic/lowenergy/src/G4ecpssrCrossSection.cc@ 1192

Last change on this file since 1192 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
File size: 11.5 KB
RevLine 
[968]1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//$Id: G4ecpssrCrossSection.cc,v 1.5 2008/12/18 13:01:32 gunter Exp $
[1007]27// GEANT4 tag $Name: geant4-09-02 $
[968]28//
29// Author: Haifa Ben Abdelouahed
30//
31//
32// History:
33// -----------
34// 21 Apr 2008 H. Ben Abdelouahed 1st implementation
35// 21 Apr 2008 MGP Major revision according to a design iteration
36//
37// -------------------------------------------------------------------
38// Class description:
39// Low Energy Electromagnetic Physics, Cross section, p ionisation, K shell
40// Further documentation available from http://www.ge.infn.it/geant4/lowE
41// -------------------------------------------------------------------
42
43
44#include "globals.hh"
45#include "G4ecpssrCrossSection.hh"
46#include "G4AtomicTransitionManager.hh"
47#include "G4NistManager.hh"
48#include "G4Proton.hh"
49#include "G4Alpha.hh"
50#include <math.h>
51
52G4ecpssrCrossSection::G4ecpssrCrossSection()
53{ }
54
55G4ecpssrCrossSection::~G4ecpssrCrossSection()
56{ }
57
58//---------------------------------this "ExpIntFunction" function allows fast evaluation of the n order exponential integral function En(x)------
59
60G4double G4ecpssrCrossSection::ExpIntFunction(G4int n,G4double x)
61
62{
63 G4int i;
64 G4int ii;
65 G4int nm1;
66 G4double a;
67 G4double b;
68 G4double c;
69 G4double d;
70 G4double del;
71 G4double fact;
72 G4double h;
73 G4double psi;
74 G4double ans = 0;
75 const G4double euler= 0.5772156649;
76 const G4int maxit= 100;
77 const G4double fpmin = 1.0e-30;
78 const G4double eps = 1.0e-7;
79 nm1=n-1;
80 if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1)))
81 G4cout << "bad arguments in ExpIntFunction" << G4endl;
82 else {
83 if (n==0) ans=std::exp(-x)/x;
84 else {
85 if (x==0.0) ans=1.0/nm1;
86 else {
87 if (x > 1.0) {
88 b=x+n;
89 c=1.0/fpmin;
90 d=1.0/b;
91 h=d;
92 for (i=1;i<=maxit;i++) {
93 a=-i*(nm1+i);
94 b +=2.0;
95 d=1.0/(a*d+b);
96 c=b+a/c;
97 del=c*d;
98 h *=del;
99 if (std::fabs(del-1.0) < eps) {
100 ans=h*std::exp(-x);
101 return ans;
102 }
103 }
104 } else {
105 ans = (nm1!=0 ? 1.0/nm1 : -std::log(x)-euler);
106 fact=1.0;
107 for (i=1;i<=maxit;i++) {
108 fact *=-x/i;
109 if (i !=nm1) del = -fact/(i-nm1);
110 else {
111 psi = -euler;
112 for (ii=1;ii<=nm1;ii++) psi +=1.0/ii;
113 del=fact*(-std::log(x)+psi);
114 }
115 ans += del;
116 if (std::fabs(del) < std::fabs(ans)*eps) return ans;
117 }
118 }
119 }
120 }
121 }
122return ans;
123}
124//-----------------------------------------------------------------------------------------------------------
125
126
127G4double G4ecpssrCrossSection::CalculateCrossSection(G4int zTarget,G4int zIncident, G4double energyIncident)
128
129 //this K-CrossSection calculation method is done according to W.Brandt and G.Lapicki, Phys.Rev.A23(1981)//
130
131{
132
133 G4NistManager* massManager = G4NistManager::Instance();
134
135 G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
136
137 G4double massIncident;
138
139 if (zIncident == 1)
140 {
141 G4Proton* aProtone = G4Proton::Proton();
142
143 massIncident = aProtone->GetPDGMass();
144 }
145 else
146 {
147 if (zIncident == 2)
148 {
149 G4Alpha* aAlpha = G4Alpha::Alpha();
150
151 massIncident = aAlpha->GetPDGMass();
152 }
153 else
154 {
155 G4cout << "we can treat only Proton or Alpha incident particles " << G4endl;
156 massIncident =0.;
157 }
158 }
159
160 G4double kBindingEnergy = transitionManager->Shell(zTarget,0)->BindingEnergy();
161
162 G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
163
164 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;//the mass of the system (projectile, target)
165
166 const G4double zkshell= 0.3;
167
168 G4double screenedzTarget = zTarget-zkshell; // screenedzTarget is the screened nuclear charge of the target
169
170 const G4double rydbergMeV= 13.6e-6;
171
172 G4double tetaK = kBindingEnergy/((screenedzTarget*screenedzTarget)*rydbergMeV); //tetaK denotes the reduced binding energy of the electron
173
174 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ;
175
176 G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.); //sigma0 is the initial cross section of K shell at stable state
177
178 //---------------------------------------------------------------------------------------------------------------------
179
180 G4double velocity = CalculateVelocity( zTarget, zIncident, energyIncident); //is the scaled velocity parameter of the system
181
182 //---------------------------------------------------------------------------------------------------------------------
183
184 const G4double kAnalyticalApproximation= 1.5;
185
186 G4double x = kAnalyticalApproximation/velocity;
187
188 G4double electrIonizationEnergy;
189
190 if ( x<0.035)
191 {
192 electrIonizationEnergy= 0.75*pi*(std::log(1./(x*x))-1.);
193 }
194 else
195 {
196 if ( x<3.)
197 {
198 electrIonizationEnergy =std::exp(-2.*x)/(0.031+(0.213*std::pow(x,0.5))+(0.005*x)-(0.069*std::pow(x,3./2.))+(0.324*x*x));
199 }
200
201 else
202 {
203 electrIonizationEnergy =2.*std::exp(-2.*x)/std::pow(x,1.6); }
204 }
205
206 G4double hFunction =(electrIonizationEnergy*2.)/(tetaK*std::pow(velocity,3)); //hFunction represents the correction for polarization effet
207
208 G4double gFunction = (1.+(9.*velocity)+(31.*velocity*velocity)+(98.*std::pow(velocity,3.))+(12.*std::pow(velocity,4.))+(25.*std::pow(velocity,5.))
209 +(4.2*std::pow(velocity,6.))+(0.515*std::pow(velocity,7.)))/std::pow(1.+velocity,9.); //gFunction represents the correction for binding effet
210
211 //-----------------------------------------------------------------------------------------------------------------------------
212
213 G4double sigmaPSS = 1.+(((2.*zIncident)/(screenedzTarget*tetaK))*(gFunction-hFunction)); //describes the perturbed stationnairy state of the affected atomic electon
214
215 //----------------------------------------------------------------------------------------------------------------------------
216
217 const G4double cNaturalUnit= 1/fine_structure_const; // it's the speed of light according to Atomic-Unit-System
218
219 G4double ykFormula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocity/sigmaPSS);
220
221 G4double relativityCorrection = std::pow((1.+(1.1*ykFormula*ykFormula)),0.5)+ykFormula;// the relativistic correction parameter
222
223 G4double reducedVelocity = velocity*std::pow(relativityCorrection,0.5); // presents the reduced collision velocity parameter
224
225 G4double universalFunction = (std::pow(2.,9.)/45.)*std::pow(reducedVelocity/sigmaPSS,8.)*std::pow((1.+(1.72*(reducedVelocity/sigmaPSS)*(reducedVelocity/sigmaPSS))),-4.);// is the reduced universal cross section
226
227 //----------------------------------------------------------------------------------------------------------------------
228
229 G4double sigmaPSSR = (sigma0/(sigmaPSS*tetaK))*universalFunction; //sigmaPSSR is the straight-line K-shell ionization cross section
230
231 //-----------------------------------------------------------------------------------------------------------------------
232
233 G4double pssDeltaK = (4./(systemMass*sigmaPSS*tetaK))*(sigmaPSS/velocity)*(sigmaPSS/velocity);
234
235 G4double energyLoss = std::pow(1-pssDeltaK,0.5); //energyLoss incorporates the straight-line energy-loss
236
237 G4double energyLossFunction = (std::pow(2.,-9)/8.)*((((9.*energyLoss)-1.)*std::pow(1.+energyLoss,9.))+(((9.*energyLoss)+1.)*std::pow(1.-energyLoss,9.)));//energy loss function
238
239 //----------------------------------------------------------------------------------------------------------------------------------------------
240
241 G4double coulombDeflection = (4.*pi*zIncident/systemMass)*std::pow(tetaK*sigmaPSS,-2.)*std::pow(velocity/sigmaPSS,-3.)*(zTarget/screenedzTarget); //incorporates Coulomb deflection parameter
242
243 G4double cParameter = 2.*coulombDeflection/(energyLoss*(energyLoss+1.));
244
245
246 G4double coulombDeflectionFunction = 9.*ExpIntFunction(10,cParameter); //this function describes Coulomb-deflection effect
247
248 //--------------------------------------------------------------------------------------------------------------------------------------------------
249
250
251 G4double crossSection = energyLossFunction* coulombDeflectionFunction*sigmaPSSR; //this ECPSSR cross section is estimated at perturbed-stationnairy-state(PSS)
252 //and it's reduced by the energy-loss(E),the Coulomb deflection(C),
253 //and the relativity(R) effects
254
255 //--------------------------------------------------------------------------------------------------------------------------------------------------
256
257 return crossSection;
258}
259
260G4double G4ecpssrCrossSection::CalculateVelocity(G4int zTarget, G4int zIncident, G4double energyIncident)
261
262{
263
264 G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
265
266 G4double kBindingEnergy = (transitionManager->Shell(zTarget,0)->BindingEnergy())/MeV;
267
268 G4double massIncident;
269
270 if (zIncident == 1)
271 {
272 G4Proton* aProtone = G4Proton::Proton();
273
274 massIncident = aProtone->GetPDGMass();
275 }
276 else
277 {
278 if (zIncident == 2)
279 {
280 G4Alpha* aAlpha = G4Alpha::Alpha();
281
282 massIncident = aAlpha->GetPDGMass();
283 }
284 else
285 {
286 G4cout << "we can treat only Proton or Alpha incident particles " << G4endl;
287 massIncident =0.;
288 }
289 }
290
291 const G4double zkshell= 0.3;
292
293 G4double screenedzTarget = zTarget- zkshell;
294
295 const G4double rydbergMeV= 13.6e-6;
296
297G4double tetaK = kBindingEnergy/(screenedzTarget*screenedzTarget*rydbergMeV);
298
299G4double velocity =(2./(tetaK*screenedzTarget))*std::pow(((energyIncident*electron_mass_c2)/(massIncident*rydbergMeV)),0.5);
300
301 return velocity;
302}
303
Note: See TracBrowser for help on using the repository browser.