source: trunk/source/processes/electromagnetic/lowenergy/src/G4ecpssrLiCrossSection.cc @ 1252

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

nvx fichiers dans CVS
  • Property svn:executable set to *
File size: 33.0 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: G4ecpssrLiCrossSection.cc,v 1.7 2009/11/11 09:14:53 mantero Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
28//
29// Author: Haifa Ben Abdelouahed
30//         
31//
32// History:
33// -----------
34//  23 Apr 2008   H. Ben Abdelouahed   1st implementation
35//  28 Apr 2008   MGP        Major revision according to a design iteration
36//  29 Apr 2009   ALF Updated Desing for Integration
37//  02 May 2009   ALF + Haifa L1,L2,L3 Extensions
38//  11 Nov 2009   ALF code cleaning for the Dec release
39//
40// -------------------------------------------------------------------
41// Class description:
42// Low Energy Electromagnetic Physics, Cross section, p and alpha ionisation, L shell
43// Further documentation available from http://www.ge.infn.it/geant4/lowE
44// -------------------------------------------------------------------
45
46
47#include "globals.hh"
48#include "G4ecpssrLiCrossSection.hh"
49#include "G4AtomicTransitionManager.hh"
50#include "G4NistManager.hh"
51#include "G4Proton.hh"
52#include "G4Alpha.hh"
53#include <math.h>
54
55G4ecpssrLiCrossSection::G4ecpssrLiCrossSection()
56{ 
57
58
59 // Storing FLi data needed for 0.2 to 3.0  velocities region
60
61    char *path = getenv("G4LEDATA");
62 
63    if (!path)
64    G4Exception("G4ecpssrLCrossSection::CalculateCrossSection: G4LEDDATA environment variable not set");
65
66     std::ostringstream fileName1;
67     std::ostringstream fileName2;
68    fileName1 << path << "/pixe/uf/FL1.dat";
69    fileName2 << path << "/pixe/uf/FL2.dat";
70   
71    std::ifstream FL1(fileName1.str().c_str());
72    std::ifstream FL2(fileName1.str().c_str());
73   
74   
75    if (!FL1) G4Exception("G4ecpssrLCrossSection::CalculateCrossSection: error opening FL1 data file");
76    if (!FL2) G4Exception("G4ecpssrLCrossSection::CalculateCrossSection: error opening FL2 data file");
77 
78 dummyVec.push_back(0.);
79
80    while(!FL1.eof())
81    {
82        double x1;
83        double y1;
84       
85        FL1>>x1>>y1;
86
87        //  Mandatory vector initialization
88        if (x1 != dummyVec.back()) 
89        { 
90          dummyVec.push_back(x1); 
91          aVecMap[x1].push_back(-1.);
92        }
93         
94        FL1>>FL1Data[x1][y1];
95
96        if (y1 != aVecMap[x1].back()) aVecMap[x1].push_back(y1);
97    }
98 while(!FL2.eof())
99    {
100        double x2;
101        double y2;
102       
103        FL2>>x2>>y2;
104
105        //  Mandatory vector initialization
106        if (x2 != dummyVec.back()) 
107        { 
108          dummyVec.push_back(x2); 
109          aVecMap[x2].push_back(-1.);
110        }
111         
112        FL2>>FL2Data[x2][y2];
113
114        if (y2 != aVecMap[x2].back()) aVecMap[x2].push_back(y2);   
115    }
116
117
118}
119
120G4ecpssrLiCrossSection::~G4ecpssrLiCrossSection()
121{ }
122
123//---------------------------------this "ExpIntFunction" function allows fast evaluation of the n order exponential integral function En(x)------
124
125G4double G4ecpssrLiCrossSection::ExpIntFunction(G4int n,G4double x)
126
127{
128  G4int i;
129  G4int ii;
130  G4int nm1;
131  G4double a;
132  G4double b;
133  G4double c;
134  G4double d;
135  G4double del;
136  G4double fact;
137  G4double h;
138  G4double psi;
139  G4double ans = 0;
140  const G4double euler= 0.5772156649;
141  const G4int maxit= 100;
142  const G4double fpmin = 1.0e-30;
143  const G4double eps = 1.0e-7;
144  nm1=n-1;
145  if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1)))
146  G4cout << "bad arguments in ExpIntFunction" << G4endl;
147  else {
148       if (n==0) ans=exp(-x)/x;
149        else {
150           if (x==0.0) ans=1.0/nm1;
151              else {
152                   if (x > 1.0) {
153                                b=x+n;
154                                c=1.0/fpmin;
155                                d=1.0/b;
156                                h=d;
157                                for (i=1;i<=maxit;i++) {
158                                  a=-i*(nm1+i);
159                                  b +=2.0;
160                                  d=1.0/(a*d+b);
161                                  c=b+a/c;
162                                  del=c*d;
163                                  h *=del;
164                                      if (fabs(del-1.0) < eps) {
165                                        ans=h*exp(-x);
166                                        return ans;
167                                      }
168                                }
169                   } else {
170                     ans = (nm1!=0 ? 1.0/nm1 : -log(x)-euler);
171                     fact=1.0;
172                     for (i=1;i<=maxit;i++) {
173                       fact *=-x/i;
174                       if (i !=nm1) del = -fact/(i-nm1);
175                       else {
176                         psi = -euler;
177                         for (ii=1;ii<=nm1;ii++) psi +=1.0/ii;
178                         del=fact*(-log(x)+psi);
179                       }
180                       ans += del;
181                       if (fabs(del) < fabs(ans)*eps) return ans;
182                     }
183                   }
184              }
185        }
186  }
187return ans;
188}
189//-----------------------------------------------------------------------------------------------------------
190
191 
192G4double G4ecpssrLiCrossSection::CalculateL1CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
193 
194 //this L-CrossSection calculation method is done according to W.Brandt and G.Lapicki, Phys.Rev.A23(1981)//                                                     
195
196{
197
198  G4NistManager* massManager = G4NistManager::Instance();   
199
200  G4AtomicTransitionManager*  transitionManager =  G4AtomicTransitionManager::Instance();
201
202
203  G4int zIncident = 0; 
204  G4Proton* aProtone = G4Proton::Proton();
205  G4Alpha* aAlpha = G4Alpha::Alpha();
206
207  if (massIncident == aProtone->GetPDGMass() )
208  {
209
210   zIncident = (G4int)((aProtone->GetPDGCharge())/eplus); 
211
212   //G4cout << "zincident:" << zIncident << G4endl;
213  }
214  else
215    {
216      if (massIncident == aAlpha->GetPDGMass())
217        {
218         
219          zIncident  =(G4int) ((aAlpha->GetPDGCharge())/eplus); 
220         
221          //G4cout << "zincident:" << zIncident << G4endl;
222        }
223      else
224        { 
225          G4cout << "we can treat only Proton or Alpha incident particles " << G4endl;
226          massIncident =0.;
227        }
228    }
229
230
231  G4double l1BindingEnergy = transitionManager->Shell(zTarget,1)->BindingEnergy();
232 
233
234
235
236
237         
238  G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
239 
240  G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;//the mass of the system (projectile, target)
241
242  const G4double zlshell= 4.15;
243
244  G4double screenedzTarget = zTarget-zlshell;         // screenedzTarget is the screened nuclear charge of the target
245
246  const G4double rydbergMeV= 13.6056923e-6;
247
248  const G4double nl= 2.;        // nl is the quantum number of the L shell
249     
250  G4double tetal1 = (l1BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV);  //tetal1 denotes the reduced L1-shell-binding-energy of the electron
251
252
253
254
255
256
257  const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ;       
258 
259  G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*pow(screenedzTarget,-4.);     //sigma0 is the initial cross section of L shell at stable state
260
261  //---------------------------------------------------------------------------------------------------------------------
262
263  G4double velocityl1 = CalculateVelocity(1, zTarget, massIncident, energyIncident);   //is the scaled velocity parameter of the system 
264
265
266
267 
268  //---------------------------------------------------------------------------------------------------------------------
269 
270  const G4double l1AnalyticalApproximation= 1.5; 
271
272
273
274  G4double x1 = nl*l1AnalyticalApproximation/velocityl1;
275
276
277
278
279
280  //-----------------------------------------x of l1 sub shell--------------------------------------
281 
282 G4double electrIonizationEnergyl1;
283                                       
284                                       
285  if ( x1<0.035 && x1>= 0.) 
286    {
287      electrIonizationEnergyl1= 0.75*pi*(log(1./(x1*x1))-1.); 
288    }
289  else 
290    { 
291      if ( x1<3.&& x1>=0.035) 
292        { 
293          electrIonizationEnergyl1 =exp(-2.*x1)/(0.031+(0.213*pow(x1,0.5))+(0.005*x1)-(0.069*pow(x1,3./2.))+(0.324*x1*x1));
294        }
295     
296      else 
297        {
298          if ( x1<=11.&& x1>=3.) {
299            electrIonizationEnergyl1 =2.*exp(-2.*x1)/pow(x1,1.6);
300          }
301          else {
302            electrIonizationEnergyl1 =0.;
303              }
304        }
305    }
306
307
308 
309 
310
311  //-------------------------------------------------------- h and g functions for l1 -------------------------------------------------
312
313
314  G4double hFunctionl1 =(electrIonizationEnergyl1*2.*nl)/(tetal1*pow(velocityl1,3)); //hFunction represents the correction for polarization effet
315   
316 
317  G4double gFunctionl1 = (1.+(9.*velocityl1)+(31.*velocityl1*velocityl1)+(49.*pow(velocityl1,3.))+(162.*pow(velocityl1,4.))+(63.*pow(velocityl1,5.))
318                          +(18.*pow(velocityl1,6.))+(1.97*pow(velocityl1,7.)))/pow(1.+velocityl1,9.); //gFunction represents the correction for binding effet
319
320
321
322
323  //-----------------------------------------------------------------------------------------------------------------------------
324
325  G4double sigmaPSS_l1 = 1.+(((2.*zIncident)/(screenedzTarget*tetal1))*(gFunctionl1-hFunctionl1)); //describes the perturbed stationnairy state of the affected atomic electon
326 
327 
328
329 
330
331
332
333 //----------------------------------------------------------------------------------------------------------------------------
334  //const G4double cNaturalUnit= 1/fine_structure_const;  // it's the speed of light according to Atomic-Unit-System
335  //--------------------------------------------------------------------------------------------------------------
336
337  //G4double yl1Formula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(nl*(velocityl1/sigmaPSS_l1));
338
339 
340
341  //-----------------------------------------------Relativity effect correction L1 -------------------------------------------------------------
342
343  //G4double relativityCorrectionl1 = pow((1.+(1.1*yl1Formula*yl1Formula)),0.5)+yl1Formula;// the relativistic correction parameter
344  //G4double reducedVelocityl1 = velocityl1*pow(relativityCorrectionl1,0.5);  // presents the reduced collision velocity parameter
345
346
347
348
349
350
351
352  //-------------------------------------------------------------------------------------------------------------------
353   //------------------------------------------------------------UNIVERSAL FUNCTION ---------------------------------
354  //------------------------------------------------------------------------------------------------------------
355                                                  // is the reduced universal cross section
356
357  G4double universalFunction_l1 ;
358
359 
360 
361    //-------------------------------------------------------------------------------------------------------------------
362   //------------------------------------------------------------ LIMITS OF ECPSSR MODEL ---------------------------
363  //------------------------------------------------------------------------------------------------------------
364
365
366  G4double L1etaOverTheta2 = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget)/(sigmaPSS_l1*tetal1)/(sigmaPSS_l1*tetal1);
367
368
369     universalFunction_l1 = FunctionFL1((sigmaPSS_l1*tetal1), L1etaOverTheta2);
370       
371
372
373
374     //-----------------------------------------------------------------------------------------------------------------------
375    //---------------------------------------------------------------------------   PSSR Li CROSS SECTION  -------------
376   //-----------------------------------------------------------------------------------------------------------------------
377
378   G4double sigmaPSSR_l1 = (sigma0/(sigmaPSS_l1*tetal1))*universalFunction_l1; //sigmaPSSR is the straight-line L-shell ionization cross section
379
380
381
382   //----------------------------------------------------------------------------------------------------------------------
383 
384  G4double pssDeltal1 = (4./(systemMass*sigmaPSS_l1*tetal1))*(sigmaPSS_l1/velocityl1)*(sigmaPSS_l1/velocityl1);
385
386
387
388  G4double energyLossl1 = pow(1-pssDeltal1,0.5); //energyLoss incorporates the straight-line energy-loss
389
390
391
392
393
394     //-----------------------------------------------------------------------------------------------------------------------
395    //----------------------------------------------------------------------------ENERGY LOSS CORRECTION-------------
396   //-----------------------------------------------------------------------------------------------------------------------
397
398/*
399 G4double energyLossFunction_L1 = (pow(2.,-9)/8.)*((((9.*energyLossl1)-1.)*pow(1.+energyLossl1,9.))+(((9.*energyLossl1)+1.)*pow(1.-energyLossl1,9.)));
400  G4double energyLossFunction_L2 = (pow(2.,-11)/10.)*((((11.*energyLossl2)-1.)*pow(1.+energyLossl2,11.))+(((11.*energyLossl2)+1.)*pow(1.-energyLossl2,11.)));
401  G4double energyLossFunction_L3 = (pow(2.,-11)/10.)*((((11.*energyLossl3)-1.)*pow(1.+energyLossl3,11.))+(((11.*energyLossl3)+1.)*pow(1.-energyLossl3,11.)));
402*/
403
404
405     //-----------------------------------------------------------------------------------------------------------------------
406    //----------------------------------------------------------------------------COULOMB DEFLECTION CORRECTION-------------
407   //-----------------------------------------------------------------------------------------------------------------------
408
409
410  G4double coulombDeflectionl1 = (4.*pi*zIncident/systemMass)*pow(tetal1*sigmaPSS_l1,-2.)*pow(velocityl1/sigmaPSS_l1,-3.)*(zTarget/screenedzTarget); //incorporates Coulomb deflection parameter
411
412
413
414
415  G4double cParameterl1 = 2.*coulombDeflectionl1/(energyLossl1*(energyLossl1+1.));
416
417
418
419  G4double coulombDeflectionFunction_l1 = 9.*ExpIntFunction(10,cParameterl1);                         //this function describes Coulomb-deflection effect
420
421
422
423      //--------------------------------------------------------------------------------------------------------------------------------------------------
424     //-----------------------------------------------------------------------------------------------------------------------
425    //---------------------------------------------------------------------------   ECPSSR Li CROSS SECTION  -------------
426   //-----------------------------------------------------------------------------------------------------------------------
427  //--------------------------------------------------------------------------------------------------------------------------------------------------
428 
429  /*
430   G4double crossSection_L1 = energyLossFunction_L1 * coulombDeflectionFunction_l1 * sigmaPSSR_l1;  //this ECPSSR cross section is estimated at perturbed-stationnairy-state(PSS)
431   G4double crossSection_L2 = energyLossFunction_L2 * coulombDeflectionFunction_l2 * sigmaPSSR_l2;         //and it's reduced by the energy-loss(E),the Coulomb deflection(C),
432   G4double crossSection_L3 = energyLossFunction_L3 * coulombDeflectionFunction_l3 * sigmaPSSR_l3;              //and the relativity(R) effects
433  */
434
435   G4double crossSection_L1 =  coulombDeflectionFunction_l1 * sigmaPSSR_l1; 
436
437 
438 if (crossSection_L1 >= 0) {
439    return crossSection_L1;
440  }
441  else {return 0;}
442}
443
444G4double G4ecpssrLiCrossSection::CalculateL2CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
445 
446 //this L-CrossSection calculation method is done according to W.Brandt and G.Lapicki, Phys.Rev.A23(1981)//                                                     
447
448{
449
450
451  G4NistManager* massManager = G4NistManager::Instance();   
452
453  G4AtomicTransitionManager*  transitionManager =  G4AtomicTransitionManager::Instance();
454
455
456  G4int zIncident = 0; 
457  G4Proton* aProtone = G4Proton::Proton();
458  G4Alpha* aAlpha = G4Alpha::Alpha();
459
460  if (massIncident == aProtone->GetPDGMass() )
461  {
462
463   zIncident =(G4int) ((aProtone->GetPDGCharge())/eplus); 
464
465   //   G4cout << "zincident:" << zIncident << G4endl;
466    }
467  else
468    {
469      if (massIncident == aAlpha->GetPDGMass())
470        {
471         
472          zIncident  = (G4int) ((aAlpha->GetPDGCharge())/eplus); 
473         
474          //      G4cout << "zincident:" << zIncident << G4endl;
475        }
476      else
477        { 
478          G4cout << "we can treat only Proton or Alpha incident particles " << G4endl;
479          massIncident =0.;
480        }
481    }
482
483  G4double l2BindingEnergy = transitionManager->Shell(zTarget,2)->BindingEnergy();
484
485  G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
486  G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;//the mass of the system (projectile, target)
487  const G4double zlshell= 4.15;
488  G4double screenedzTarget = zTarget-zlshell;         // screenedzTarget is the screened nuclear charge of the target
489  const G4double rydbergMeV= 13.6056923e-6;
490  const G4double nl= 2.;        // nl is the quantum number of the L shell
491
492  G4double tetal2 = (l2BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV); 
493
494  const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ;       
495  G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*pow(screenedzTarget,-4.);     //sigma0 is the initial cross section of L shell at stable state
496
497  G4double velocityl2 = CalculateVelocity(2,  zTarget, massIncident, energyIncident);
498
499  const G4double l2AnalyticalApproximation= 1.25;
500
501  G4double x2 = nl*l2AnalyticalApproximation/velocityl2;
502
503  //---------------------------------------------------------------x of l2 sub shell------------------------------------------
504 
505G4double electrIonizationEnergyl2;
506                                       
507                                       
508  if ( x2<0.035 && x2>= 0.) 
509    {
510      electrIonizationEnergyl2= 0.75*pi*(log(1./(x2*x2))-1.); 
511    }
512  else 
513    { 
514      if ( x2<3. && x2 >=0.035) 
515        { 
516          electrIonizationEnergyl2 =exp(-2.*x2)/(0.031+(0.213*pow(x2,0.5))+(0.005*x2)-(0.069*pow(x2,3./2.))+(0.324*x2*x2));
517        }
518     
519      else 
520        {
521         
522          if ( x2<=11.&& x2>=3.) {
523            electrIonizationEnergyl2 =2.*exp(-2.*x2)/pow(x2,1.6);
524          }
525          else {
526            electrIonizationEnergyl2=0.;
527          }
528         
529        }
530    }
531 
532  //----------------------------------------------------------------------------- h and g functions for l2  ----------------------------
533 
534 
535  G4double hFunctionl2 =(electrIonizationEnergyl2*2.*nl)/(tetal2*pow(velocityl2,3)); //hFunction represents the correction for polarization effet
536 
537  G4double gFunctionl2 = (1.+(10.*velocityl2)+(45.*velocityl2*velocityl2)+(102.*pow(velocityl2,3.))+(331.*pow(velocityl2,4.))+(6.7*pow(velocityl2,5.))
538                          +(58.*pow(velocityl2,6.))+(7.8*pow(velocityl2,7.))+ (0.888*pow(velocityl2,8.)) )/pow(1.+velocityl2,10.); //gFunction represents the correction for binding effet
539 
540 
541  G4double sigmaPSS_l2 = 1.+(((2.*zIncident)/(screenedzTarget*tetal2))*(gFunctionl2-hFunctionl2)); 
542 
543  //----------------------------------------------------------------------------------------------------------------------------
544  //const G4double cNaturalUnit= 1/fine_structure_const;  // it's the speed of light according to Atomic-Unit-System
545  //--------------------------------------------------------------------------------------------------------------
546 
547  //G4double yl2Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl2/sigmaPSS_l2);
548 
549  //-----------------------------------------------Relativity effect correction L2 --------------------------------------------------------------
550 
551  //G4double relativityCorrectionl2 = pow((1.+(1.1*yl2Formula*yl2Formula)),0.5)+yl2Formula;// the relativistic correction parameter
552  //  G4double reducedVelocityl2 = velocityl2*pow(relativityCorrectionl2,0.5);  // presents the reduced collision velocity parameter
553 
554  G4double L2etaOverTheta2 = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget)/(sigmaPSS_l2*tetal2)/(sigmaPSS_l2*tetal2);
555 
556 
557  G4double universalFunction_l2 ;
558 
559 
560  universalFunction_l2 = FunctionFL2((sigmaPSS_l2*tetal2), L2etaOverTheta2);
561 
562  G4double sigmaPSSR_l2 = (sigma0/(sigmaPSS_l2*tetal2))*universalFunction_l2;
563  G4double pssDeltal2 = (4./(systemMass*sigmaPSS_l2*tetal2))*(sigmaPSS_l2/velocityl2)*(sigmaPSS_l2/velocityl2);
564  G4double energyLossl2 = pow(1-pssDeltal2,0.5); 
565 
566  G4double coulombDeflectionl2 = (4.*pi*zIncident/systemMass)*pow(tetal2*sigmaPSS_l2,-2.)*pow(velocityl2/sigmaPSS_l2,-3.)*(zTarget/screenedzTarget); //incorporates Coulomb deflection parameter
567  G4double cParameterl2 = 2.*coulombDeflectionl2/(energyLossl2*(energyLossl2+1.));
568  G4double coulombDeflectionFunction_l2 = 11.*ExpIntFunction(12,cParameterl2);
569 
570  G4double crossSection_L2 = coulombDeflectionFunction_l2 * sigmaPSSR_l2 ;         
571 
572 
573  if (crossSection_L2 >= 0) {
574    return crossSection_L2;
575  }
576  else {return 0;}
577 
578 
579}
580
581
582G4double G4ecpssrLiCrossSection::CalculateL3CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
583 
584 //this L-CrossSection calculation method is done according to W.Brandt and G.Lapicki, Phys.Rev.A23(1981)//                                                     
585
586{
587
588
589  G4NistManager* massManager = G4NistManager::Instance();   
590
591  G4AtomicTransitionManager*  transitionManager =  G4AtomicTransitionManager::Instance();
592
593
594  G4int zIncident = 0; 
595  G4Proton* aProtone = G4Proton::Proton();
596  G4Alpha* aAlpha = G4Alpha::Alpha();
597
598  if (massIncident == aProtone->GetPDGMass() )
599  {
600
601   zIncident = (G4int) ((aProtone->GetPDGCharge())/eplus); 
602
603   //   G4cout << "zincident:" << zIncident << G4endl;
604    }
605  else
606    {
607      if (massIncident == aAlpha->GetPDGMass())
608        {
609         
610          zIncident  =(G4int) ((aAlpha->GetPDGCharge())/eplus); 
611         
612          //      G4cout << "zincident:" << zIncident << G4endl;
613        }
614      else
615        { 
616          G4cout << "we can treat only Proton or Alpha incident particles " << G4endl;
617          massIncident =0.;
618        }
619    }
620
621  G4double l3BindingEnergy = transitionManager->Shell(zTarget,3)->BindingEnergy();
622       
623  G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
624 
625  G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;//the mass of the system (projectile, target)
626
627  const G4double zlshell= 4.15;
628
629  G4double screenedzTarget = zTarget-zlshell;         // screenedzTarget is the screened nuclear charge of the target
630
631  const G4double rydbergMeV= 13.6056923e-6;
632
633  const G4double nl= 2.;        // nl is the quantum number of the L shell
634     
635  G4double tetal3 = (l3BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV); 
636
637
638
639  const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ;       
640 
641  G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*pow(screenedzTarget,-4.);     //sigma0 is the initial cross section of L shell at stable state
642
643  //---------------------------------------------------------------------------------------------------------------------
644
645  G4double velocityl3 = CalculateVelocity(3, zTarget, massIncident, energyIncident);
646
647  const G4double l3AnalyticalApproximation= 1.25;
648
649  G4double x3 = nl*l3AnalyticalApproximation/velocityl3;
650
651
652  //--------------------------------------------------------------------- x of l3 sub shell---------------------------------
653
654
655 G4double electrIonizationEnergyl3;
656                                     
657                                       
658  if ( x3<0.035 && x3>=0. ) 
659    {
660      electrIonizationEnergyl3= 0.75*pi*(log(1./(x3*x3))-1.); 
661    }
662  else 
663    { 
664      if ( x3<3. && x3 >= 0.035) 
665        { 
666          electrIonizationEnergyl3 =exp(-2.*x3)/(0.031+(0.213*pow(x3,0.5))+(0.005*x3)-(0.069*pow(x3,3./2.))+(0.324*x3*x3));
667        }
668     
669      else 
670        {
671          if ( x3<=11.&& x3>=3.) {
672            electrIonizationEnergyl3 =2.*exp(-2.*x3)/pow(x3,1.6);
673          }
674          else {
675            electrIonizationEnergyl3=0.;
676          }
677         
678        }
679    }
680
681 //------------------------------------------------------------ h and g function  for l3  ---------------------------------------------
682
683
684  G4double hFunctionl3 =(electrIonizationEnergyl3*2.*nl)/(tetal3*pow(velocityl3,3)); //hFunction represents the correction for polarization effet
685   
686 
687  G4double gFunctionl3 = (1.+(10.*velocityl3)+(45.*velocityl3*velocityl3)+(102.*pow(velocityl3,3.))+(331.*pow(velocityl3,4.))+(6.7*pow(velocityl3,5.))
688                          +(58.*pow(velocityl3,6.))+(7.8*pow(velocityl3,7.))+ (0.888*pow(velocityl3,8.)) )/pow(1.+velocityl3,10.); //gFunction represents the correction for binding effet
689  //-----------------------------------------------------------------------------------------------------------------------------
690
691  G4double sigmaPSS_l3 = 1.+(((2.*zIncident)/(screenedzTarget*tetal3))*(gFunctionl3-hFunctionl3)); 
692
693
694  //----------------------------------------------------------------------------------------------------------------------------
695  //const G4double cNaturalUnit= 1/fine_structure_const;  // it's the speed of light according to Atomic-Unit-System
696  //--------------------------------------------------------------------------------------------------------------
697
698  //G4double yl3Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl3/sigmaPSS_l3);
699
700 //-----------------------------------------------Relativity effect correction L3 --------------------------------------------------------------
701
702 //G4double relativityCorrectionl3 = pow((1.+(1.1*yl3Formula*yl3Formula)),0.5)+yl3Formula;// the relativistic correction parameter
703 //G4double reducedVelocityl3 = velocityl3*pow(relativityCorrectionl3,0.5);  // presents the reduced collision velocity parameter
704
705  G4double universalFunction_l3 ;
706
707  G4double L3etaOverTheta2 = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget)/(sigmaPSS_l3*tetal3)/(sigmaPSS_l3*tetal3);
708 
709  universalFunction_l3 = 2*FunctionFL2((sigmaPSS_l3*tetal3), L3etaOverTheta2);
710 
711 
712  G4double sigmaPSSR_l3 = (sigma0/(sigmaPSS_l3*tetal3))*universalFunction_l3;
713 
714  G4double pssDeltal3 = (4./(systemMass*sigmaPSS_l3*tetal3))*(sigmaPSS_l3/velocityl3)*(sigmaPSS_l3/velocityl3);
715
716  G4double energyLossl3 = pow(1-pssDeltal3,0.5); 
717
718  G4double coulombDeflectionl3 = (4.*pi*zIncident/systemMass)*pow(tetal3*sigmaPSS_l3,-2.)*pow(velocityl3/sigmaPSS_l3,-3.)*(zTarget/screenedzTarget); //incorporates Coulomb deflection parameter
719
720  G4double cParameterl3 = 2.*coulombDeflectionl3/(energyLossl3*(energyLossl3+1.));
721 
722  G4double coulombDeflectionFunction_l3 = 11.*ExpIntFunction(12,cParameterl3);
723 
724  G4double crossSection_L3 = coulombDeflectionFunction_l3 * sigmaPSSR_l3;             
725 
726  if (crossSection_L3 >= 0) {
727    return crossSection_L3;
728  }
729  else {return 0;}
730 
731}
732
733
734
735G4double G4ecpssrLiCrossSection::CalculateVelocity(G4int subShell, G4int zTarget, G4double massIncident,  G4double energyIncident) 
736                                             
737{ 
738
739  G4AtomicTransitionManager*  transitionManager =  G4AtomicTransitionManager::Instance();
740
741  G4double liBindingEnergy = transitionManager->Shell(zTarget,subShell)->BindingEnergy();
742
743 
744  G4Proton* aProtone = G4Proton::Proton();
745  G4Alpha* aAlpha = G4Alpha::Alpha();
746 
747  if (!(massIncident == aProtone->GetPDGMass() || massIncident == aAlpha->GetPDGMass()))
748    {
749      G4cout << "we can treat only Proton or Alpha incident particles " << G4endl;
750      return 0;
751    }
752 
753 
754  const G4double zlshell= 4.15;     
755 
756  G4double screenedzTarget = zTarget- zlshell;                                 
757 
758  const G4double rydbergMeV= 13.6e-6;
759 
760  const G4double nl= 2.;        // nl is the quantum number of the L shell
761 
762  G4double tetali = (liBindingEnergy*nl*nl)/(screenedzTarget*screenedzTarget*rydbergMeV);     
763
764   G4double velocity =(2.*nl/(tetali*screenedzTarget))*pow(((energyIncident*electron_mass_c2)/(massIncident*rydbergMeV)),0.5);
765 
766  return velocity;
767}
768
769
770//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
771
772
773G4double G4ecpssrLiCrossSection::FunctionFL1(G4double k, G4double theta)                                                         
774{
775
776  G4double sigma = 0.;
777  G4double valueT1 = 0;
778  G4double valueT2 = 0;
779  G4double valueE21 = 0;
780  G4double valueE22 = 0;
781  G4double valueE12 = 0;
782  G4double valueE11 = 0;
783  G4double xs11 = 0;   
784  G4double xs12 = 0; 
785  G4double xs21 = 0; 
786  G4double xs22 = 0; 
787
788  // PROTECTION TO ALLOW INTERPOLATION AT MINIMUM AND MAXIMUM EtaK/Theta2 values
789  // (in particular for FK computation at 95 for high velocity formula)
790 
791  if (
792  theta==9.5e-2 ||
793  theta==9.5e-1 ||
794  theta==9.5e+00 ||
795  theta==9.5e+01
796  ) theta=theta-1e-12;
797
798  if (
799  theta==1.e-2 ||
800  theta==1.e-1 ||
801  theta==1.e+00 ||
802  theta==1.e+01
803  ) theta=theta+1e-12;
804
805  // END PROTECTION
806 
807  {
808    std::vector<double>::iterator t2 = std::upper_bound(dummyVec.begin(),dummyVec.end(), k);
809    std::vector<double>::iterator t1 = t2-1;
810 
811    std::vector<double>::iterator e12 = std::upper_bound(aVecMap[(*t1)].begin(),aVecMap[(*t1)].end(), theta);
812    std::vector<double>::iterator e11 = e12-1; 
813         
814    std::vector<double>::iterator e22 = std::upper_bound(aVecMap[(*t2)].begin(),aVecMap[(*t2)].end(), theta);
815    std::vector<double>::iterator e21 = e22-1;
816               
817    valueT1  =*t1;
818    valueT2  =*t2;
819    valueE21 =*e21;
820    valueE22 =*e22;
821    valueE12 =*e12;
822    valueE11 =*e11;
823
824    xs11 = FL1Data[valueT1][valueE11];
825    xs12 = FL1Data[valueT1][valueE12];
826    xs21 = FL1Data[valueT2][valueE21];
827    xs22 = FL1Data[valueT2][valueE22];
828
829}
830     
831  G4double xsProduct = xs11 * xs12 * xs21 * xs22;
832 
833  if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
834     
835  if (xsProduct != 0.)
836  {
837    sigma = QuadInterpolator(  valueE11, valueE12, 
838                               valueE21, valueE22, 
839                               xs11, xs12, 
840                               xs21, xs22, 
841                               valueT1, valueT2, 
842                               k, theta );
843  }
844
845  return sigma;
846}
847
848
849
850//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
851
852
853G4double G4ecpssrLiCrossSection::FunctionFL2(G4double k, G4double theta)                                                         
854{
855
856  G4double sigma = 0.;
857  G4double valueT1 = 0;
858  G4double valueT2 = 0;
859  G4double valueE21 = 0;
860  G4double valueE22 = 0;
861  G4double valueE12 = 0;
862  G4double valueE11 = 0;
863  G4double xs11 = 0;   
864  G4double xs12 = 0; 
865  G4double xs21 = 0; 
866  G4double xs22 = 0; 
867
868  // PROTECTION TO ALLOW INTERPOLATION AT MINIMUM AND MAXIMUM EtaK/Theta2 values
869  // (in particular for FK computation at 95 for high velocity formula)
870 
871  if (
872  theta==9.5e-2 ||
873  theta==9.5e-1 ||
874  theta==9.5e+00 ||
875  theta==9.5e+01
876  ) theta=theta-1e-12;
877
878  if (
879  theta==1.e-2 ||
880  theta==1.e-1 ||
881  theta==1.e+00 ||
882  theta==1.e+01
883  ) theta=theta+1e-12;
884
885  // END PROTECTION
886 
887  {
888    std::vector<double>::iterator t2 = std::upper_bound(dummyVec.begin(),dummyVec.end(), k);
889    std::vector<double>::iterator t1 = t2-1;
890 
891    std::vector<double>::iterator e12 = std::upper_bound(aVecMap[(*t1)].begin(),aVecMap[(*t1)].end(), theta);
892    std::vector<double>::iterator e11 = e12-1; 
893         
894    std::vector<double>::iterator e22 = std::upper_bound(aVecMap[(*t2)].begin(),aVecMap[(*t2)].end(), theta);
895    std::vector<double>::iterator e21 = e22-1;
896               
897    valueT1  =*t1;
898    valueT2  =*t2;
899    valueE21 =*e21;
900    valueE22 =*e22;
901    valueE12 =*e12;
902    valueE11 =*e11;
903
904    xs11 = FL2Data[valueT1][valueE11];
905    xs12 = FL2Data[valueT1][valueE12];
906    xs21 = FL2Data[valueT2][valueE21];
907    xs22 = FL2Data[valueT2][valueE22];
908
909
910}
911     
912  G4double xsProduct = xs11 * xs12 * xs21 * xs22;
913 
914  if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
915     
916  if (xsProduct != 0.)
917  {
918    sigma = QuadInterpolator(  valueE11, valueE12, 
919                               valueE21, valueE22, 
920                               xs11, xs12, 
921                               xs21, xs22, 
922                               valueT1, valueT2, 
923                               k, theta );
924  }
925
926  return sigma;
927}
928
929
930
931//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
932
933G4double G4ecpssrLiCrossSection::LinLogInterpolate(G4double e1, 
934                                                        G4double e2, 
935                                                        G4double e, 
936                                                        G4double xs1, 
937                                                        G4double xs2)
938{
939  G4double d1 = std::log(xs1);
940  G4double d2 = std::log(xs2);
941  G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
942  return value;
943}
944
945//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
946
947G4double G4ecpssrLiCrossSection::LogLogInterpolate(G4double e1, 
948                                                        G4double e2, 
949                                                        G4double e, 
950                                                        G4double xs1, 
951                                                        G4double xs2)
952{
953  G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
954  G4double b = std::log10(xs2) - a*std::log10(e2);
955  G4double sigma = a*std::log10(e) + b;
956  G4double value = (std::pow(10.,sigma));
957  return value;
958}
959
960//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
961
962G4double G4ecpssrLiCrossSection::QuadInterpolator(G4double e11, G4double e12, 
963                                                       G4double e21, G4double e22, 
964                                                       G4double xs11, G4double xs12, 
965                                                       G4double xs21, G4double xs22, 
966                                                       G4double t1, G4double t2, 
967                                                       G4double t, G4double e)
968{
969// Log-Log
970/*
971  G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
972  G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
973  G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
974*/
975
976// Lin-Log
977  G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
978  G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
979  G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
980  return value;
981}
982
Note: See TracBrowser for help on using the repository browser.