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