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