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