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

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

update to geant4.9.2
File size: 11.5 KB
Line 
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 $
27// GEANT4 tag $Name: geant4-09-02 $
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.