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