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