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