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