G4ecpssrBaseLixsModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 
00027 #include <cmath>
00028 #include <iostream>
00029 
00030 #include "G4ecpssrBaseLixsModel.hh"
00031 
00032 #include "globals.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "G4AtomicTransitionManager.hh"
00036 #include "G4NistManager.hh"
00037 #include "G4Proton.hh"
00038 #include "G4Alpha.hh"
00039 #include "G4LinLogInterpolation.hh"
00040 
00041 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00042 
00043 G4ecpssrBaseLixsModel::G4ecpssrBaseLixsModel()
00044 {
00045     verboseLevel=0;
00046 
00047     // Storing FLi data needed for 0.2 to 3.0  velocities region
00048 
00049     char *path = getenv("G4LEDATA");
00050 
00051     if (!path) {      
00052       G4Exception("G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()","em0006", FatalException ,"G4LEDATA environment variable not set");
00053       return;
00054     }
00055     std::ostringstream fileName1;
00056     std::ostringstream fileName2;
00057     
00058     fileName1 << path << "/pixe/uf/FL1.dat";
00059     fileName2 << path << "/pixe/uf/FL2.dat";
00060 
00061     // Reading of FL1.dat
00062 
00063     std::ifstream FL1(fileName1.str().c_str());
00064     if (!FL1) G4Exception("G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()","em0003",FatalException, "error opening FL1 data file");
00065   
00066     dummyVec1.push_back(0.);
00067 
00068     while(!FL1.eof())
00069     {
00070         double x1;
00071         double y1;
00072 
00073         FL1>>x1>>y1;
00074 
00075         //  Mandatory vector initialization
00076         if (x1 != dummyVec1.back())
00077         {
00078           dummyVec1.push_back(x1);
00079           aVecMap1[x1].push_back(-1.);
00080         }
00081 
00082         FL1>>FL1Data[x1][y1];
00083 
00084         if (y1 != aVecMap1[x1].back()) aVecMap1[x1].push_back(y1);
00085     }
00086 
00087     // Reading of FL2.dat
00088     
00089     std::ifstream FL2(fileName2.str().c_str());
00090     if (!FL2) G4Exception("G4ecpssrLCrossSection::G4ecpssrBaseLixsModel()","em0003", FatalException," error opening FL2 data file");
00091     
00092     dummyVec2.push_back(0.);
00093 
00094     while(!FL2.eof())
00095     {
00096         double x2;
00097         double y2;
00098 
00099         FL2>>x2>>y2;
00100 
00101         //  Mandatory vector initialization
00102         if (x2 != dummyVec2.back())
00103         {
00104           dummyVec2.push_back(x2);
00105           aVecMap2[x2].push_back(-1.);
00106         }
00107 
00108         FL2>>FL2Data[x2][y2];
00109 
00110         if (y2 != aVecMap2[x2].back()) aVecMap2[x2].push_back(y2);
00111     }
00112 
00113 }
00114 
00115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00116 
00117 G4ecpssrBaseLixsModel::~G4ecpssrBaseLixsModel()
00118 {}
00119 
00120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00121 
00122 G4double G4ecpssrBaseLixsModel::ExpIntFunction(G4int n,G4double x)
00123 
00124 {
00125 // this function allows fast evaluation of the n order exponential integral function En(x)
00126 
00127   G4int i;
00128   G4int ii;
00129   G4int nm1;
00130   G4double a;
00131   G4double b;
00132   G4double c;
00133   G4double d;
00134   G4double del;
00135   G4double fact;
00136   G4double h;
00137   G4double psi;
00138   G4double ans = 0;
00139   const G4double euler= 0.5772156649;
00140   const G4int maxit= 100;
00141   const G4double fpmin = 1.0e-30;
00142   const G4double eps = 1.0e-7;
00143   nm1=n-1;
00144   if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1)))
00145   G4cout << "*** WARNING in G4ecpssrBaseLixsModel::ExpIntFunction: bad arguments in ExpIntFunction" 
00146          << G4endl;
00147   else {
00148        if (n==0) ans=std::exp(-x)/x;
00149         else {
00150            if (x==0.0) ans=1.0/nm1;
00151               else {
00152                    if (x > 1.0) {
00153                                 b=x+n;
00154                                 c=1.0/fpmin;
00155                                 d=1.0/b;
00156                                 h=d;
00157                                 for (i=1;i<=maxit;i++) {
00158                                   a=-i*(nm1+i);
00159                                   b +=2.0;
00160                                   d=1.0/(a*d+b);
00161                                   c=b+a/c;
00162                                   del=c*d;
00163                                   h *=del;
00164                                       if (std::fabs(del-1.0) < eps) {
00165                                         ans=h*std::exp(-x);
00166                                         return ans;
00167                                       }
00168                                 }
00169                    } else {
00170                      ans = (nm1!=0 ? 1.0/nm1 : -std::log(x)-euler);
00171                      fact=1.0;
00172                      for (i=1;i<=maxit;i++) {
00173                        fact *=-x/i;
00174                        if (i !=nm1) del = -fact/(i-nm1);
00175                        else {
00176                          psi = -euler;
00177                          for (ii=1;ii<=nm1;ii++) psi +=1.0/ii;
00178                          del=fact*(-std::log(x)+psi);
00179                        }
00180                        ans += del;
00181                        if (std::fabs(del) < std::fabs(ans)*eps) return ans;
00182                      }
00183                    }
00184               }
00185         }
00186   }
00187 return ans;
00188 }
00189 
00190 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00191 
00192 G4double G4ecpssrBaseLixsModel::CalculateL1CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
00193 {
00194 
00195   if (zTarget <=4) return 0.;
00196 
00197  //this L1-CrossSection calculation method is done according to Werner Brandt and Grzegorz Lapicki, Phys.Rev.A20 N2 (1979),
00198  //and using data tables of O. Benka et al. At.Data Nucl.Data Tables Vol.22 No.3 (1978).
00199 
00200   G4NistManager* massManager = G4NistManager::Instance();
00201 
00202   G4AtomicTransitionManager*  transitionManager =  G4AtomicTransitionManager::Instance();
00203 
00204   G4double  zIncident = 0;
00205   G4Proton* aProtone = G4Proton::Proton();
00206   G4Alpha* aAlpha = G4Alpha::Alpha();
00207 
00208   if (massIncident == aProtone->GetPDGMass() )
00209 
00210    zIncident = (aProtone->GetPDGCharge())/eplus;
00211 
00212   else
00213     {
00214       if (massIncident == aAlpha->GetPDGMass())
00215 
00216           zIncident  = (aAlpha->GetPDGCharge())/eplus;
00217 
00218       else
00219         {
00220           G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateL1CrossSection : Proton or Alpha incident particles only. " << G4endl;
00221           G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
00222           return 0;
00223         }
00224     }
00225 
00226   G4double l1BindingEnergy = transitionManager->Shell(zTarget,1)->BindingEnergy(); //Observed binding energy of L1-subshell
00227 
00228   G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
00229 
00230   G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2; //Mass of the system (projectile, target)
00231 
00232   const G4double zlshell= 4.15;
00233   // *** see Benka, ADANDT 22, p 223
00234 
00235   G4double screenedzTarget = zTarget-zlshell; //Effective nuclear charge as seen by electrons in L1-sub shell
00236 
00237   const G4double rydbergMeV= 13.6056923e-6;
00238 
00239   const G4double nl= 2.;
00240   // *** see Benka, ADANDT 22, p 220, f3
00241 
00242   G4double tetal1 = (l1BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV); //Screening parameter
00243   // *** see Benka, ADANDT 22, p 220, f3
00244 
00245     if (verboseLevel>0) G4cout << "  tetal1=" <<  tetal1<< G4endl;
00246 
00247   G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
00248   // *** also called etaS
00249   // *** see Benka, ADANDT 22, p 220, f3
00250 
00251   const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ; //Bohr radius of hydrogen
00252 
00253   G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
00254   // *** see Benka, ADANDT 22, p 220, f2, for protons
00255   // *** see Basbas, Phys Rev A7, p 1000
00256 
00257   G4double velocityl1 = CalculateVelocity(1, zTarget, massIncident, energyIncident); // Scaled velocity
00258 
00259     if (verboseLevel>0) G4cout << "  velocityl1=" << velocityl1<< G4endl;
00260 
00261   const G4double l1AnalyticalApproximation= 1.5;
00262   G4double x1 =(nl*l1AnalyticalApproximation)/velocityl1;
00263   // *** 1.5 is cK = cL1 (it is 1.25 for L2 & L3)
00264   // *** see Brandt, Phys Rev A20, p 469, f16 in expression of h
00265   
00266    if (verboseLevel>0) G4cout << "  x1=" << x1<< G4endl;
00267 
00268   G4double electrIonizationEnergyl1=0.;
00269   // *** see Basbas, Phys Rev A17, p1665, f27
00270   // *** see Brandt, Phys Rev A20, p469
00271   // *** see Liu, Comp Phys Comm 97, p325, f A5
00272 
00273   if ( x1<=0.035)  electrIonizationEnergyl1= 0.75*pi*(std::log(1./(x1*x1))-1.);
00274   else
00275     {
00276       if ( x1<=3.)
00277         electrIonizationEnergyl1 =std::exp(-2.*x1)/(0.031+(0.213*std::pow(x1,0.5))+(0.005*x1)-(0.069*std::pow(x1,3./2.))+(0.324*x1*x1));
00278       else
00279         {if ( x1<=11.) electrIonizationEnergyl1 =2.*std::exp(-2.*x1)/std::pow(x1,1.6);}
00280     }
00281 
00282   G4double hFunctionl1 =(electrIonizationEnergyl1*2.*nl)/(tetal1*std::pow(velocityl1,3)); //takes into account the polarization effect
00283   // *** see Brandt, Phys Rev A20, p 469, f16
00284 
00285     if (verboseLevel>0) G4cout << "  hFunctionl1=" << hFunctionl1<< G4endl;
00286 
00287   G4double gFunctionl1 = (1.+(9.*velocityl1)+(31.*velocityl1*velocityl1)+(49.*std::pow(velocityl1,3.))+(162.*std::pow(velocityl1,4.))+(63.*std::pow(velocityl1,5.))+(18.*std::pow(velocityl1,6.))+(1.97*std::pow(velocityl1,7.)))/std::pow(1.+velocityl1,9.);//takes into account the reduced binding effect
00288   // *** see Brandt, Phys Rev A20, p 469, f19
00289 
00290     if (verboseLevel>0) G4cout << "  gFunctionl1=" << gFunctionl1<< G4endl;
00291 
00292   G4double sigmaPSS_l1 = 1.+(((2.*zIncident)/(screenedzTarget*tetal1))*(gFunctionl1-hFunctionl1)); //Binding-polarization factor
00293   // *** also called dzeta
00294   // *** also called epsilon
00295   // *** see Basbas, Phys Rev A17, p1667, f45
00296 
00297     if (verboseLevel>0) G4cout << "sigmaPSS_l1 =" << sigmaPSS_l1<< G4endl;
00298 
00299   const G4double cNaturalUnit= 137.;
00300   
00301   G4double yl1Formula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(nl*velocityl1/sigmaPSS_l1);
00302   // *** also called yS
00303   // *** see Brandt, Phys Rev A20, p467, f6
00304   // *** see Brandt, Phys Rev A23, p1728
00305       
00306   G4double l1relativityCorrection = std::pow((1.+(1.1*yl1Formula*yl1Formula)),0.5)+yl1Formula; // Relativistic correction parameter
00307   // *** also called mRS
00308   // *** see Brandt, Phys Rev A20, p467, f6
00309   
00310   //G4double reducedVelocity_l1 = velocityl1*std::pow(l1relativityCorrection,0.5); //Reduced velocity parameter
00311 
00312   G4double L1etaOverTheta2;
00313 
00314   G4double  universalFunction_l1 = 0.;
00315 
00316   G4double sigmaPSSR_l1;
00317 
00318   // low velocity formula
00319   // *****************  
00320   if ( velocityl1 <20.  )
00321   {
00322 
00323     L1etaOverTheta2 =(reducedEnergy* l1relativityCorrection)/((tetal1*sigmaPSS_l1)*(tetal1*sigmaPSS_l1));
00324     // *** 1) RELATIVISTIC CORRECTION ADDED
00325     // *** 2) sigma_PSS_l1 ADDED
00326     // *** reducedEnergy is etaS, l1relativityCorrection is mRS
00327     // *** see Phys Rev A20, p468, top
00328     
00329     if ( ((tetal1*sigmaPSS_l1) >=0.2) && ((tetal1*sigmaPSS_l1) <=2.6670) && (L1etaOverTheta2>=0.1e-3) && (L1etaOverTheta2<=0.866e2) )
00330   
00331       universalFunction_l1 = FunctionFL1((tetal1*sigmaPSS_l1),L1etaOverTheta2);
00332 
00333       if (verboseLevel>0) G4cout << "at low velocity range, universalFunction_l1  =" << universalFunction_l1 << G4endl;
00334 
00335     sigmaPSSR_l1 = (sigma0/(tetal1*sigmaPSS_l1))*universalFunction_l1;// Plane-wave Born -Aproximation L1-subshell ionisation Cross Section
00336   // *** see Benka, ADANDT 22, p220, f1
00337 
00338       if (verboseLevel>0) G4cout << "  at low velocity range, sigma PWBA L1 CS  = " << sigmaPSSR_l1<< G4endl;
00339 
00340    }
00341 
00342   else
00343 
00344   {
00345 
00346     L1etaOverTheta2 = reducedEnergy/(tetal1*tetal1);
00347     // Medium & high velocity
00348     // *** 1) NO RELATIVISTIC CORRECTION
00349     // *** 2) NO sigma_PSS_l1
00350     // *** see Benka, ADANDT 22, p223
00351 
00352     if ( (tetal1 >=0.2) && (tetal1 <=2.6670) && (L1etaOverTheta2>=0.1e-3) && (L1etaOverTheta2<=0.866e2) ) 
00353     
00354       universalFunction_l1 = FunctionFL1(tetal1,L1etaOverTheta2);
00355 
00356     if (verboseLevel>0) G4cout << "at medium and high velocity range, universalFunction_l1  =" << universalFunction_l1 << G4endl;
00357 
00358     sigmaPSSR_l1 = (sigma0/tetal1)*universalFunction_l1;// Plane-wave Born -Aproximation L1-subshell ionisation Cross Section
00359   // *** see Benka, ADANDT 22, p220, f1
00360 
00361     if (verboseLevel>0) G4cout << "  sigma PWBA L1 CS at medium and high velocity range = " << sigmaPSSR_l1<< G4endl;    
00362   }
00363  
00364   G4double pssDeltal1 = (4./(systemMass *sigmaPSS_l1*tetal1))*(sigmaPSS_l1/velocityl1)*(sigmaPSS_l1/velocityl1);
00365   // *** also called dzeta*delta
00366   // *** see Brandt, Phys Rev A23, p1727, f B2
00367 
00368     if (verboseLevel>0) G4cout << "  pssDeltal1=" << pssDeltal1<< G4endl;
00369 
00370   if (pssDeltal1>1) return 0.;
00371 
00372   G4double energyLossl1 = std::pow(1-pssDeltal1,0.5);
00373   // *** also called z
00374   // *** see Brandt, Phys Rev A23, p1727, after f B2
00375 
00376     if (verboseLevel>0) G4cout << "  energyLossl1=" << energyLossl1<< G4endl;
00377 
00378   G4double coulombDeflectionl1 = 
00379    (8.*pi*zIncident/systemMass)*std::pow(tetal1*sigmaPSS_l1,-2.)*std::pow(velocityl1/sigmaPSS_l1,-3.)*(zTarget/screenedzTarget);
00380   // *** see Brandt, Phys Rev A20, v2s and f2 and B2 
00381   // *** with factor n2 compared to Brandt, Phys Rev A23, p1727, f B3
00382 
00383   G4double cParameterl1 =2.* coulombDeflectionl1/(energyLossl1*(energyLossl1+1.));
00384   // *** see Brandt, Phys Rev A23, p1727, f B4
00385 
00386   G4double coulombDeflectionFunction_l1 = 9.*ExpIntFunction(10,cParameterl1); //Coulomb-deflection effect correction
00387 
00388     if (verboseLevel>0) G4cout << "  coulombDeflectionFunction_l1 =" << coulombDeflectionFunction_l1 << G4endl;
00389 
00390   G4double crossSection_L1 = coulombDeflectionFunction_l1 * sigmaPSSR_l1;
00391 
00392   //ECPSSR L1 -subshell cross section is estimated at perturbed-stationnairy-state(PSS)
00393   //and reduced by the energy-loss(E),the Coulomb deflection(C),and the relativity(R) effects
00394 
00395     if (verboseLevel>0) G4cout << "  crossSection_L1 =" << crossSection_L1 << G4endl;
00396 
00397   if (crossSection_L1 >= 0) 
00398   {
00399     return crossSection_L1 * barn;
00400   }
00401   
00402   else {return 0;}
00403 }
00404 
00405 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00406 
00407 G4double G4ecpssrBaseLixsModel::CalculateL2CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
00408 
00409 {
00410   if (zTarget <=13 ) return 0.;
00411 
00412   // this L2-CrossSection calculation method is done according to Werner Brandt and Grzegorz Lapicki, Phys.Rev.A20 N2 (1979),
00413   // and using data tables of O. Benka et al. At.Data Nucl.Data Tables Vol.22 No.3 (1978).
00414 
00415   G4NistManager* massManager = G4NistManager::Instance();
00416 
00417   G4AtomicTransitionManager*  transitionManager =  G4AtomicTransitionManager::Instance();
00418 
00419   G4double  zIncident = 0;
00420   
00421   G4Proton* aProtone = G4Proton::Proton();
00422   G4Alpha* aAlpha = G4Alpha::Alpha();
00423 
00424   if (massIncident == aProtone->GetPDGMass() )
00425 
00426    zIncident = (aProtone->GetPDGCharge())/eplus;
00427 
00428   else
00429     {
00430       if (massIncident == aAlpha->GetPDGMass())
00431 
00432           zIncident  = (aAlpha->GetPDGCharge())/eplus;
00433 
00434       else
00435         {
00436           G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateL2CrossSection : Proton or Alpha incident particles only. " << G4endl;
00437           G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
00438           return 0;
00439         }
00440     }
00441 
00442   G4double l2BindingEnergy = transitionManager->Shell(zTarget,2)->BindingEnergy(); //Observed binding energy of L2-subshell
00443 
00444   G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
00445 
00446   G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2; //Mass of the system (projectile, target)
00447 
00448   const G4double zlshell= 4.15;
00449 
00450   G4double screenedzTarget = zTarget-zlshell; //Effective nuclear charge as seen by electrons in L2-subshell
00451 
00452   const G4double rydbergMeV= 13.6056923e-6;
00453 
00454   const G4double nl= 2.;
00455 
00456   G4double tetal2 = (l2BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV); //Screening parameter
00457 
00458     if (verboseLevel>0) G4cout << "  tetal2=" <<  tetal2<< G4endl;
00459 
00460   G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget); 
00461 
00462   const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ; //Bohr radius of hydrogen
00463 
00464   G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
00465 
00466   G4double velocityl2 = CalculateVelocity(2,  zTarget, massIncident, energyIncident); // Scaled velocity
00467 
00468     if (verboseLevel>0) G4cout << "  velocityl2=" << velocityl2<< G4endl;
00469 
00470   const G4double l23AnalyticalApproximation= 1.25;
00471 
00472   G4double x2 = (nl*l23AnalyticalApproximation)/velocityl2;
00473     
00474     if (verboseLevel>0) G4cout << "  x2=" << x2<< G4endl;
00475 
00476   G4double electrIonizationEnergyl2=0.;
00477 
00478   if ( x2<=0.035)  electrIonizationEnergyl2= 0.75*pi*(std::log(1./(x2*x2))-1.);
00479   else
00480     {
00481       if ( x2<=3.)
00482         electrIonizationEnergyl2 =std::exp(-2.*x2)/(0.031+(0.213*std::pow(x2,0.5))+(0.005*x2)-(0.069*std::pow(x2,3./2.))+(0.324*x2*x2));
00483       else
00484         {if ( x2<=11.) electrIonizationEnergyl2 =2.*std::exp(-2.*x2)/std::pow(x2,1.6);  }
00485     }
00486 
00487   G4double hFunctionl2 =(electrIonizationEnergyl2*2.*nl)/(tetal2*std::pow(velocityl2,3)); //takes into account  the polarization effect
00488 
00489    if (verboseLevel>0) G4cout << "  hFunctionl2=" << hFunctionl2<< G4endl;
00490 
00491   G4double gFunctionl2 = (1.+(10.*velocityl2)+(45.*velocityl2*velocityl2)+(102.*std::pow(velocityl2,3.))+(331.*std::pow(velocityl2,4.))+(6.7*std::pow(velocityl2,5.))+(58.*std::pow(velocityl2,6.))+(7.8*std::pow(velocityl2,7.))+ (0.888*std::pow(velocityl2,8.)) )/std::pow(1.+velocityl2,10.);
00492   //takes into account the reduced binding effect
00493 
00494     if (verboseLevel>0) G4cout << "  gFunctionl2=" << gFunctionl2<< G4endl;
00495 
00496   G4double sigmaPSS_l2 = 1.+(((2.*zIncident)/(screenedzTarget*tetal2))*(gFunctionl2-hFunctionl2)); //Binding-polarization factor
00497 
00498     if (verboseLevel>0) G4cout << "  sigmaPSS_l2=" << sigmaPSS_l2<< G4endl;
00499 
00500   const G4double cNaturalUnit= 137.;
00501      
00502   G4double yl2Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl2/sigmaPSS_l2);
00503   
00504   G4double l2relativityCorrection = std::pow((1.+(1.1*yl2Formula*yl2Formula)),0.5)+yl2Formula; // Relativistic correction parameter
00505     
00506   G4double L2etaOverTheta2;
00507 
00508   G4double  universalFunction_l2 = 0.;
00509 
00510   G4double sigmaPSSR_l2 ;
00511  
00512   if ( velocityl2 < 20. )
00513   {
00514 
00515     L2etaOverTheta2 = (reducedEnergy*l2relativityCorrection)/((sigmaPSS_l2*tetal2)*(sigmaPSS_l2*tetal2));
00516 
00517     if ( (tetal2*sigmaPSS_l2>=0.2) && (tetal2*sigmaPSS_l2<=2.6670) && (L2etaOverTheta2>=0.1e-3) && (L2etaOverTheta2<=0.866e2) )
00518 
00519       universalFunction_l2 = FunctionFL2((tetal2*sigmaPSS_l2),L2etaOverTheta2);
00520 
00521     sigmaPSSR_l2 = (sigma0/(tetal2*sigmaPSS_l2))*universalFunction_l2;
00522 
00523     if (verboseLevel>0) G4cout << "  sigma PWBA L2 CS at low velocity range = " << sigmaPSSR_l2<< G4endl;
00524   }    
00525   else
00526   { 
00527     
00528     L2etaOverTheta2 = reducedEnergy /(tetal2*tetal2);
00529 
00530     if ( (tetal2>=0.2) && (tetal2<=2.6670) && (L2etaOverTheta2>=0.1e-3) && (L2etaOverTheta2<=0.866e2) )
00531 
00532       universalFunction_l2 = FunctionFL2((tetal2),L2etaOverTheta2);
00533    
00534     sigmaPSSR_l2 = (sigma0/tetal2)*universalFunction_l2;
00535 
00536       if (verboseLevel>0) G4cout << "  sigma PWBA L2 CS at medium and high velocity range = " << sigmaPSSR_l2<< G4endl;
00537 
00538   }
00539  
00540   G4double pssDeltal2 = (4./(systemMass*sigmaPSS_l2*tetal2))*(sigmaPSS_l2/velocityl2)*(sigmaPSS_l2/velocityl2);
00541 
00542   if (pssDeltal2>1) return 0.;
00543 
00544   G4double energyLossl2 = std::pow(1-pssDeltal2,0.5);
00545  
00546     if (verboseLevel>0) G4cout << "  energyLossl2=" << energyLossl2<< G4endl;
00547 
00548   G4double coulombDeflectionl2 
00549     =(8.*pi*zIncident/systemMass)*std::pow(tetal2*sigmaPSS_l2,-2.)*std::pow(velocityl2/sigmaPSS_l2,-3.)*(zTarget/screenedzTarget);
00550 
00551   G4double cParameterl2 = 2.*coulombDeflectionl2/(energyLossl2*(energyLossl2+1.));
00552 
00553   G4double coulombDeflectionFunction_l2 = 11.*ExpIntFunction(12,cParameterl2); //Coulomb-deflection effect correction
00554   // *** see Brandt, Phys Rev A10, p477, f25
00555   
00556     if (verboseLevel>0) G4cout << "  coulombDeflectionFunction_l2 =" << coulombDeflectionFunction_l2 << G4endl;
00557 
00558   G4double crossSection_L2 = coulombDeflectionFunction_l2 * sigmaPSSR_l2;
00559   //ECPSSR L2 -subshell cross section is estimated at perturbed-stationnairy-state(PSS)
00560   //and reduced by the energy-loss(E),the Coulomb deflection(C),and the relativity(R) effects
00561 
00562     if (verboseLevel>0) G4cout << "  crossSection_L2 =" << crossSection_L2 << G4endl;
00563 
00564   if (crossSection_L2 >= 0) 
00565   {
00566      return crossSection_L2 * barn;
00567   }
00568   else {return 0;}
00569 }
00570 
00571 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00572 
00573 
00574 G4double G4ecpssrBaseLixsModel::CalculateL3CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
00575 
00576 {
00577   if (zTarget <=13) return 0.;
00578 
00579   //this L3-CrossSection calculation method is done according to Werner Brandt and Grzegorz Lapicki, Phys.Rev.A20 N2 (1979),
00580   //and using data tables of O. Benka et al. At.Data Nucl.Data Tables Vol.22 No.3 (1978).
00581 
00582   G4NistManager* massManager = G4NistManager::Instance();
00583 
00584   G4AtomicTransitionManager*  transitionManager =  G4AtomicTransitionManager::Instance();
00585 
00586   G4double  zIncident = 0;
00587 
00588   G4Proton* aProtone = G4Proton::Proton();
00589   G4Alpha* aAlpha = G4Alpha::Alpha();
00590 
00591   if (massIncident == aProtone->GetPDGMass() )
00592 
00593    zIncident = (aProtone->GetPDGCharge())/eplus;
00594 
00595   else
00596     {
00597       if (massIncident == aAlpha->GetPDGMass())
00598 
00599           zIncident  = (aAlpha->GetPDGCharge())/eplus;
00600 
00601       else
00602         {
00603           G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateL3CrossSection : Proton or Alpha incident particles only. " << G4endl;
00604           G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
00605           return 0;
00606         }
00607     }
00608 
00609   G4double l3BindingEnergy = transitionManager->Shell(zTarget,3)->BindingEnergy();
00610 
00611   G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
00612 
00613   G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;//Mass of the system (projectile, target)
00614 
00615   const G4double zlshell= 4.15;
00616 
00617   G4double screenedzTarget = zTarget-zlshell;//Effective nuclear charge as seen by electrons in L3-subshell
00618 
00619   const G4double rydbergMeV= 13.6056923e-6;
00620 
00621   const G4double nl= 2.;
00622 
00623   G4double tetal3 = (l3BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV);//Screening parameter
00624 
00625     if (verboseLevel>0) G4cout << "  tetal3=" <<  tetal3<< G4endl;
00626 
00627   G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
00628 
00629   const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ;//Bohr radius of hydrogen
00630 
00631   G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);
00632 
00633   G4double velocityl3 = CalculateVelocity(3, zTarget, massIncident, energyIncident);// Scaled velocity
00634 
00635     if (verboseLevel>0) G4cout << "  velocityl3=" << velocityl3<< G4endl;
00636 
00637   const G4double l23AnalyticalApproximation= 1.25;
00638 
00639   G4double x3 = (nl*l23AnalyticalApproximation)/velocityl3;
00640 
00641     if (verboseLevel>0) G4cout << "  x3=" << x3<< G4endl;
00642 
00643   G4double electrIonizationEnergyl3=0.;
00644 
00645   if ( x3<=0.035)  electrIonizationEnergyl3= 0.75*pi*(std::log(1./(x3*x3))-1.);
00646     else
00647     {
00648       if ( x3<=3.) electrIonizationEnergyl3 =std::exp(-2.*x3)/(0.031+(0.213*std::pow(x3,0.5))+(0.005*x3)-(0.069*std::pow(x3,3./2.))+(0.324*x3*x3));
00649       else
00650         {
00651           if ( x3<=11.) electrIonizationEnergyl3 =2.*std::exp(-2.*x3)/std::pow(x3,1.6);}
00652     }
00653 
00654   G4double hFunctionl3 =(electrIonizationEnergyl3*2.*nl)/(tetal3*std::pow(velocityl3,3));//takes into account the polarization effect
00655 
00656     if (verboseLevel>0) G4cout << "  hFunctionl3=" << hFunctionl3<< G4endl;
00657 
00658   G4double gFunctionl3 = (1.+(10.*velocityl3)+(45.*velocityl3*velocityl3)+(102.*std::pow(velocityl3,3.))+(331.*std::pow(velocityl3,4.))+(6.7*std::pow(velocityl3,5.))+(58.*std::pow(velocityl3,6.))+(7.8*std::pow(velocityl3,7.))+ (0.888*std::pow(velocityl3,8.)) )/std::pow(1.+velocityl3,10.);
00659   //takes into account the reduced binding effect
00660 
00661     if (verboseLevel>0) G4cout << "  gFunctionl3=" << gFunctionl3<< G4endl;
00662 
00663   G4double sigmaPSS_l3 = 1.+(((2.*zIncident)/(screenedzTarget*tetal3))*(gFunctionl3-hFunctionl3));//Binding-polarization factor
00664 
00665     if (verboseLevel>0) G4cout << "sigmaPSS_l3 =" << sigmaPSS_l3<< G4endl;
00666 
00667   const G4double cNaturalUnit= 137.;
00668         
00669   G4double yl3Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl3/sigmaPSS_l3);
00670         
00671   G4double l3relativityCorrection = std::pow((1.+(1.1*yl3Formula*yl3Formula)),0.5)+yl3Formula; // Relativistic correction parameter
00672   
00673   G4double L3etaOverTheta2;
00674  
00675   G4double  universalFunction_l3 = 0.;
00676  
00677   G4double sigmaPSSR_l3;
00678 
00679   if ( velocityl3 < 20. )
00680   {
00681 
00682     L3etaOverTheta2 = (reducedEnergy* l3relativityCorrection)/((sigmaPSS_l3*tetal3)*(sigmaPSS_l3*tetal3));
00683 
00684     if ( (tetal3*sigmaPSS_l3>=0.2) && (tetal3*sigmaPSS_l3<=2.6670) && (L3etaOverTheta2>=0.1e-3) && (L3etaOverTheta2<=0.866e2) )
00685 
00686       universalFunction_l3 = 2.*FunctionFL2((tetal3*sigmaPSS_l3), L3etaOverTheta2  );
00687 
00688     sigmaPSSR_l3 = (sigma0/(tetal3*sigmaPSS_l3))*universalFunction_l3;
00689 
00690     if (verboseLevel>0) G4cout << "  sigma PWBA L3 CS at low velocity range = " << sigmaPSSR_l3<< G4endl;
00691 
00692   }
00693 
00694   else 
00695 
00696   {
00697 
00698     L3etaOverTheta2 = reducedEnergy/(tetal3*tetal3);
00699 
00700     if ( (tetal3>=0.2) && (tetal3<=2.6670) && (L3etaOverTheta2>=0.1e-3) && (L3etaOverTheta2<=0.866e2) ) 
00701      
00702       universalFunction_l3 = 2.*FunctionFL2(tetal3, L3etaOverTheta2  );
00703 
00704     sigmaPSSR_l3 = (sigma0/tetal3)*universalFunction_l3;
00705 
00706       if (verboseLevel>0) G4cout << "  sigma PWBA L3 CS at medium and high velocity range = " << sigmaPSSR_l3<< G4endl;
00707   }
00708   
00709   G4double pssDeltal3 = (4./(systemMass*sigmaPSS_l3*tetal3))*(sigmaPSS_l3/velocityl3)*(sigmaPSS_l3/velocityl3);
00710 
00711     if (verboseLevel>0) G4cout << "  pssDeltal3=" << pssDeltal3<< G4endl;
00712 
00713   if (pssDeltal3>1) return 0.;
00714 
00715   G4double energyLossl3 = std::pow(1-pssDeltal3,0.5);
00716 
00717   if (verboseLevel>0) G4cout << "  energyLossl3=" << energyLossl3<< G4endl;
00718 
00719   G4double coulombDeflectionl3 = 
00720     (8.*pi*zIncident/systemMass)*std::pow(tetal3*sigmaPSS_l3,-2.)*std::pow(velocityl3/sigmaPSS_l3,-3.)*(zTarget/screenedzTarget);
00721 
00722   G4double cParameterl3 = 2.*coulombDeflectionl3/(energyLossl3*(energyLossl3+1.));
00723 
00724   G4double coulombDeflectionFunction_l3 = 11.*ExpIntFunction(12,cParameterl3);//Coulomb-deflection effect correction
00725   // *** see Brandt, Phys Rev A10, p477, f25
00726 
00727     if (verboseLevel>0) G4cout << "  coulombDeflectionFunction_l3 =" << coulombDeflectionFunction_l3 << G4endl;
00728 
00729   G4double crossSection_L3 =  coulombDeflectionFunction_l3 * sigmaPSSR_l3;
00730   //ECPSSR L3 -subshell cross section is estimated at perturbed-stationnairy-state(PSS)
00731   //and reduced by the energy-loss(E),the Coulomb deflection(C),and the relativity(R) effects
00732 
00733     if (verboseLevel>0) G4cout << "  crossSection_L3 =" << crossSection_L3 << G4endl;
00734  
00735   if (crossSection_L3 >= 0) 
00736   {
00737     return crossSection_L3 * barn;
00738   }
00739   else {return 0;}
00740 }
00741 
00742 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00743 
00744 G4double G4ecpssrBaseLixsModel::CalculateVelocity(G4int subShell, G4int zTarget, G4double massIncident,  G4double energyIncident)
00745 
00746 {
00747 
00748   G4AtomicTransitionManager*  transitionManager =  G4AtomicTransitionManager::Instance();
00749 
00750   G4double liBindingEnergy = transitionManager->Shell(zTarget,subShell)->BindingEnergy();
00751 
00752   G4Proton* aProtone = G4Proton::Proton();
00753   G4Alpha* aAlpha = G4Alpha::Alpha();
00754 
00755   if (!((massIncident == aProtone->GetPDGMass()) || (massIncident == aAlpha->GetPDGMass())))
00756     {
00757       G4cout << "*** WARNING in G4ecpssrBaseLixsModel::CalculateVelocity : Proton or Alpha incident particles only. " << G4endl;
00758       G4cout << massIncident << ", " << aAlpha->GetPDGMass() << " (alpha)" << aProtone->GetPDGMass() << " (proton)" << G4endl;
00759       return 0;
00760     }
00761 
00762   const G4double zlshell= 4.15;
00763 
00764   G4double screenedzTarget = zTarget- zlshell;
00765 
00766   const G4double rydbergMeV= 13.6056923e-6;
00767 
00768   const G4double nl= 2.;
00769 
00770   G4double tetali = (liBindingEnergy*nl*nl)/(screenedzTarget*screenedzTarget*rydbergMeV);
00771 
00772   G4double reducedEnergy = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
00773 
00774   G4double velocity = 2.*nl*std::pow(reducedEnergy,0.5)/tetali;
00775   // *** see Brandt, Phys Rev A10, p10, f4
00776   
00777   return velocity;
00778 }
00779 
00780 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00781 
00782 G4double G4ecpssrBaseLixsModel::FunctionFL1(G4double k, G4double theta)
00783 {
00784 
00785   G4double sigma = 0.;
00786   G4double valueT1 = 0;
00787   G4double valueT2 = 0;
00788   G4double valueE21 = 0;
00789   G4double valueE22 = 0;
00790   G4double valueE12 = 0;
00791   G4double valueE11 = 0;
00792   G4double xs11 = 0;
00793   G4double xs12 = 0;
00794   G4double xs21 = 0;
00795   G4double xs22 = 0;
00796 
00797   // PROTECTION TO ALLOW INTERPOLATION AT MINIMUM AND MAXIMUM Eta/Theta2 values
00798 
00799   if (
00800        theta==8.66e-4 ||
00801        theta==8.66e-3 ||
00802        theta==8.66e-2 ||
00803        theta==8.66e-1 ||
00804        theta==8.66e+00 ||
00805        theta==8.66e+01
00806   ) theta=theta-1e-12;
00807 
00808   if (
00809        theta==1.e-4 ||
00810        theta==1.e-3 ||
00811        theta==1.e-2 ||
00812        theta==1.e-1 ||
00813        theta==1.e+00 ||
00814        theta==1.e+01
00815   ) theta=theta+1e-12;
00816 
00817   // END PROTECTION
00818 
00819     std::vector<double>::iterator t2 = std::upper_bound(dummyVec1.begin(),dummyVec1.end(), k);
00820     std::vector<double>::iterator t1 = t2-1;
00821 
00822     std::vector<double>::iterator e12 = std::upper_bound(aVecMap1[(*t1)].begin(),aVecMap1[(*t1)].end(), theta);
00823     std::vector<double>::iterator e11 = e12-1;
00824 
00825     std::vector<double>::iterator e22 = std::upper_bound(aVecMap1[(*t2)].begin(),aVecMap1[(*t2)].end(), theta);
00826     std::vector<double>::iterator e21 = e22-1;
00827 
00828     valueT1  =*t1;
00829     valueT2  =*t2;
00830     valueE21 =*e21;
00831     valueE22 =*e22;
00832     valueE12 =*e12;
00833     valueE11 =*e11;
00834 
00835     xs11 = FL1Data[valueT1][valueE11];
00836     xs12 = FL1Data[valueT1][valueE12];
00837     xs21 = FL1Data[valueT2][valueE21];
00838     xs22 = FL1Data[valueT2][valueE22];
00839 
00840   if (verboseLevel>0) 
00841     G4cout 
00842     << valueT1 << " "
00843     << valueT2 << " "
00844     << valueE11 << " "
00845     << valueE12 << " "
00846     << valueE21 << " "
00847     << valueE22 << " "
00848     << xs11 << " "
00849     << xs12 << " "
00850     << xs21 << " "
00851     << xs22 << " "
00852     << G4endl;
00853 
00854   G4double xsProduct = xs11 * xs12 * xs21 * xs22;
00855 
00856   if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
00857 
00858   if (xsProduct != 0.)
00859   {
00860     sigma = QuadInterpolator(  valueE11, valueE12,
00861                                valueE21, valueE22,
00862                                xs11, xs12,
00863                                xs21, xs22,
00864                                valueT1, valueT2,
00865                                k, theta );
00866   }
00867 
00868   return sigma;
00869 }
00870 
00871 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00872 
00873 G4double G4ecpssrBaseLixsModel::FunctionFL2(G4double k, G4double theta)
00874 {
00875 
00876   G4double sigma = 0.;
00877   G4double valueT1 = 0;
00878   G4double valueT2 = 0;
00879   G4double valueE21 = 0;
00880   G4double valueE22 = 0;
00881   G4double valueE12 = 0;
00882   G4double valueE11 = 0;
00883   G4double xs11 = 0;
00884   G4double xs12 = 0;
00885   G4double xs21 = 0;
00886   G4double xs22 = 0;
00887 
00888   // PROTECTION TO ALLOW INTERPOLATION AT MINIMUM AND MAXIMUM Eta/Theta2 values
00889 
00890   if (
00891        theta==8.66e-4 ||
00892        theta==8.66e-3 ||
00893        theta==8.66e-2 ||
00894        theta==8.66e-1 ||
00895        theta==8.66e+00 ||
00896        theta==8.66e+01
00897   ) theta=theta-1e-12;
00898 
00899   if (
00900        theta==1.e-4 ||
00901        theta==1.e-3 ||
00902        theta==1.e-2 ||
00903        theta==1.e-1 ||
00904        theta==1.e+00 ||
00905        theta==1.e+01
00906   ) theta=theta+1e-12;
00907 
00908   // END PROTECTION
00909 
00910     std::vector<double>::iterator t2 = std::upper_bound(dummyVec2.begin(),dummyVec2.end(), k);
00911     std::vector<double>::iterator t1 = t2-1;
00912 
00913     std::vector<double>::iterator e12 = std::upper_bound(aVecMap2[(*t1)].begin(),aVecMap2[(*t1)].end(), theta);
00914     std::vector<double>::iterator e11 = e12-1;
00915 
00916     std::vector<double>::iterator e22 = std::upper_bound(aVecMap2[(*t2)].begin(),aVecMap2[(*t2)].end(), theta);
00917     std::vector<double>::iterator e21 = e22-1;
00918 
00919     valueT1  =*t1;
00920     valueT2  =*t2;
00921     valueE21 =*e21;
00922     valueE22 =*e22;
00923     valueE12 =*e12;
00924     valueE11 =*e11;
00925 
00926     xs11 = FL2Data[valueT1][valueE11];
00927     xs12 = FL2Data[valueT1][valueE12];
00928     xs21 = FL2Data[valueT2][valueE21];
00929     xs22 = FL2Data[valueT2][valueE22];
00930 
00931   if (verboseLevel>0) 
00932     G4cout 
00933     << valueT1 << " "
00934     << valueT2 << " "
00935     << valueE11 << " "
00936     << valueE12 << " "
00937     << valueE21 << " "
00938     << valueE22 << " "
00939     << xs11 << " "
00940     << xs12 << " "
00941     << xs21 << " "
00942     << xs22 << " "
00943     << G4endl;
00944 
00945   G4double xsProduct = xs11 * xs12 * xs21 * xs22;
00946 
00947   if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
00948 
00949   if (xsProduct != 0.)
00950   {
00951     sigma = QuadInterpolator(  valueE11, valueE12,
00952                                valueE21, valueE22,
00953                                xs11, xs12,
00954                                xs21, xs22,
00955                                valueT1, valueT2,
00956                                k, theta );
00957   }
00958 
00959   return sigma;
00960 }
00961 
00962 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00963 
00964 G4double G4ecpssrBaseLixsModel::LinLinInterpolate(G4double e1,
00965                                                         G4double e2,
00966                                                         G4double e,
00967                                                         G4double xs1,
00968                                                         G4double xs2)
00969 {
00970   G4double value = xs1 + (xs2 - xs1)*(e - e1)/ (e2 - e1);
00971   return value;
00972 }
00973 
00974 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00975 
00976 G4double G4ecpssrBaseLixsModel::LinLogInterpolate(G4double e1,
00977                                                         G4double e2,
00978                                                         G4double e,
00979                                                         G4double xs1,
00980                                                         G4double xs2)
00981 {
00982   G4double d1 = std::log(xs1);
00983   G4double d2 = std::log(xs2);
00984   G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
00985   return value;
00986 }
00987 
00988 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00989 
00990 G4double G4ecpssrBaseLixsModel::LogLogInterpolate(G4double e1,
00991                                                         G4double e2,
00992                                                         G4double e,
00993                                                         G4double xs1,
00994                                                         G4double xs2)
00995 {
00996   G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
00997   G4double b = std::log10(xs2) - a*std::log10(e2);
00998   G4double sigma = a*std::log10(e) + b;
00999   G4double value = (std::pow(10.,sigma));
01000   return value;
01001 }
01002 
01003 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
01004 
01005 G4double G4ecpssrBaseLixsModel::QuadInterpolator(G4double e11, G4double e12,
01006                                                        G4double e21, G4double e22,
01007                                                        G4double xs11, G4double xs12,
01008                                                        G4double xs21, G4double xs22,
01009                                                        G4double t1, G4double t2,
01010                                                        G4double t, G4double e)
01011 {
01012 // Log-Log
01013   G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
01014   G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
01015   G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
01016 
01017 /*
01018 // Lin-Log
01019   G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
01020   G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
01021   G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
01022 */
01023 
01024 /*
01025 // Lin-Lin
01026   G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
01027   G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
01028   G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
01029 */
01030   return value;
01031 
01032 }
01033 

Generated on Mon May 27 17:48:05 2013 for Geant4 by  doxygen 1.4.7