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