Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4ecpssrBaseLixsModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/lowenergy/src/G4ecpssrBaseLixsModel.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4ecpssrBaseLixsModel.cc (Version 4.0.p2)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26                                                   
 27 #include <iostream>                               
 28 #include "G4ecpssrBaseLixsModel.hh"               
 29 #include "globals.hh"                             
 30 #include "G4PhysicalConstants.hh"                 
 31 #include "G4SystemOfUnits.hh"                     
 32 #include "G4AtomicTransitionManager.hh"           
 33 #include "G4NistManager.hh"                       
 34 #include "G4Proton.hh"                            
 35 #include "G4Alpha.hh"                             
 36 #include "G4LinLogInterpolation.hh"               
 37 #include "G4Exp.hh"                               
 38                                                   
 39 //....oooOO0OOooo........oooOO0OOooo........oo    
 40                                                   
 41 G4ecpssrBaseLixsModel::G4ecpssrBaseLixsModel()    
 42 {                                                 
 43   verboseLevel=0;                                 
 44                                                   
 45   // Storing FLi data needed for 0.2 to 3.0  v    
 46   const char* path = G4FindDataDir("G4LEDATA")    
 47                                                   
 48   if (!path) {                                    
 49     G4Exception("G4ecpssrLCrossSection::G4ecps    
 50     return;                                       
 51   }                                               
 52   std::ostringstream fileName1;                   
 53   std::ostringstream fileName2;                   
 54                                                   
 55   fileName1 << path << "/pixe/uf/FL1.dat";        
 56   fileName2 << path << "/pixe/uf/FL2.dat";        
 57                                                   
 58   // Reading of FL1.dat                           
 59   std::ifstream FL1(fileName1.str().c_str());     
 60   if (!FL1) G4Exception("G4ecpssrLCrossSection    
 61                                                   
 62   dummyVec1.push_back(0.);                        
 63                                                   
 64   while(!FL1.eof())                               
 65     {                                             
 66       double x1;                                  
 67       double y1;                                  
 68                                                   
 69       FL1>>x1>>y1;                                
 70                                                   
 71       //  Mandatory vector initialization         
 72       if (x1 != dummyVec1.back())                 
 73         {                                         
 74           dummyVec1.push_back(x1);                
 75           aVecMap1[x1].push_back(-1.);            
 76         }                                         
 77                                                   
 78       FL1>>FL1Data[x1][y1];                       
 79                                                   
 80       if (y1 != aVecMap1[x1].back()) aVecMap1[    
 81     }                                             
 82                                                   
 83   // Reading of FL2.dat                           
 84                                                   
 85   std::ifstream FL2(fileName2.str().c_str());     
 86   if (!FL2) G4Exception("G4ecpssrLCrossSection    
 87                                                   
 88   dummyVec2.push_back(0.);                        
 89                                                   
 90   while(!FL2.eof())                               
 91     {                                             
 92       double x2;                                  
 93       double y2;                                  
 94                                                   
 95       FL2>>x2>>y2;                                
 96                                                   
 97       //  Mandatory vector initialization         
 98       if (x2 != dummyVec2.back())                 
 99         {                                         
100           dummyVec2.push_back(x2);                
101           aVecMap2[x2].push_back(-1.);            
102         }                                         
103                                                   
104       FL2>>FL2Data[x2][y2];                       
105                                                   
106       if (y2 != aVecMap2[x2].back()) aVecMap2[    
107     }                                             
108 }                                                 
109                                                   
110 //....oooOO0OOooo........oooOO0OOooo........oo    
111                                                   
112 G4ecpssrBaseLixsModel::~G4ecpssrBaseLixsModel(    
113 {}                                                
114                                                   
115 //....oooOO0OOooo........oooOO0OOooo........oo    
116                                                   
117 G4double G4ecpssrBaseLixsModel::ExpIntFunction    
118                                                   
119 {                                                 
120 // this function allows fast evaluation of the    
121   G4int i;                                        
122   G4int ii;                                       
123   G4int nm1;                                      
124   G4double a;                                     
125   G4double b;                                     
126   G4double c;                                     
127   G4double d;                                     
128   G4double del;                                   
129   G4double fact;                                  
130   G4double h;                                     
131   G4double psi;                                   
132   G4double ans = 0;                               
133   static const G4double euler= 0.5772156649;      
134   static const G4int maxit= 100;                  
135   static const G4double fpmin = 1.0e-30;          
136   static const G4double eps = 1.0e-7;             
137   nm1=n-1;                                        
138   if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1    
139   G4cout << "*** WARNING in G4ecpssrBaseLixsMo    
140          << G4endl;                               
141   else {                                          
142     if (n==0) ans=G4Exp(-x)/x;                    
143     else {                                        
144       if (x==0.0) ans=1.0/nm1;                    
145       else {                                      
146   if (x > 1.0) {                                  
147     b=x+n;                                        
148     c=1.0/fpmin;                                  
149     d=1.0/b;                                      
150     h=d;                                          
151     for (i=1;i<=maxit;i++) {                      
152       a=-i*(nm1+i);                               
153       b +=2.0;                                    
154       d=1.0/(a*d+b);                              
155       c=b+a/c;                                    
156       del=c*d;                                    
157       h *=del;                                    
158       if (std::fabs(del-1.0) < eps) {             
159         ans=h*G4Exp(-x);                          
160         return ans;                               
161       }                                           
162     }                                             
163   } else {                                        
164     ans = (nm1!=0 ? 1.0/nm1 : -std::log(x)-eul    
165     fact=1.0;                                     
166     for (i=1;i<=maxit;i++) {                      
167       fact *=-x/i;                                
168       if (i !=nm1) del = -fact/(i-nm1);           
169       else {                                      
170         psi = -euler;                             
171         for (ii=1;ii<=nm1;ii++) psi +=1.0/ii;     
172         del=fact*(-std::log(x)+psi);              
173       }                                           
174       ans += del;                                 
175       if (std::fabs(del) < std::fabs(ans)*eps)    
176     }                                             
177   }                                               
178       }                                           
179     }                                             
180   }                                               
181   return ans;                                     
182 }                                                 
183                                                   
184 //....oooOO0OOooo........oooOO0OOooo........oo    
185                                                   
186 G4double G4ecpssrBaseLixsModel::CalculateL1Cro    
187 {                                                 
188                                                   
189   if (zTarget <=4) return 0.;                     
190                                                   
191   //this L1-CrossSection calculation method is    
192   //and using data tables of O. Benka et al. A    
193                                                   
194   G4NistManager* massManager = G4NistManager::    
195                                                   
196   G4AtomicTransitionManager*  transitionManage    
197                                                   
198   G4double  zIncident = 0;                        
199   G4Proton* aProtone = G4Proton::Proton();        
200   G4Alpha* aAlpha = G4Alpha::Alpha();             
201                                                   
202   if (massIncident == aProtone->GetPDGMass() )    
203     {                                             
204       zIncident = (aProtone->GetPDGCharge())/e    
205     }                                             
206   else                                            
207     {                                             
208       if (massIncident == aAlpha->GetPDGMass()    
209   {                                               
210     zIncident  = (aAlpha->GetPDGCharge())/eplu    
211   }                                               
212       else                                        
213   {                                               
214     G4cout << "*** WARNING in G4ecpssrBaseLixs    
215     G4cout << massIncident << ", " << aAlpha->    
216     return 0;                                     
217   }                                               
218     }                                             
219                                                   
220   G4double l1BindingEnergy = transitionManager    
221   G4double massTarget = (massManager->GetAtomi    
222   G4double systemMass =((massIncident*massTarg    
223   static const G4double zlshell= 4.15;            
224   // *** see Benka, ADANDT 22, p 223              
225   G4double screenedzTarget = zTarget-zlshell;     
226   static const G4double rydbergMeV= 13.6056923    
227                                                   
228   static const G4double nl= 2.;                   
229   // *** see Benka, ADANDT 22, p 220, f3          
230   G4double tetal1 = (l1BindingEnergy*nl*nl)/((    
231   // *** see Benka, ADANDT 22, p 220, f3          
232                                                   
233   if (verboseLevel>0) G4cout << "  tetal1=" <<    
234                                                   
235   G4double reducedEnergy = (energyIncident*ele    
236   // *** also called etaS                         
237   // *** see Benka, ADANDT 22, p 220, f3          
238                                                   
239   static const G4double bohrPow2Barn=(Bohr_rad    
240                                                   
241   G4double sigma0 = 8.*pi*(zIncident*zIncident    
242   // *** see Benka, ADANDT 22, p 220, f2, for     
243   // *** see Basbas, Phys Rev A7, p 1000          
244                                                   
245   G4double velocityl1 = CalculateVelocity(1, z    
246                                                   
247   if (verboseLevel>0) G4cout << "  velocityl1=    
248                                                   
249   static const G4double l1AnalyticalApproximat    
250   G4double x1 =(nl*l1AnalyticalApproximation)/    
251   // *** 1.5 is cK = cL1 (it is 1.25 for L2 &     
252   // *** see Brandt, Phys Rev A20, p 469, f16     
253                                                   
254   if (verboseLevel>0) G4cout << "  x1=" << x1<    
255                                                   
256   G4double electrIonizationEnergyl1=0.;           
257   // *** see Basbas, Phys Rev A17, p1665, f27     
258   // *** see Brandt, Phys Rev A20, p469           
259   // *** see Liu, Comp Phys Comm 97, p325, f A    
260                                                   
261   if ( x1<=0.035)  electrIonizationEnergyl1= 0    
262   else                                            
263     {                                             
264       if ( x1<=3.)                                
265         electrIonizationEnergyl1 =G4Exp(-2.*x1    
266       else                                        
267   {if ( x1<=11.) electrIonizationEnergyl1 =2.*    
268     }                                             
269                                                   
270   G4double hFunctionl1 =(electrIonizationEnerg    
271   // *** see Brandt, Phys Rev A20, p 469, f16     
272                                                   
273   if (verboseLevel>0) G4cout << "  hFunctionl1    
274                                                   
275   G4double gFunctionl1 = (1.+(9.*velocityl1)+(    
276   // *** see Brandt, Phys Rev A20, p 469, f19     
277                                                   
278   if (verboseLevel>0) G4cout << "  gFunctionl1    
279                                                   
280   G4double sigmaPSS_l1 = 1.+(((2.*zIncident)/(    
281   // *** also called dzeta                        
282   // *** also called epsilon                      
283   // *** see Basbas, Phys Rev A17, p1667, f45     
284                                                   
285   if (verboseLevel>0) G4cout << "sigmaPSS_l1 =    
286                                                   
287   const G4double cNaturalUnit= 137.;              
288                                                   
289   G4double yl1Formula=0.4*(screenedzTarget/cNa    
290   // *** also called yS                           
291   // *** see Brandt, Phys Rev A20, p467, f6       
292   // *** see Brandt, Phys Rev A23, p1728          
293                                                   
294   G4double l1relativityCorrection = std::pow((    
295   // *** also called mRS                          
296   // *** see Brandt, Phys Rev A20, p467, f6       
297                                                   
298   //G4double reducedVelocity_l1 = velocityl1*s    
299                                                   
300   G4double L1etaOverTheta2;                       
301                                                   
302   G4double  universalFunction_l1 = 0.;            
303                                                   
304   G4double sigmaPSSR_l1;                          
305                                                   
306   // low velocity formula                         
307   // *****************                            
308   if ( velocityl1 <20.  )                         
309   {                                               
310                                                   
311     L1etaOverTheta2 =(reducedEnergy* l1relativ    
312     // *** 1) RELATIVISTIC CORRECTION ADDED       
313     // *** 2) sigma_PSS_l1 ADDED                  
314     // *** reducedEnergy is etaS, l1relativity    
315     // *** see Phys Rev A20, p468, top            
316                                                   
317     if ( ((tetal1*sigmaPSS_l1) >=0.2) && ((tet    
318       universalFunction_l1 = FunctionFL1((teta    
319                                                   
320     if (verboseLevel>0) G4cout << "at low velo    
321                                                   
322     sigmaPSSR_l1 = (sigma0/(tetal1*sigmaPSS_l1    
323   // *** see Benka, ADANDT 22, p220, f1           
324                                                   
325     if (verboseLevel>0) G4cout << "  at low ve    
326   }                                               
327   else                                            
328   {                                               
329     L1etaOverTheta2 = reducedEnergy/(tetal1*te    
330     // Medium & high velocity                     
331     // *** 1) NO RELATIVISTIC CORRECTION          
332     // *** 2) NO sigma_PSS_l1                     
333     // *** see Benka, ADANDT 22, p223             
334                                                   
335     if ( (tetal1 >=0.2) && (tetal1 <=2.6670) &    
336       universalFunction_l1 = FunctionFL1(tetal    
337                                                   
338     if (verboseLevel>0) G4cout << "at medium a    
339                                                   
340     sigmaPSSR_l1 = (sigma0/tetal1)*universalFu    
341   // *** see Benka, ADANDT 22, p220, f1           
342                                                   
343     if (verboseLevel>0) G4cout << "  sigma PWB    
344   }                                               
345                                                   
346   G4double pssDeltal1 = (4./(systemMass *sigma    
347   // *** also called dzeta*delta                  
348   // *** see Brandt, Phys Rev A23, p1727, f B2    
349                                                   
350   if (verboseLevel>0) G4cout << "  pssDeltal1=    
351                                                   
352   if (pssDeltal1>1) return 0.;                    
353                                                   
354   G4double energyLossl1 = std::pow(1-pssDeltal    
355   // *** also called z                            
356   // *** see Brandt, Phys Rev A23, p1727, afte    
357                                                   
358   if (verboseLevel>0) G4cout << "  energyLossl    
359                                                   
360   G4double coulombDeflectionl1 =                  
361    (8.*pi*zIncident/systemMass)*std::pow(tetal    
362   // *** see Brandt, Phys Rev A20, v2s and f2     
363   // *** with factor n2 compared to Brandt, Ph    
364                                                   
365   G4double cParameterl1 =2.* coulombDeflection    
366   // *** see Brandt, Phys Rev A23, p1727, f B4    
367                                                   
368   G4double coulombDeflectionFunction_l1 = 9.*E    
369                                                   
370   if (verboseLevel>0) G4cout << "  coulombDefl    
371                                                   
372   G4double crossSection_L1 = coulombDeflection    
373                                                   
374   //ECPSSR L1 -subshell cross section is estim    
375   //and reduced by the energy-loss(E),the Coul    
376                                                   
377   if (verboseLevel>0) G4cout << "  crossSectio    
378                                                   
379   if (crossSection_L1 >= 0)                       
380   {                                               
381     return crossSection_L1 * barn;                
382   }                                               
383                                                   
384   else {return 0;}                                
385 }                                                 
386                                                   
387 //....oooOO0OOooo........oooOO0OOooo........oo    
388                                                   
389 G4double G4ecpssrBaseLixsModel::CalculateL2Cro    
390                                                   
391 {                                                 
392   if (zTarget <=13 ) return 0.;                   
393                                                   
394   // this L2-CrossSection calculation method i    
395   // and using data tables of O. Benka et al.     
396                                                   
397   G4NistManager* massManager = G4NistManager::    
398                                                   
399   G4AtomicTransitionManager*  transitionManage    
400                                                   
401   G4double  zIncident = 0;                        
402                                                   
403   G4Proton* aProtone = G4Proton::Proton();        
404   G4Alpha* aAlpha = G4Alpha::Alpha();             
405                                                   
406   if (massIncident == aProtone->GetPDGMass() )    
407    zIncident = (aProtone->GetPDGCharge())/eplu    
408                                                   
409   else                                            
410     {                                             
411       if (massIncident == aAlpha->GetPDGMass()    
412   zIncident  = (aAlpha->GetPDGCharge())/eplus;    
413                                                   
414       else                                        
415   {                                               
416     G4cout << "*** WARNING in G4ecpssrBaseLixs    
417     G4cout << massIncident << ", " << aAlpha->    
418     return 0;                                     
419   }                                               
420     }                                             
421                                                   
422   G4double l2BindingEnergy = transitionManager    
423                                                   
424   G4double massTarget = (massManager->GetAtomi    
425                                                   
426   G4double systemMass =((massIncident*massTarg    
427                                                   
428   const G4double zlshell= 4.15;                   
429                                                   
430   G4double screenedzTarget = zTarget-zlshell;     
431                                                   
432   const G4double rydbergMeV= 13.6056923e-6;       
433                                                   
434   const G4double nl= 2.;                          
435                                                   
436   G4double tetal2 = (l2BindingEnergy*nl*nl)/((    
437                                                   
438   if (verboseLevel>0) G4cout << "  tetal2=" <<    
439                                                   
440   G4double reducedEnergy = (energyIncident*ele    
441                                                   
442   const G4double bohrPow2Barn=(Bohr_radius*Boh    
443                                                   
444   G4double sigma0 = 8.*pi*(zIncident*zIncident    
445                                                   
446   G4double velocityl2 = CalculateVelocity(2,      
447                                                   
448   if (verboseLevel>0) G4cout << "  velocityl2=    
449                                                   
450   const G4double l23AnalyticalApproximation= 1    
451                                                   
452   G4double x2 = (nl*l23AnalyticalApproximation    
453                                                   
454   if (verboseLevel>0) G4cout << "  x2=" << x2<    
455                                                   
456   G4double electrIonizationEnergyl2=0.;           
457                                                   
458   if ( x2<=0.035)  electrIonizationEnergyl2= 0    
459   else                                            
460     {                                             
461       if ( x2<=3.)                                
462         electrIonizationEnergyl2 =G4Exp(-2.*x2    
463       else                                        
464         {if ( x2<=11.) electrIonizationEnergyl    
465     }                                             
466                                                   
467   G4double hFunctionl2 =(electrIonizationEnerg    
468                                                   
469   if (verboseLevel>0) G4cout << "  hFunctionl2    
470                                                   
471   G4double gFunctionl2 = (1.+(10.*velocityl2)+    
472   //takes into account the reduced binding eff    
473                                                   
474   if (verboseLevel>0) G4cout << "  gFunctionl2    
475                                                   
476   G4double sigmaPSS_l2 = 1.+(((2.*zIncident)/(    
477                                                   
478   if (verboseLevel>0) G4cout << "  sigmaPSS_l2    
479                                                   
480   const G4double cNaturalUnit= 137.;              
481                                                   
482   G4double yl2Formula=0.15*(screenedzTarget/cN    
483                                                   
484   G4double l2relativityCorrection = std::pow((    
485                                                   
486   G4double L2etaOverTheta2;                       
487                                                   
488   G4double  universalFunction_l2 = 0.;            
489                                                   
490   G4double sigmaPSSR_l2 ;                         
491                                                   
492   if ( velocityl2 < 20. )                         
493   {                                               
494                                                   
495     L2etaOverTheta2 = (reducedEnergy*l2relativ    
496                                                   
497     if ( (tetal2*sigmaPSS_l2>=0.2) && (tetal2*    
498       universalFunction_l2 = FunctionFL2((teta    
499                                                   
500     sigmaPSSR_l2 = (sigma0/(tetal2*sigmaPSS_l2    
501                                                   
502     if (verboseLevel>0) G4cout << "  sigma PWB    
503   }                                               
504   else                                            
505   {                                               
506                                                   
507     L2etaOverTheta2 = reducedEnergy /(tetal2*t    
508                                                   
509     if ( (tetal2>=0.2) && (tetal2<=2.6670) &&     
510       universalFunction_l2 = FunctionFL2((teta    
511                                                   
512     sigmaPSSR_l2 = (sigma0/tetal2)*universalFu    
513                                                   
514     if (verboseLevel>0) G4cout << "  sigma PWB    
515                                                   
516   }                                               
517                                                   
518   G4double pssDeltal2 = (4./(systemMass*sigmaP    
519                                                   
520   if (pssDeltal2>1) return 0.;                    
521                                                   
522   G4double energyLossl2 = std::pow(1-pssDeltal    
523                                                   
524     if (verboseLevel>0) G4cout << "  energyLos    
525                                                   
526   G4double coulombDeflectionl2                    
527     =(8.*pi*zIncident/systemMass)*std::pow(tet    
528                                                   
529   G4double cParameterl2 = 2.*coulombDeflection    
530                                                   
531   G4double coulombDeflectionFunction_l2 = 11.*    
532   // *** see Brandt, Phys Rev A10, p477, f25      
533                                                   
534   if (verboseLevel>0) G4cout << "  coulombDefl    
535                                                   
536   G4double crossSection_L2 = coulombDeflection    
537   //ECPSSR L2 -subshell cross section is estim    
538   //and reduced by the energy-loss(E),the Coul    
539                                                   
540   if (verboseLevel>0) G4cout << "  crossSectio    
541                                                   
542   if (crossSection_L2 >= 0)                       
543   {                                               
544      return crossSection_L2 * barn;               
545   }                                               
546   else {return 0;}                                
547 }                                                 
548                                                   
549 //....oooOO0OOooo........oooOO0OOooo........oo    
550                                                   
551                                                   
552 G4double G4ecpssrBaseLixsModel::CalculateL3Cro    
553                                                   
554 {                                                 
555   if (zTarget <=13) return 0.;                    
556                                                   
557   //this L3-CrossSection calculation method is    
558   //and using data tables of O. Benka et al. A    
559                                                   
560   G4NistManager* massManager = G4NistManager::    
561                                                   
562   G4AtomicTransitionManager*  transitionManage    
563                                                   
564   G4double  zIncident = 0;                        
565                                                   
566   G4Proton* aProtone = G4Proton::Proton();        
567   G4Alpha* aAlpha = G4Alpha::Alpha();             
568                                                   
569   if (massIncident == aProtone->GetPDGMass() )    
570                                                   
571    zIncident = (aProtone->GetPDGCharge())/eplu    
572                                                   
573   else                                            
574     {                                             
575       if (massIncident == aAlpha->GetPDGMass()    
576                                                   
577     zIncident  = (aAlpha->GetPDGCharge())/eplu    
578                                                   
579       else                                        
580   {                                               
581     G4cout << "*** WARNING in G4ecpssrBaseLixs    
582     G4cout << massIncident << ", " << aAlpha->    
583     return 0;                                     
584   }                                               
585     }                                             
586                                                   
587   G4double l3BindingEnergy = transitionManager    
588                                                   
589   G4double massTarget = (massManager->GetAtomi    
590                                                   
591   G4double systemMass =((massIncident*massTarg    
592                                                   
593   const G4double zlshell= 4.15;                   
594                                                   
595   G4double screenedzTarget = zTarget-zlshell;/    
596                                                   
597   const G4double rydbergMeV= 13.6056923e-6;       
598                                                   
599   const G4double nl= 2.;                          
600                                                   
601   G4double tetal3 = (l3BindingEnergy*nl*nl)/((    
602                                                   
603     if (verboseLevel>0) G4cout << "  tetal3="     
604                                                   
605   G4double reducedEnergy = (energyIncident*ele    
606                                                   
607   const G4double bohrPow2Barn=(Bohr_radius*Boh    
608                                                   
609   G4double sigma0 = 8.*pi*(zIncident*zIncident    
610                                                   
611   G4double velocityl3 = CalculateVelocity(3, z    
612                                                   
613     if (verboseLevel>0) G4cout << "  velocityl    
614                                                   
615   const G4double l23AnalyticalApproximation= 1    
616                                                   
617   G4double x3 = (nl*l23AnalyticalApproximation    
618                                                   
619     if (verboseLevel>0) G4cout << "  x3=" << x    
620                                                   
621   G4double electrIonizationEnergyl3=0.;           
622                                                   
623   if ( x3<=0.035)  electrIonizationEnergyl3= 0    
624     else                                          
625     {                                             
626       if ( x3<=3.) electrIonizationEnergyl3 =G    
627       else                                        
628   {                                               
629     if ( x3<=11.) electrIonizationEnergyl3 =2.    
630     }                                             
631                                                   
632   G4double hFunctionl3 =(electrIonizationEnerg    
633                                                   
634     if (verboseLevel>0) G4cout << "  hFunction    
635                                                   
636   G4double gFunctionl3 = (1.+(10.*velocityl3)+    
637   //takes into account the reduced binding eff    
638                                                   
639     if (verboseLevel>0) G4cout << "  gFunction    
640                                                   
641   G4double sigmaPSS_l3 = 1.+(((2.*zIncident)/(    
642                                                   
643     if (verboseLevel>0) G4cout << "sigmaPSS_l3    
644                                                   
645   const G4double cNaturalUnit= 137.;              
646                                                   
647   G4double yl3Formula=0.15*(screenedzTarget/cN    
648                                                   
649   G4double l3relativityCorrection = std::pow((    
650                                                   
651   G4double L3etaOverTheta2;                       
652                                                   
653   G4double  universalFunction_l3 = 0.;            
654                                                   
655   G4double sigmaPSSR_l3;                          
656                                                   
657   if ( velocityl3 < 20. )                         
658   {                                               
659                                                   
660     L3etaOverTheta2 = (reducedEnergy* l3relati    
661                                                   
662     if ( (tetal3*sigmaPSS_l3>=0.2) && (tetal3*    
663                                                   
664       universalFunction_l3 = 2.*FunctionFL2((t    
665                                                   
666     sigmaPSSR_l3 = (sigma0/(tetal3*sigmaPSS_l3    
667                                                   
668     if (verboseLevel>0) G4cout << "  sigma PWB    
669                                                   
670   }                                               
671                                                   
672   else                                            
673                                                   
674   {                                               
675                                                   
676     L3etaOverTheta2 = reducedEnergy/(tetal3*te    
677                                                   
678     if ( (tetal3>=0.2) && (tetal3<=2.6670) &&     
679                                                   
680       universalFunction_l3 = 2.*FunctionFL2(te    
681                                                   
682     sigmaPSSR_l3 = (sigma0/tetal3)*universalFu    
683                                                   
684       if (verboseLevel>0) G4cout << "  sigma P    
685   }                                               
686                                                   
687   G4double pssDeltal3 = (4./(systemMass*sigmaP    
688                                                   
689     if (verboseLevel>0) G4cout << "  pssDeltal    
690                                                   
691   if (pssDeltal3>1) return 0.;                    
692                                                   
693   G4double energyLossl3 = std::pow(1-pssDeltal    
694                                                   
695   if (verboseLevel>0) G4cout << "  energyLossl    
696                                                   
697   G4double coulombDeflectionl3 =                  
698     (8.*pi*zIncident/systemMass)*std::pow(teta    
699                                                   
700   G4double cParameterl3 = 2.*coulombDeflection    
701                                                   
702   G4double coulombDeflectionFunction_l3 = 11.*    
703   // *** see Brandt, Phys Rev A10, p477, f25      
704                                                   
705     if (verboseLevel>0) G4cout << "  coulombDe    
706                                                   
707   G4double crossSection_L3 =  coulombDeflectio    
708   //ECPSSR L3 -subshell cross section is estim    
709   //and reduced by the energy-loss(E),the Coul    
710                                                   
711     if (verboseLevel>0) G4cout << "  crossSect    
712                                                   
713   if (crossSection_L3 >= 0)                       
714   {                                               
715     return crossSection_L3 * barn;                
716   }                                               
717   else {return 0;}                                
718 }                                                 
719                                                   
720 //....oooOO0OOooo........oooOO0OOooo........oo    
721                                                   
722 G4double G4ecpssrBaseLixsModel::CalculateVeloc    
723                                                   
724 {                                                 
725                                                   
726   G4AtomicTransitionManager*  transitionManage    
727                                                   
728   G4double liBindingEnergy = transitionManager    
729                                                   
730   G4Proton* aProtone = G4Proton::Proton();        
731   G4Alpha* aAlpha = G4Alpha::Alpha();             
732                                                   
733   if (!((massIncident == aProtone->GetPDGMass(    
734     {                                             
735       G4cout << "*** WARNING in G4ecpssrBaseLi    
736       G4cout << massIncident << ", " << aAlpha    
737       return 0;                                   
738     }                                             
739                                                   
740   constexpr G4double zlshell= 4.15;               
741                                                   
742   G4double screenedzTarget = zTarget- zlshell;    
743                                                   
744   constexpr G4double rydbergMeV= 13.6056923e-6    
745                                                   
746   constexpr G4double nl= 2.;                      
747                                                   
748   G4double tetali = (liBindingEnergy*nl*nl)/(s    
749                                                   
750   G4double reducedEnergy = (energyIncident*ele    
751                                                   
752   G4double velocity = 2.*nl*std::pow(reducedEn    
753   // *** see Brandt, Phys Rev A10, p10, f4        
754                                                   
755   return velocity;                                
756 }                                                 
757                                                   
758 //....oooOO0OOooo........oooOO0OOooo........oo    
759                                                   
760 G4double G4ecpssrBaseLixsModel::FunctionFL1(G4    
761 {                                                 
762   G4double sigma = 0.;                            
763   G4double valueT1 = 0;                           
764   G4double valueT2 = 0;                           
765   G4double valueE21 = 0;                          
766   G4double valueE22 = 0;                          
767   G4double valueE12 = 0;                          
768   G4double valueE11 = 0;                          
769   G4double xs11 = 0;                              
770   G4double xs12 = 0;                              
771   G4double xs21 = 0;                              
772   G4double xs22 = 0;                              
773                                                   
774   // PROTECTION TO ALLOW INTERPOLATION AT MINI    
775                                                   
776   if (                                            
777        theta==8.66e-4 ||                          
778        theta==8.66e-3 ||                          
779        theta==8.66e-2 ||                          
780        theta==8.66e-1 ||                          
781        theta==8.66e+00 ||                         
782        theta==8.66e+01                            
783   ) theta=theta-1e-12;                            
784                                                   
785   if (                                            
786        theta==1.e-4 ||                            
787        theta==1.e-3 ||                            
788        theta==1.e-2 ||                            
789        theta==1.e-1 ||                            
790        theta==1.e+00 ||                           
791        theta==1.e+01                              
792   ) theta=theta+1e-12;                            
793                                                   
794   // END PROTECTION                               
795                                                   
796   auto t2 = std::upper_bound(dummyVec1.begin()    
797   auto t1 = t2-1;                                 
798                                                   
799   auto e12 = std::upper_bound(aVecMap1[(*t1)].    
800   auto e11 = e12-1;                               
801                                                   
802   auto e22 = std::upper_bound(aVecMap1[(*t2)].    
803   auto e21 = e22-1;                               
804                                                   
805   valueT1  =*t1;                                  
806   valueT2  =*t2;                                  
807   valueE21 =*e21;                                 
808   valueE22 =*e22;                                 
809   valueE12 =*e12;                                 
810   valueE11 =*e11;                                 
811                                                   
812   xs11 = FL1Data[valueT1][valueE11];              
813   xs12 = FL1Data[valueT1][valueE12];              
814   xs21 = FL1Data[valueT2][valueE21];              
815   xs22 = FL1Data[valueT2][valueE22];              
816                                                   
817   if (verboseLevel>0)                             
818     G4cout                                        
819     << valueT1 << " "                             
820     << valueT2 << " "                             
821     << valueE11 << " "                            
822     << valueE12 << " "                            
823     << valueE21 << " "                            
824     << valueE22 << " "                            
825     << xs11 << " "                                
826     << xs12 << " "                                
827     << xs21 << " "                                
828     << xs22 << " "                                
829     << G4endl;                                    
830                                                   
831   G4double xsProduct = xs11 * xs12 * xs21 * xs    
832                                                   
833   if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0)     
834                                                   
835   if (xsProduct != 0.)                            
836   {                                               
837     sigma = QuadInterpolator(  valueE11, value    
838                  valueE21, valueE22,              
839              xs11, xs12,                          
840              xs21, xs22,                          
841              valueT1, valueT2,                    
842              k, theta );                          
843   }                                               
844                                                   
845   return sigma;                                   
846 }                                                 
847                                                   
848 //....oooOO0OOooo........oooOO0OOooo........oo    
849                                                   
850 G4double G4ecpssrBaseLixsModel::FunctionFL2(G4    
851 {                                                 
852                                                   
853   G4double sigma = 0.;                            
854   G4double valueT1 = 0;                           
855   G4double valueT2 = 0;                           
856   G4double valueE21 = 0;                          
857   G4double valueE22 = 0;                          
858   G4double valueE12 = 0;                          
859   G4double valueE11 = 0;                          
860   G4double xs11 = 0;                              
861   G4double xs12 = 0;                              
862   G4double xs21 = 0;                              
863   G4double xs22 = 0;                              
864                                                   
865   // PROTECTION TO ALLOW INTERPOLATION AT MINI    
866                                                   
867   if (                                            
868        theta==8.66e-4 ||                          
869        theta==8.66e-3 ||                          
870        theta==8.66e-2 ||                          
871        theta==8.66e-1 ||                          
872        theta==8.66e+00 ||                         
873        theta==8.66e+01                            
874   ) theta=theta-1e-12;                            
875                                                   
876   if (                                            
877        theta==1.e-4 ||                            
878        theta==1.e-3 ||                            
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   auto t2 = std::upper_bound(dummyVec2.begin()    
888   auto t1 = t2-1;                                 
889   auto e12 = std::upper_bound(aVecMap2[(*t1)].    
890   auto e11 = e12-1;                               
891   auto e22 = std::upper_bound(aVecMap2[(*t2)].    
892   auto e21 = e22-1;                               
893                                                   
894   valueT1  =*t1;                                  
895   valueT2  =*t2;                                  
896   valueE21 =*e21;                                 
897   valueE22 =*e22;                                 
898   valueE12 =*e12;                                 
899   valueE11 =*e11;                                 
900                                                   
901   xs11 = FL2Data[valueT1][valueE11];              
902   xs12 = FL2Data[valueT1][valueE12];              
903   xs21 = FL2Data[valueT2][valueE21];              
904   xs22 = FL2Data[valueT2][valueE22];              
905                                                   
906   if (verboseLevel>0)                             
907     G4cout                                        
908     << valueT1 << " "                             
909     << valueT2 << " "                             
910     << valueE11 << " "                            
911     << valueE12 << " "                            
912     << valueE21 << " "                            
913     << valueE22 << " "                            
914     << xs11 << " "                                
915     << xs12 << " "                                
916     << xs21 << " "                                
917     << xs22 << " "                                
918     << G4endl;                                    
919                                                   
920   G4double xsProduct = xs11 * xs12 * xs21 * xs    
921                                                   
922   if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0)     
923                                                   
924   if (xsProduct != 0.)                            
925   {                                               
926     sigma = QuadInterpolator(  valueE11, value    
927                  valueE21, valueE22,              
928              xs11, xs12,                          
929              xs21, xs22,                          
930              valueT1, valueT2,                    
931              k, theta );                          
932   }                                               
933                                                   
934   return sigma;                                   
935 }                                                 
936                                                   
937 //....oooOO0OOooo........oooOO0OOooo........oo    
938                                                   
939 G4double G4ecpssrBaseLixsModel::LinLinInterpol    
940                     G4double e2,                  
941                     G4double e,                   
942                     G4double xs1,                 
943                     G4double xs2)                 
944 {                                                 
945   G4double value = xs1 + (xs2 - xs1)*(e - e1)/    
946   return value;                                   
947 }                                                 
948                                                   
949 //....oooOO0OOooo........oooOO0OOooo........oo    
950                                                   
951 G4double G4ecpssrBaseLixsModel::LinLogInterpol    
952                     G4double e2,                  
953                     G4double e,                   
954                     G4double xs1,                 
955                     G4double xs2)                 
956 {                                                 
957   G4double d1 = std::log(xs1);                    
958   G4double d2 = std::log(xs2);                    
959   G4double value = G4Exp(d1 + (d2 - d1)*(e - e    
960   return value;                                   
961 }                                                 
962                                                   
963 //....oooOO0OOooo........oooOO0OOooo........oo    
964                                                   
965 G4double G4ecpssrBaseLixsModel::LogLogInterpol    
966                     G4double e2,                  
967                     G4double e,                   
968                     G4double xs1,                 
969                     G4double xs2)                 
970 {                                                 
971   G4double a = (std::log10(xs2)-std::log10(xs1    
972   G4double b = std::log10(xs2) - a*std::log10(    
973   G4double sigma = a*std::log10(e) + b;           
974   G4double value = (std::pow(10.,sigma));         
975   return value;                                   
976 }                                                 
977                                                   
978 //....oooOO0OOooo........oooOO0OOooo........oo    
979                                                   
980 G4double G4ecpssrBaseLixsModel::QuadInterpolat    
981                    G4double e21, G4double e22,    
982                    G4double xs11, G4double xs1    
983                    G4double xs21, G4double xs2    
984                    G4double t1, G4double t2,      
985                    G4double t, G4double e)        
986 {                                                 
987   // Log-Log                                      
988   G4double interpolatedvalue1 = LogLogInterpol    
989   G4double interpolatedvalue2 = LogLogInterpol    
990   G4double value = LogLogInterpolate(t1, t2, t    
991   return value;                                   
992                                                   
993 }                                                 
994