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

Last change on this file since 1201 was 1197, checked in by garnier, 16 years ago

nvx fichiers dans CVS
  • Property svn:executable set to *
File size: 33.0 KB
RevLine 
[1197]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.