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 //....oooOO0OOooo........oooOO0OOooo........oo 27 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 28 28 29 #include <cmath> 29 #include <cmath> 30 #include <iostream> 30 #include <iostream> 31 #include "G4ecpssrBaseKxsModel.hh" 31 #include "G4ecpssrBaseKxsModel.hh" 32 #include "globals.hh" 32 #include "globals.hh" 33 #include "G4PhysicalConstants.hh" 33 #include "G4PhysicalConstants.hh" 34 #include "G4SystemOfUnits.hh" 34 #include "G4SystemOfUnits.hh" 35 #include "G4AtomicTransitionManager.hh" 35 #include "G4AtomicTransitionManager.hh" 36 #include "G4NistManager.hh" 36 #include "G4NistManager.hh" 37 #include "G4Proton.hh" 37 #include "G4Proton.hh" 38 #include "G4Alpha.hh" 38 #include "G4Alpha.hh" 39 #include "G4SemiLogInterpolation.hh" 39 #include "G4SemiLogInterpolation.hh" 40 #include "G4Exp.hh" 40 #include "G4Exp.hh" 41 41 42 //....oooOO0OOooo........oooOO0OOooo........oo 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 43 43 44 G4ecpssrBaseKxsModel::G4ecpssrBaseKxsModel() 44 G4ecpssrBaseKxsModel::G4ecpssrBaseKxsModel() 45 { 45 { 46 verboseLevel=0; 46 verboseLevel=0; 47 47 48 // Storing C coefficients for high velocit 48 // Storing C coefficients for high velocity formula 49 G4String fileC1("pixe/uf/c1"); 49 G4String fileC1("pixe/uf/c1"); 50 tableC1 = new G4CrossSectionDataSet(new G4 50 tableC1 = new G4CrossSectionDataSet(new G4SemiLogInterpolation, 1.,1.); 51 51 52 G4String fileC2("pixe/uf/c2"); 52 G4String fileC2("pixe/uf/c2"); 53 tableC2 = new G4CrossSectionDataSet(new G4 53 tableC2 = new G4CrossSectionDataSet(new G4SemiLogInterpolation, 1.,1.); 54 54 55 G4String fileC3("pixe/uf/c3"); 55 G4String fileC3("pixe/uf/c3"); 56 tableC3 = new G4CrossSectionDataSet(new G4 56 tableC3 = new G4CrossSectionDataSet(new G4SemiLogInterpolation, 1.,1.); 57 57 58 // Storing FK data needed for medium veloc 58 // Storing FK data needed for medium velocities region 59 const char* path = G4FindDataDir("G4LEDATA 59 const char* path = G4FindDataDir("G4LEDATA"); 60 60 61 if (!path) { 61 if (!path) { 62 G4Exception("G4ecpssrBaseKxsModel::G4ecp 62 G4Exception("G4ecpssrBaseKxsModel::G4ecpssrBaseKxsModel()", "em0006", FatalException,"G4LEDATA environment variable not set" ); 63 return; 63 return; 64 } 64 } 65 65 66 std::ostringstream fileName; 66 std::ostringstream fileName; 67 fileName << path << "/pixe/uf/FK.dat"; 67 fileName << path << "/pixe/uf/FK.dat"; 68 std::ifstream FK(fileName.str().c_str()); 68 std::ifstream FK(fileName.str().c_str()); 69 69 70 if (!FK) 70 if (!FK) 71 G4Exception("G4ecpssrBaseKxsModel::G4ecp 71 G4Exception("G4ecpssrBaseKxsModel::G4ecpssrBaseKxsModel()", "em0003", FatalException,"error opening FK data file" ); 72 72 73 dummyVec.push_back(0.); 73 dummyVec.push_back(0.); 74 74 75 while(!FK.eof()) 75 while(!FK.eof()) 76 { 76 { 77 double x; 77 double x; 78 double y; 78 double y; 79 79 80 FK>>x>>y; 80 FK>>x>>y; 81 81 82 // Mandatory vector initialization 82 // Mandatory vector initialization 83 if (x != dummyVec.back()) 83 if (x != dummyVec.back()) 84 { 84 { 85 dummyVec.push_back(x); 85 dummyVec.push_back(x); 86 aVecMap[x].push_back(-1.); 86 aVecMap[x].push_back(-1.); 87 } 87 } 88 88 89 FK>>FKData[x][y]; 89 FK>>FKData[x][y]; 90 90 91 if (y != aVecMap[x].back()) aVecMap[x] 91 if (y != aVecMap[x].back()) aVecMap[x].push_back(y); 92 92 93 } 93 } 94 94 95 tableC1->LoadData(fileC1); 95 tableC1->LoadData(fileC1); 96 tableC2->LoadData(fileC2); 96 tableC2->LoadData(fileC2); 97 tableC3->LoadData(fileC3); 97 tableC3->LoadData(fileC3); 98 98 99 } 99 } 100 100 101 //....oooOO0OOooo........oooOO0OOooo........oo 101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 102 102 103 void print (G4double elem) 103 void print (G4double elem) 104 { 104 { 105 G4cout << elem << " "; 105 G4cout << elem << " "; 106 } 106 } 107 //....oooOO0OOooo........oooOO0OOooo........oo 107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 108 108 109 G4ecpssrBaseKxsModel::~G4ecpssrBaseKxsModel() 109 G4ecpssrBaseKxsModel::~G4ecpssrBaseKxsModel() 110 { 110 { 111 delete tableC1; 111 delete tableC1; 112 delete tableC2; 112 delete tableC2; 113 delete tableC3; 113 delete tableC3; 114 } 114 } 115 115 116 //....oooOO0OOooo........oooOO0OOooo........oo 116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 117 117 118 G4double G4ecpssrBaseKxsModel::ExpIntFunction( 118 G4double G4ecpssrBaseKxsModel::ExpIntFunction(G4int n,G4double x) 119 119 120 { 120 { 121 // this "ExpIntFunction" function allows fast 121 // this "ExpIntFunction" function allows fast evaluation of the n order exponential integral function En(x) 122 G4int i; 122 G4int i; 123 G4int ii; 123 G4int ii; 124 G4int nm1; 124 G4int nm1; 125 G4double a; 125 G4double a; 126 G4double b; 126 G4double b; 127 G4double c; 127 G4double c; 128 G4double d; 128 G4double d; 129 G4double del; 129 G4double del; 130 G4double fact; 130 G4double fact; 131 G4double h; 131 G4double h; 132 G4double psi; 132 G4double psi; 133 G4double ans = 0; 133 G4double ans = 0; 134 const G4double euler= 0.5772156649; 134 const G4double euler= 0.5772156649; 135 const G4int maxit= 100; 135 const G4int maxit= 100; 136 const G4double fpmin = 1.0e-30; 136 const G4double fpmin = 1.0e-30; 137 const G4double eps = 1.0e-7; 137 const G4double eps = 1.0e-7; 138 nm1=n-1; 138 nm1=n-1; 139 if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1 139 if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1))) { 140 G4cout << "*** WARNING in G4ecpssrBaseKxsM 140 G4cout << "*** WARNING in G4ecpssrBaseKxsModel::ExpIntFunction: bad arguments in ExpIntFunction" << G4endl; 141 G4cout << n << ", " << x << G4endl; 141 G4cout << n << ", " << x << G4endl; 142 } 142 } 143 else { 143 else { 144 if (n==0) ans=G4Exp(-x)/x; 144 if (n==0) ans=G4Exp(-x)/x; 145 else { 145 else { 146 if (x==0.0) ans=1.0/nm1; 146 if (x==0.0) ans=1.0/nm1; 147 else { 147 else { 148 if (x > 1.0) { 148 if (x > 1.0) { 149 b=x+n; 149 b=x+n; 150 c=1.0/fpmin; 150 c=1.0/fpmin; 151 d=1.0/b; 151 d=1.0/b; 152 h=d; 152 h=d; 153 for (i=1;i<=maxit;i++) { 153 for (i=1;i<=maxit;i++) { 154 a=-i*(nm1+i); 154 a=-i*(nm1+i); 155 b +=2.0; 155 b +=2.0; 156 d=1.0/(a*d+b); 156 d=1.0/(a*d+b); 157 c=b+a/c; 157 c=b+a/c; 158 del=c*d; 158 del=c*d; 159 h *=del; 159 h *=del; 160 if (std::fabs(del-1.0) < eps) { 160 if (std::fabs(del-1.0) < eps) { 161 ans=h*G4Exp(-x); 161 ans=h*G4Exp(-x); 162 return ans; 162 return ans; 163 } 163 } 164 } 164 } 165 } else { 165 } else { 166 ans = (nm1!=0 ? 1.0/nm1 : -std::log(x 166 ans = (nm1!=0 ? 1.0/nm1 : -std::log(x)-euler); 167 fact=1.0; 167 fact=1.0; 168 for (i=1;i<=maxit;i++) { 168 for (i=1;i<=maxit;i++) { 169 fact *=-x/i; 169 fact *=-x/i; 170 if (i !=nm1) del = -fact/(i-nm1); 170 if (i !=nm1) del = -fact/(i-nm1); 171 else { 171 else { 172 psi = -euler; 172 psi = -euler; 173 for (ii=1;ii<=nm1;ii++) psi +=1.0/ii; 173 for (ii=1;ii<=nm1;ii++) psi +=1.0/ii; 174 del=fact*(-std::log(x)+psi); 174 del=fact*(-std::log(x)+psi); 175 } 175 } 176 ans += del; 176 ans += del; 177 if (std::fabs(del) < std::fabs(ans) 177 if (std::fabs(del) < std::fabs(ans)*eps) return ans; 178 } 178 } 179 } 179 } 180 } 180 } 181 } 181 } 182 } 182 } 183 return ans; 183 return ans; 184 } 184 } 185 185 186 //....oooOO0OOooo........oooOO0OOooo........oo 186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 187 187 188 G4double G4ecpssrBaseKxsModel::CalculateCrossS 188 G4double G4ecpssrBaseKxsModel::CalculateCrossSection(G4int zTarget,G4double massIncident, G4double energyIncident) 189 189 190 { 190 { 191 // this K-CrossSection calculation method is 191 // this K-CrossSection calculation method is done according to W.Brandt and G.Lapicki, Phys.Rev.A23(1981)// 192 G4NistManager* massManager = G4NistManager:: 192 G4NistManager* massManager = G4NistManager::Instance(); 193 193 194 G4AtomicTransitionManager* transitionManage 194 G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance(); 195 195 196 G4double zIncident = 0; 196 G4double zIncident = 0; 197 G4Proton* aProtone = G4Proton::Proton(); 197 G4Proton* aProtone = G4Proton::Proton(); 198 G4Alpha* aAlpha = G4Alpha::Alpha(); 198 G4Alpha* aAlpha = G4Alpha::Alpha(); 199 199 200 if (massIncident == aProtone->GetPDGMass() ) 200 if (massIncident == aProtone->GetPDGMass() ) 201 { 201 { 202 zIncident = (aProtone->GetPDGCharge())/eplu 202 zIncident = (aProtone->GetPDGCharge())/eplus; 203 } 203 } 204 else 204 else 205 { 205 { 206 if (massIncident == aAlpha->GetPDGMass() 206 if (massIncident == aAlpha->GetPDGMass()) 207 { 207 { 208 zIncident = (aAlpha->GetPDGCharge())/eplu 208 zIncident = (aAlpha->GetPDGCharge())/eplus; 209 } 209 } 210 else 210 else 211 { 211 { 212 G4cout << "*** WARNING in G4ecpssrBaseKxsM 212 G4cout << "*** WARNING in G4ecpssrBaseKxsModel::CalculateCrossSection : we can treat only Proton or Alpha incident particles " << G4endl; 213 return 0; 213 return 0; 214 } 214 } 215 } 215 } 216 216 217 if (verboseLevel>0) G4cout << " massIncid 217 if (verboseLevel>0) G4cout << " massIncident=" << massIncident<< G4endl; 218 218 219 G4double kBindingEnergy = transitionManager- 219 G4double kBindingEnergy = transitionManager->Shell(zTarget,0)->BindingEnergy(); 220 220 221 if (verboseLevel>0) G4cout << " kBindingE 221 if (verboseLevel>0) G4cout << " kBindingEnergy=" << kBindingEnergy/eV<< G4endl; 222 222 223 G4double massTarget = (massManager->GetAtomi 223 G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2; 224 224 225 if (verboseLevel>0) G4cout << " massTarge 225 if (verboseLevel>0) G4cout << " massTarget=" << massTarget<< G4endl; 226 226 227 G4double systemMass =((massIncident*massTarg 227 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2; //the mass of the system (projectile, target) 228 228 229 if (verboseLevel>0) G4cout << " systemMas 229 if (verboseLevel>0) G4cout << " systemMass=" << systemMass<< G4endl; 230 230 231 constexpr G4double zkshell= 0.3; 231 constexpr G4double zkshell= 0.3; 232 // *** see Brandt, Phys Rev A23, p 1727 232 // *** see Brandt, Phys Rev A23, p 1727 233 233 234 G4double screenedzTarget = zTarget-zkshell; 234 G4double screenedzTarget = zTarget-zkshell; // screenedzTarget is the screened nuclear charge of the target 235 // *** see Brandt, Phys Rev A23, p 1727 235 // *** see Brandt, Phys Rev A23, p 1727 236 236 237 constexpr G4double rydbergMeV= 13.6056923e-6 237 constexpr G4double rydbergMeV= 13.6056923e-6; 238 238 239 G4double tetaK = kBindingEnergy/((screenedzT 239 G4double tetaK = kBindingEnergy/((screenedzTarget*screenedzTarget)*rydbergMeV); //tetaK denotes the reduced binding energy of the electron 240 // *** see Rice, ADANDT 20, p 504, f 2 240 // *** see Rice, ADANDT 20, p 504, f 2 241 241 242 if (verboseLevel>0) G4cout << " tetaK=" < 242 if (verboseLevel>0) G4cout << " tetaK=" << tetaK<< G4endl; 243 243 244 G4double velocity =(2./(tetaK*screenedzTarge 244 G4double velocity =(2./(tetaK*screenedzTarget))*std::pow(((energyIncident*electron_mass_c2)/(massIncident*rydbergMeV)),0.5); 245 // *** also called xiK 245 // *** also called xiK 246 // *** see Brandt, Phys Rev A23, p 1727 246 // *** see Brandt, Phys Rev A23, p 1727 247 // *** see Basbas, Phys Rev A17, p 1656, f4 247 // *** see Basbas, Phys Rev A17, p 1656, f4 248 248 249 if (verboseLevel>0) G4cout << " velocity= 249 if (verboseLevel>0) G4cout << " velocity=" << velocity<< G4endl; 250 250 251 const G4double bohrPow2Barn=(Bohr_radius*Boh 251 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ; 252 252 253 if (verboseLevel>0) G4cout << " bohrPow2B 253 if (verboseLevel>0) G4cout << " bohrPow2Barn=" << bohrPow2Barn<< G4endl; 254 254 255 G4double sigma0 = 8.*pi*(zIncident*zIncident 255 G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.); //sigma0 is the initial cross section of K shell at stable state 256 // *** see Benka, ADANDT 22, p 220, f2, for 256 // *** see Benka, ADANDT 22, p 220, f2, for protons 257 // *** see Basbas, Phys Rev A7, p 1000 257 // *** see Basbas, Phys Rev A7, p 1000 258 258 259 if (verboseLevel>0) G4cout << " sigma0=" << 259 if (verboseLevel>0) G4cout << " sigma0=" << sigma0<< G4endl; 260 260 261 const G4double kAnalyticalApproximation= 1.5 261 const G4double kAnalyticalApproximation= 1.5; 262 G4double x = kAnalyticalApproximation/veloci 262 G4double x = kAnalyticalApproximation/velocity; 263 // *** see Brandt, Phys Rev A23, p 1727 263 // *** see Brandt, Phys Rev A23, p 1727 264 // *** see Brandt, Phys Rev A20, p 469, f16 264 // *** see Brandt, Phys Rev A20, p 469, f16 in expression of h 265 265 266 if (verboseLevel>0) G4cout << " x=" << x< 266 if (verboseLevel>0) G4cout << " x=" << x<< G4endl; 267 267 268 G4double electrIonizationEnergy; 268 G4double electrIonizationEnergy; 269 // *** see Basbas, Phys Rev A17, p1665, f27 269 // *** see Basbas, Phys Rev A17, p1665, f27 270 // *** see Brandt, Phys Rev A20, p469 270 // *** see Brandt, Phys Rev A20, p469 271 // *** see Liu, Comp Phys Comm 97, p325, f A 271 // *** see Liu, Comp Phys Comm 97, p325, f A5 272 272 273 if ((0.< x) && (x <= 0.035)) 273 if ((0.< x) && (x <= 0.035)) 274 { 274 { 275 electrIonizationEnergy= 0.75*pi*(std::lo 275 electrIonizationEnergy= 0.75*pi*(std::log(1./(x*x))-1.); 276 } 276 } 277 else 277 else 278 { 278 { 279 if ( (0.035 < x) && (x <=3.)) 279 if ( (0.035 < x) && (x <=3.)) 280 { 280 { 281 electrIonizationEnergy =G4Exp(-2.*x)/(0.03 281 electrIonizationEnergy =G4Exp(-2.*x)/(0.031+(0.213*std::pow(x,0.5))+(0.005*x)-(0.069*std::pow(x,3./2.))+(0.324*x*x)); 282 } 282 } 283 283 284 else 284 else 285 { 285 { 286 if ( (3.< x) && (x<=11.)) 286 if ( (3.< x) && (x<=11.)) 287 { 287 { 288 electrIonizationEnergy =2.*G4Exp(-2.*x) 288 electrIonizationEnergy =2.*G4Exp(-2.*x)/std::pow(x,1.6); 289 } 289 } 290 290 291 else electrIonizationEnergy =0.; 291 else electrIonizationEnergy =0.; 292 } 292 } 293 } 293 } 294 294 295 if (verboseLevel>0) G4cout << " electrIon 295 if (verboseLevel>0) G4cout << " electrIonizationEnergy=" << electrIonizationEnergy<< G4endl; 296 296 297 G4double hFunction =(electrIonizationEnergy* 297 G4double hFunction =(electrIonizationEnergy*2.)/(tetaK*std::pow(velocity,3)); //hFunction represents the correction for polarization effet 298 // *** see Brandt, Phys Rev A20, p 469, f16 298 // *** see Brandt, Phys Rev A20, p 469, f16 299 299 300 if (verboseLevel>0) G4cout << " hFunction 300 if (verboseLevel>0) G4cout << " hFunction=" << hFunction<< G4endl; 301 301 302 G4double gFunction = (1.+(9.*velocity)+(31.* 302 G4double gFunction = (1.+(9.*velocity)+(31.*velocity*velocity)+(98.*std::pow(velocity,3.))+(12.*std::pow(velocity,4.))+(25.*std::pow(velocity,5.)) 303 +(4.2*std::pow(velocity,6.))+(0.515*std: 303 +(4.2*std::pow(velocity,6.))+(0.515*std::pow(velocity,7.)))/std::pow(1.+velocity,9.); //gFunction represents the correction for binding effet 304 // *** see Brandt, Phys Rev A20, p 469, f19 304 // *** see Brandt, Phys Rev A20, p 469, f19 305 305 306 if (verboseLevel>0) G4cout << " gFunction 306 if (verboseLevel>0) G4cout << " gFunction=" << gFunction<< G4endl; 307 307 308 //------------------------------------------ 308 //----------------------------------------------------------------------------------------------------------------------------- 309 309 310 G4double sigmaPSS = 1.+(((2.*zIncident)/(scr 310 G4double sigmaPSS = 1.+(((2.*zIncident)/(screenedzTarget*tetaK))*(gFunction-hFunction)); //describes the perturbed stationnairy state of the affected atomic electon 311 // *** also called dzeta 311 // *** also called dzeta 312 // *** also called epsilon 312 // *** also called epsilon 313 // *** see Basbas, Phys Rev A17, p1667, f45 313 // *** see Basbas, Phys Rev A17, p1667, f45 314 314 315 if (verboseLevel>0) G4cout << " sigmaPSS= 315 if (verboseLevel>0) G4cout << " sigmaPSS=" << sigmaPSS<< G4endl; 316 316 317 if (verboseLevel>0) G4cout << " sigmaPSS* 317 if (verboseLevel>0) G4cout << " sigmaPSS*tetaK=" << sigmaPSS*tetaK<< G4endl; 318 318 319 //------------------------------------------ 319 //---------------------------------------------------------------------------------------------------------------------------- 320 320 321 const G4double cNaturalUnit= 1/fine_structur 321 const G4double cNaturalUnit= 1/fine_structure_const; // it's the speed of light according to Atomic-Unit-System 322 322 323 if (verboseLevel>0) G4cout << " cNaturalU 323 if (verboseLevel>0) G4cout << " cNaturalUnit=" << cNaturalUnit<< G4endl; 324 324 325 G4double ykFormula=0.4*(screenedzTarget/cNat 325 G4double ykFormula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocity/sigmaPSS); 326 // *** also called yS 326 // *** also called yS 327 // *** see Brandt, Phys Rev A20, p467, f6 327 // *** see Brandt, Phys Rev A20, p467, f6 328 // *** see Brandt, Phys Rev A23, p1728 328 // *** see Brandt, Phys Rev A23, p1728 329 329 330 if (verboseLevel>0) G4cout << " ykFormula 330 if (verboseLevel>0) G4cout << " ykFormula=" << ykFormula<< G4endl; 331 331 332 G4double relativityCorrection = std::pow((1. 332 G4double relativityCorrection = std::pow((1.+(1.1*ykFormula*ykFormula)),0.5)+ykFormula;// the relativistic correction parameter 333 // *** also called mRS 333 // *** also called mRS 334 // *** see Brandt, Phys Rev A20, p467, f6 334 // *** see Brandt, Phys Rev A20, p467, f6 335 335 336 if (verboseLevel>0) G4cout << " relativit 336 if (verboseLevel>0) G4cout << " relativityCorrection=" << relativityCorrection<< G4endl; 337 337 338 G4double reducedVelocity = velocity*std::pow 338 G4double reducedVelocity = velocity*std::pow(relativityCorrection,0.5); // presents the reduced collision velocity parameter 339 // *** also called xiR 339 // *** also called xiR 340 // *** see Brandt, Phys Rev A20, p468, f7 340 // *** see Brandt, Phys Rev A20, p468, f7 341 // *** see Brandt, Phys Rev A23, p1728 341 // *** see Brandt, Phys Rev A23, p1728 342 342 343 if (verboseLevel>0) G4cout << " reducedVe 343 if (verboseLevel>0) G4cout << " reducedVelocity=" << reducedVelocity<< G4endl; 344 344 345 G4double etaOverTheta2 = (energyIncident*ele 345 G4double etaOverTheta2 = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget) 346 /(sigmaPSS*tetaK)/( 346 /(sigmaPSS*tetaK)/(sigmaPSS*tetaK); 347 // *** see Benka, ADANDT 22, p220, f4 for et 347 // *** see Benka, ADANDT 22, p220, f4 for eta 348 // then we use sigmaPSS*tetaK == epsilon*tet 348 // then we use sigmaPSS*tetaK == epsilon*tetaK 349 349 350 if (verboseLevel>0) G4cout << " etaOverTh 350 if (verboseLevel>0) G4cout << " etaOverTheta2=" << etaOverTheta2<< G4endl; 351 351 352 G4double universalFunction = 0; 352 G4double universalFunction = 0; 353 353 354 // low velocity formula 354 // low velocity formula 355 // ***************** 355 // ***************** 356 if ( velocity < 1. ) 356 if ( velocity < 1. ) 357 // OR 357 // OR 358 //if ( reducedVelocity/sigmaPSS < 1.) 358 //if ( reducedVelocity/sigmaPSS < 1.) 359 // *** see Brandt, Phys Rev A23, p1727 359 // *** see Brandt, Phys Rev A23, p1727 360 // *** reducedVelocity/sigmaPSS is also call 360 // *** reducedVelocity/sigmaPSS is also called xiR/dzeta 361 // ***************** 361 // ***************** 362 { 362 { 363 if (verboseLevel>0) G4cout << " Notice 363 if (verboseLevel>0) G4cout << " Notice : FK is computed from low velocity formula" << G4endl; 364 364 365 universalFunction = (std::pow(2.,9.)/45. 365 universalFunction = (std::pow(2.,9.)/45.)*std::pow(reducedVelocity/sigmaPSS,8.)*std::pow((1.+(1.72*(reducedVelocity/sigmaPSS)*(reducedVelocity/sigmaPSS))),-4.);// is the reduced universal cross section 366 // *** see Brandt, Phys Rev A23, p1728 366 // *** see Brandt, Phys Rev A23, p1728 367 367 368 if (verboseLevel>0) G4cout << " univers 368 if (verboseLevel>0) G4cout << " universalFunction by Brandt 1981 =" << universalFunction<< G4endl; 369 369 370 } 370 } 371 else 371 else 372 { 372 { 373 373 374 if ( etaOverTheta2 > 86.6 && (sigmaPSS*tet 374 if ( etaOverTheta2 > 86.6 && (sigmaPSS*tetaK) > 0.4 && (sigmaPSS*tetaK) < 2.9996 ) 375 { 375 { 376 // High and medium energies. Method from 376 // High and medium energies. Method from Rice ADANDT 20, p506, 1977 on tables from Benka 1978 377 377 378 if (verboseLevel>0) G4cout << " Notice 378 if (verboseLevel>0) G4cout << " Notice : FK is computed from high velocity formula" << G4endl; 379 379 380 if (verboseLevel>0) G4cout << " sigmaPS 380 if (verboseLevel>0) G4cout << " sigmaPSS*tetaK=" << sigmaPSS*tetaK << G4endl; 381 381 382 G4double C1= tableC1->FindValue(sigmaPSS 382 G4double C1= tableC1->FindValue(sigmaPSS*tetaK); 383 G4double C2= tableC2->FindValue(sigmaPSS 383 G4double C2= tableC2->FindValue(sigmaPSS*tetaK); 384 G4double C3= tableC3->FindValue(sigmaPSS 384 G4double C3= tableC3->FindValue(sigmaPSS*tetaK); 385 385 386 if (verboseLevel>0) G4cout << " C1=" 386 if (verboseLevel>0) G4cout << " C1=" << C1 << G4endl; 387 if (verboseLevel>0) G4cout << " C2=" 387 if (verboseLevel>0) G4cout << " C2=" << C2 << G4endl; 388 if (verboseLevel>0) G4cout << " C3=" 388 if (verboseLevel>0) G4cout << " C3=" << C3 << G4endl; 389 389 390 G4double etaK = (energyIncident*electron 390 G4double etaK = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget); 391 // *** see Benka, ADANDT 22, p220, f4 fo 391 // *** see Benka, ADANDT 22, p220, f4 for eta 392 392 393 if (verboseLevel>0) G4cout << " etaK= 393 if (verboseLevel>0) G4cout << " etaK=" << etaK << G4endl; 394 394 395 G4double etaT = (sigmaPSS*tetaK)*(sigmaP 395 G4double etaT = (sigmaPSS*tetaK)*(sigmaPSS*tetaK)*(86.6); // at any theta, the largest tabulated etaOverTheta2 is 86.6 396 // *** see Rice, ADANDT 20, p506 396 // *** see Rice, ADANDT 20, p506 397 397 398 if (verboseLevel>0) G4cout << " etaT= 398 if (verboseLevel>0) G4cout << " etaT=" << etaT << G4endl; 399 399 400 G4double fKT = FunctionFK((sigmaPSS*teta 400 G4double fKT = FunctionFK((sigmaPSS*tetaK),86.6)*(etaT/(sigmaPSS*tetaK)); 401 // *** see Rice, ADANDT 20, p506 401 // *** see Rice, ADANDT 20, p506 402 402 403 if (FunctionFK((sigmaPSS*tetaK),86.6)<=0 403 if (FunctionFK((sigmaPSS*tetaK),86.6)<=0.) 404 { 404 { 405 G4cout << 405 G4cout << 406 "*** WARNING in G4ecpssrBaseKxsModel:: 406 "*** WARNING in G4ecpssrBaseKxsModel::CalculateCrossSection : unable to interpolate FK function in high velocity region ! ***" << G4endl; 407 return 0; 407 return 0; 408 } 408 } 409 409 410 if (verboseLevel>0) G4cout << " Funct 410 if (verboseLevel>0) G4cout << " FunctionFK=" << FunctionFK((sigmaPSS*tetaK),86.6) << G4endl; 411 411 412 if (verboseLevel>0) G4cout << " fKT=" 412 if (verboseLevel>0) G4cout << " fKT=" << fKT << G4endl; 413 413 414 G4double GK = C2/(4*etaK) + C3/(32*etaK* 414 G4double GK = C2/(4*etaK) + C3/(32*etaK*etaK); 415 415 416 if (verboseLevel>0) G4cout << " GK=" << GK 416 if (verboseLevel>0) G4cout << " GK=" << GK << G4endl; 417 417 418 G4double GT = C2/(4*etaT) + C3/(32*etaT* 418 G4double GT = C2/(4*etaT) + C3/(32*etaT*etaT); 419 419 420 if (verboseLevel>0) G4cout << " GT=" 420 if (verboseLevel>0) G4cout << " GT=" << GT << G4endl; 421 421 422 G4double DT = fKT - C1*std::log(etaT) + 422 G4double DT = fKT - C1*std::log(etaT) + GT; 423 423 424 if (verboseLevel>0) G4cout << " DT=" 424 if (verboseLevel>0) G4cout << " DT=" << DT << G4endl; 425 425 426 G4double fKK = C1*std::log(etaK) + DT - 426 G4double fKK = C1*std::log(etaK) + DT - GK; 427 427 428 if (verboseLevel>0) G4cout << " fKK=" 428 if (verboseLevel>0) G4cout << " fKK=" << fKK << G4endl; 429 429 430 G4double universalFunction3= fKK/(etaK/t 430 G4double universalFunction3= fKK/(etaK/tetaK); 431 // *** see Rice, ADANDT 20, p505, f7 431 // *** see Rice, ADANDT 20, p505, f7 432 432 433 if (verboseLevel>0) G4cout << " unive 433 if (verboseLevel>0) G4cout << " universalFunction3=" << universalFunction3 << G4endl; 434 434 435 universalFunction=universalFunction3; 435 universalFunction=universalFunction3; 436 436 437 } 437 } 438 else if ( etaOverTheta2 >= 1.e-3 && etaOve 438 else if ( etaOverTheta2 >= 1.e-3 && etaOverTheta2 <= 86.6 && (sigmaPSS*tetaK) >= 0.4 && (sigmaPSS*tetaK) <= 2.9996 ) 439 { 439 { 440 // From Benka 1978 440 // From Benka 1978 441 441 442 if (verboseLevel>0) G4cout << " Notice 442 if (verboseLevel>0) G4cout << " Notice : FK is computed from INTERPOLATED data" << G4endl; 443 443 444 G4double universalFunction2 = FunctionFK 444 G4double universalFunction2 = FunctionFK((sigmaPSS*tetaK),etaOverTheta2); 445 445 446 if (universalFunction2<=0) 446 if (universalFunction2<=0) 447 { 447 { 448 G4cout << 448 G4cout << 449 "*** WARNING : G4ecpssrBaseKxsModel::C 449 "*** WARNING : G4ecpssrBaseKxsModel::CalculateCrossSection is unable to interpolate FK function in medium velocity region ! ***" << G4endl; 450 return 0; 450 return 0; 451 } 451 } 452 452 453 if (verboseLevel>0) G4cout << " univers 453 if (verboseLevel>0) G4cout << " universalFunction2=" << universalFunction2 << " for theta=" << sigmaPSS*tetaK << " and etaOverTheta2=" << etaOverTheta2 << G4endl; 454 454 455 universalFunction=universalFunction2; 455 universalFunction=universalFunction2; 456 } 456 } 457 457 458 } 458 } 459 459 460 //------------------------------------------ 460 //---------------------------------------------------------------------------------------------------------------------- 461 461 462 G4double sigmaPSSR = (sigma0/(sigmaPSS*tetaK 462 G4double sigmaPSSR = (sigma0/(sigmaPSS*tetaK))*universalFunction; //sigmaPSSR is the straight-line K-shell ionization cross section 463 // *** see Benka, ADANDT 22, p220, f1 463 // *** see Benka, ADANDT 22, p220, f1 464 464 465 if (verboseLevel>0) G4cout << " sigmaPSSR 465 if (verboseLevel>0) G4cout << " sigmaPSSR=" << sigmaPSSR<< G4endl; 466 466 467 //------------------------------------------ 467 //----------------------------------------------------------------------------------------------------------------------- 468 468 469 G4double pssDeltaK = (4./(systemMass*sigmaPS 469 G4double pssDeltaK = (4./(systemMass*sigmaPSS*tetaK))*(sigmaPSS/velocity)*(sigmaPSS/velocity); 470 // *** also called dzetaK*deltaK 470 // *** also called dzetaK*deltaK 471 // *** see Brandt, Phys Rev A23, p1727, f B2 471 // *** see Brandt, Phys Rev A23, p1727, f B2 472 472 473 if (verboseLevel>0) G4cout << " pssDeltaK 473 if (verboseLevel>0) G4cout << " pssDeltaK=" << pssDeltaK<< G4endl; 474 474 475 if (pssDeltaK>1) return 0.; 475 if (pssDeltaK>1) return 0.; 476 476 477 G4double energyLoss = std::pow(1-pssDeltaK,0 477 G4double energyLoss = std::pow(1-pssDeltaK,0.5); //energyLoss incorporates the straight-line energy-loss 478 // *** also called zK 478 // *** also called zK 479 // *** see Brandt, Phys Rev A23, p1727, afte 479 // *** see Brandt, Phys Rev A23, p1727, after f B2 480 480 481 if (verboseLevel>0) G4cout << " energyLos 481 if (verboseLevel>0) G4cout << " energyLoss=" << energyLoss<< G4endl; 482 482 483 G4double energyLossFunction = (std::pow(2.,- 483 G4double energyLossFunction = (std::pow(2.,-9)/8.)*((((9.*energyLoss)-1.)*std::pow(1.+energyLoss,9.))+(((9.*energyLoss)+1.)*std::pow(1.-energyLoss,9.)));//energy loss function 484 // *** also called fs 484 // *** also called fs 485 // *** see Brandt, Phys Rev A23, p1718, f7 485 // *** see Brandt, Phys Rev A23, p1718, f7 486 486 487 if (verboseLevel>0) G4cout << " energyLos 487 if (verboseLevel>0) G4cout << " energyLossFunction=" << energyLossFunction<< G4endl; 488 488 489 //------------------------------------------ 489 //---------------------------------------------------------------------------------------------------------------------------------------------- 490 490 491 G4double coulombDeflection = (4.*pi*zInciden 491 G4double coulombDeflection = (4.*pi*zIncident/systemMass)*std::pow(tetaK*sigmaPSS,-2.)*std::pow(velocity/sigmaPSS,-3.)*(zTarget/screenedzTarget); //incorporates Coulomb deflection parameter 492 // *** see Brandt, Phys Rev A23, p1727, f B3 492 // *** see Brandt, Phys Rev A23, p1727, f B3 493 493 494 if (verboseLevel>0) G4cout << " cParamete 494 if (verboseLevel>0) G4cout << " cParameter-short=" << coulombDeflection<< G4endl; 495 495 496 G4double cParameter = 2.*coulombDeflection/( 496 G4double cParameter = 2.*coulombDeflection/(energyLoss*(energyLoss+1.)); 497 // *** see Brandt, Phys Rev A23, p1727, f B4 497 // *** see Brandt, Phys Rev A23, p1727, f B4 498 498 499 if (verboseLevel>0) G4cout << " cParamete 499 if (verboseLevel>0) G4cout << " cParameter-full=" << cParameter<< G4endl; 500 500 501 G4double coulombDeflectionFunction = 9.*ExpI 501 G4double coulombDeflectionFunction = 9.*ExpIntFunction(10,cParameter); //this function describes Coulomb-deflection effect 502 // *** see Brandt, Phys Rev A23, p1727 502 // *** see Brandt, Phys Rev A23, p1727 503 503 504 if (verboseLevel>0) G4cout << " ExpIntFun 504 if (verboseLevel>0) G4cout << " ExpIntFunction(10,cParameter) =" << ExpIntFunction(10,cParameter) << G4endl; 505 505 506 if (verboseLevel>0) G4cout << " coulombDe 506 if (verboseLevel>0) G4cout << " coulombDeflectionFunction =" << coulombDeflectionFunction << G4endl; 507 507 508 //------------------------------------------ 508 //-------------------------------------------------------------------------------------------------------------------------------------------------- 509 509 510 G4double crossSection = 0; 510 G4double crossSection = 0; 511 511 512 crossSection = energyLossFunction* coulombDe 512 crossSection = energyLossFunction* coulombDeflectionFunction*sigmaPSSR; //this ECPSSR cross section is estimated at perturbed-stationnairy-state(PSS) 513 513 //and it's reduced by the energy-loss(E),the Coulomb deflection(C), 514 514 //and the relativity(R) effects 515 515 516 //------------------------------------------ 516 //-------------------------------------------------------------------------------------------------------------------------------------------------- 517 517 518 if (crossSection >= 0) { 518 if (crossSection >= 0) { 519 return crossSection * barn; 519 return crossSection * barn; 520 } 520 } 521 else {return 0;} 521 else {return 0;} 522 522 523 } 523 } 524 524 525 //....oooOO0OOooo........oooOO0OOooo........oo 525 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 526 526 527 G4double G4ecpssrBaseKxsModel::FunctionFK(G4do 527 G4double G4ecpssrBaseKxsModel::FunctionFK(G4double k, G4double theta) 528 { 528 { 529 529 530 G4double sigma = 0.; 530 G4double sigma = 0.; 531 G4double valueT1 = 0; 531 G4double valueT1 = 0; 532 G4double valueT2 = 0; 532 G4double valueT2 = 0; 533 G4double valueE21 = 0; 533 G4double valueE21 = 0; 534 G4double valueE22 = 0; 534 G4double valueE22 = 0; 535 G4double valueE12 = 0; 535 G4double valueE12 = 0; 536 G4double valueE11 = 0; 536 G4double valueE11 = 0; 537 G4double xs11 = 0; 537 G4double xs11 = 0; 538 G4double xs12 = 0; 538 G4double xs12 = 0; 539 G4double xs21 = 0; 539 G4double xs21 = 0; 540 G4double xs22 = 0; 540 G4double xs22 = 0; 541 541 542 // PROTECTION TO ALLOW INTERPOLATION AT MINI 542 // PROTECTION TO ALLOW INTERPOLATION AT MINIMUM AND MAXIMUM EtaK/Theta2 values 543 // (in particular for FK computation at 8.66 543 // (in particular for FK computation at 8.66EXX for high velocity formula) 544 544 545 if ( 545 if ( 546 theta==8.66e-3 || 546 theta==8.66e-3 || 547 theta==8.66e-2 || 547 theta==8.66e-2 || 548 theta==8.66e-1 || 548 theta==8.66e-1 || 549 theta==8.66e+0 || 549 theta==8.66e+0 || 550 theta==8.66e+1 550 theta==8.66e+1 551 ) theta=theta-1e-12; 551 ) theta=theta-1e-12; 552 552 553 if ( 553 if ( 554 theta==1.e-3 || 554 theta==1.e-3 || 555 theta==1.e-2 || 555 theta==1.e-2 || 556 theta==1.e-1 || 556 theta==1.e-1 || 557 theta==1.e+00 || 557 theta==1.e+00 || 558 theta==1.e+01 558 theta==1.e+01 559 ) theta=theta+1e-12; 559 ) theta=theta+1e-12; 560 560 561 // END PROTECTION 561 // END PROTECTION 562 562 563 auto t2 = std::upper_bound(dummyVec.begin(), 563 auto t2 = std::upper_bound(dummyVec.begin(),dummyVec.end(), k); 564 auto t1 = t2-1; 564 auto t1 = t2-1; 565 565 566 auto e12 = std::upper_bound(aVecMap[(*t1)].b 566 auto e12 = std::upper_bound(aVecMap[(*t1)].begin(),aVecMap[(*t1)].end(), theta); 567 auto e11 = e12-1; 567 auto e11 = e12-1; 568 568 569 auto e22 = std::upper_bound(aVecMap[(*t2)].b 569 auto e22 = std::upper_bound(aVecMap[(*t2)].begin(),aVecMap[(*t2)].end(), theta); 570 auto e21 = e22-1; 570 auto e21 = e22-1; 571 571 572 valueT1 =*t1; 572 valueT1 =*t1; 573 valueT2 =*t2; 573 valueT2 =*t2; 574 valueE21 =*e21; 574 valueE21 =*e21; 575 valueE22 =*e22; 575 valueE22 =*e22; 576 valueE12 =*e12; 576 valueE12 =*e12; 577 valueE11 =*e11; 577 valueE11 =*e11; 578 578 579 xs11 = FKData[valueT1][valueE11]; 579 xs11 = FKData[valueT1][valueE11]; 580 xs12 = FKData[valueT1][valueE12]; 580 xs12 = FKData[valueT1][valueE12]; 581 xs21 = FKData[valueT2][valueE21]; 581 xs21 = FKData[valueT2][valueE21]; 582 xs22 = FKData[valueT2][valueE22]; 582 xs22 = FKData[valueT2][valueE22]; 583 583 584 G4double xsProduct = xs11 * xs12 * xs21 * xs 584 G4double xsProduct = xs11 * xs12 * xs21 * xs22; 585 585 586 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) 586 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.); 587 587 588 if (xsProduct != 0.) 588 if (xsProduct != 0.) 589 { 589 { 590 sigma = QuadInterpolator( valueE11, value 590 sigma = QuadInterpolator( valueE11, valueE12, 591 valueE21, valueE22, 591 valueE21, valueE22, 592 xs11, xs12, 592 xs11, xs12, 593 xs21, xs22, 593 xs21, xs22, 594 valueT1, valueT2, 594 valueT1, valueT2, 595 k, theta ); 595 k, theta ); 596 } 596 } 597 597 598 return sigma; 598 return sigma; 599 } 599 } 600 600 601 //....oooOO0OOooo........oooOO0OOooo........oo 601 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 602 602 603 G4double G4ecpssrBaseKxsModel::LinLogInterpola 603 G4double G4ecpssrBaseKxsModel::LinLogInterpolate(G4double e1, 604 G4double e2, 604 G4double e2, 605 G4double e, 605 G4double e, 606 G4double xs1, 606 G4double xs1, 607 G4double xs2) 607 G4double xs2) 608 { 608 { 609 G4double d1 = std::log(xs1); 609 G4double d1 = std::log(xs1); 610 G4double d2 = std::log(xs2); 610 G4double d2 = std::log(xs2); 611 G4double value = G4Exp(d1 + (d2 - d1)*(e - e 611 G4double value = G4Exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1)); 612 return value; 612 return value; 613 } 613 } 614 614 615 //....oooOO0OOooo........oooOO0OOooo........oo 615 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 616 616 617 G4double G4ecpssrBaseKxsModel::LogLogInterpola 617 G4double G4ecpssrBaseKxsModel::LogLogInterpolate(G4double e1, 618 G4double e2, 618 G4double e2, 619 G4double e, 619 G4double e, 620 G4double xs1, 620 G4double xs1, 621 G4double xs2) 621 G4double xs2) 622 { 622 { 623 G4double a = (std::log10(xs2)-std::log10(xs1 623 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1)); 624 G4double b = std::log10(xs2) - a*std::log10( 624 G4double b = std::log10(xs2) - a*std::log10(e2); 625 G4double sigma = a*std::log10(e) + b; 625 G4double sigma = a*std::log10(e) + b; 626 G4double value = (std::pow(10.,sigma)); 626 G4double value = (std::pow(10.,sigma)); 627 return value; 627 return value; 628 } 628 } 629 629 630 //....oooOO0OOooo........oooOO0OOooo........oo 630 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 631 631 632 G4double G4ecpssrBaseKxsModel::QuadInterpolato 632 G4double G4ecpssrBaseKxsModel::QuadInterpolator(G4double e11, G4double e12, 633 G4double e21, G4double e22, 633 G4double e21, G4double e22, 634 G4double xs11, G4double xs1 634 G4double xs11, G4double xs12, 635 G4double xs21, G4double xs2 635 G4double xs21, G4double xs22, 636 G4double t1, G4double t2, 636 G4double t1, G4double t2, 637 G4double t, G4double e) 637 G4double t, G4double e) 638 { 638 { 639 // Log-Log 639 // Log-Log 640 G4double interpolatedvalue1 = LogLogInterpol 640 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12); 641 G4double interpolatedvalue2 = LogLogInterpol 641 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22); 642 G4double value = LogLogInterpolate(t1, t2, t 642 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2); 643 643 644 return value; 644 return value; 645 } 645 } 646 646