source: trunk/source/processes/electromagnetic/lowenergy/src/G4ecpssrKCrossSection.cc@ 1197

Last change on this file since 1197 was 1197, checked in by garnier, 16 years ago

nvx fichiers dans CVS
File size: 23.2 KB
Line 
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//$Id: G4ecpssrKCrossSection.cc,v 1.7 2009/11/11 09:14:53 mantero Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
28//
29// Author: Haifa Ben Abdelouahed
30//
31//
32// History:
33// -----------
34// 21 Apr 2008 H. Ben Abdelouahed 1st implementation
35// 21 Apr 2008 MGP Major revision according to a design iteration
36// 21 Apr 2009 ALF Some correction for compatibility to G4VShellCrossSection
37// and changed name to G4ecpssrKCrossSection
38// 11 Nov 2009 ALF update and code cleaning for the Dec Release
39//
40// -------------------------------------------------------------------
41// Class description:
42// Low Energy Electromagnetic Physics, Cross section, p ionisation, K shell
43// Further documentation available from http://www.ge.infn.it/geant4/lowE
44// -------------------------------------------------------------------
45
46
47#include "globals.hh"
48#include "G4ecpssrKCrossSection.hh"
49#include "G4AtomicTransitionManager.hh"
50#include "G4NistManager.hh"
51#include "G4Proton.hh"
52#include "G4Alpha.hh"
53#include <math.h>
54#include <iostream>
55
56#include "G4SemiLogInterpolation.hh"
57
58
59G4ecpssrKCrossSection::G4ecpssrKCrossSection()
60{
61
62 // Storing FK data needed for medium velocities region
63
64 char *path = getenv("G4LEDATA");
65
66 if (!path)
67 G4Exception("G4ecpssrKCrossSection::CalculateCrossSection: G4LEDATA environment variable not set");
68
69 std::ostringstream fileName;
70 fileName << path << "/pixe/uf/FK.dat";
71 std::ifstream FK(fileName.str().c_str());
72
73 if (!FK) G4Exception("G4ecpssrKCrossSection::CalculateCrossSection: 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 // Storing C coefficients for high velocity formula
98
99 G4String fileC1("pixe/uf/c1");
100 tableC1 = new G4DNACrossSectionDataSet(new G4SemiLogInterpolation, 1.,1.);
101 tableC1->LoadData(fileC1);
102
103 G4String fileC2("pixe/uf/c2");
104 tableC2 = new G4DNACrossSectionDataSet(new G4SemiLogInterpolation, 1.,1.);
105 tableC2->LoadData(fileC2);
106
107 G4String fileC3("pixe/uf/c3");
108 tableC3 = new G4DNACrossSectionDataSet(new G4SemiLogInterpolation, 1.,1.);
109 tableC3->LoadData(fileC3);
110
111 //
112
113 verboseLevel=0;
114}
115
116//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
117
118void print (G4double elem)
119{
120 G4cout << elem << " ";
121}
122//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123
124G4ecpssrKCrossSection::~G4ecpssrKCrossSection()
125{
126
127 delete tableC1;
128 delete tableC2;
129 delete tableC3;
130
131}
132
133//---------------------------------this "ExpIntFunction" function allows fast evaluation of the n order exponential integral function En(x)------
134
135G4double G4ecpssrKCrossSection::ExpIntFunction(G4int n,G4double x)
136
137{
138 G4int i;
139 G4int ii;
140 G4int nm1;
141 G4double a;
142 G4double b;
143 G4double c;
144 G4double d;
145 G4double del;
146 G4double fact;
147 G4double h;
148 G4double psi;
149 G4double ans = 0;
150 const G4double euler= 0.5772156649;
151 const G4int maxit= 100;
152 const G4double fpmin = 1.0e-30;
153 const G4double eps = 1.0e-7;
154 nm1=n-1;
155 if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1)))
156 G4cout << "bad arguments in ExpIntFunction" << G4endl;
157 else {
158 if (n==0) ans=std::exp(-x)/x;
159 else {
160 if (x==0.0) ans=1.0/nm1;
161 else {
162 if (x > 1.0) {
163 b=x+n;
164 c=1.0/fpmin;
165 d=1.0/b;
166 h=d;
167 for (i=1;i<=maxit;i++) {
168 a=-i*(nm1+i);
169 b +=2.0;
170 d=1.0/(a*d+b);
171 c=b+a/c;
172 del=c*d;
173 h *=del;
174 if (std::fabs(del-1.0) < eps) {
175 ans=h*std::exp(-x);
176 return ans;
177 }
178 }
179 } else {
180 ans = (nm1!=0 ? 1.0/nm1 : -std::log(x)-euler);
181 fact=1.0;
182 for (i=1;i<=maxit;i++) {
183 fact *=-x/i;
184 if (i !=nm1) del = -fact/(i-nm1);
185 else {
186 psi = -euler;
187 for (ii=1;ii<=nm1;ii++) psi +=1.0/ii;
188 del=fact*(-std::log(x)+psi);
189 }
190 ans += del;
191 if (std::fabs(del) < std::fabs(ans)*eps) return ans;
192 }
193 }
194 }
195 }
196 }
197return ans;
198}
199//-----------------------------------------------------------------------------------------------------------
200
201
202
203G4double G4ecpssrKCrossSection::CalculateCrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
204
205 //this K-CrossSection calculation method is done according to W.Brandt and G.Lapicki, Phys.Rev.A23(1981)//
206
207{
208
209 //if (energyIncident < 150 *keV) {return 0;}
210
211 G4NistManager* massManager = G4NistManager::Instance();
212
213 G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
214
215 G4double zIncident = 0;
216 G4Proton* aProtone = G4Proton::Proton();
217 G4Alpha* aAlpha = G4Alpha::Alpha();
218
219 if (massIncident == aProtone->GetPDGMass() )
220 {
221
222 zIncident = (aProtone->GetPDGCharge())/eplus;
223
224 // G4cout << "zincident:" << zIncident << G4endl;
225 }
226 else
227 {
228 if (massIncident == aAlpha->GetPDGMass())
229 {
230
231 zIncident = (aAlpha->GetPDGCharge())/eplus;
232
233 // G4cout << "zincident:" << zIncident << G4endl;
234 }
235 else
236 {
237 G4cout << "we can treat only Proton or Alpha incident particles " << G4endl;
238 return 0;
239 }
240 }
241
242 ///////////////////////////////////////////
243 //from here I will substitute this version with Seb's one.
244 ///////////////////////////////////////////
245
246
247 if (verboseLevel>0) G4cout << " massIncident=" << massIncident<< G4endl;
248
249 G4double kBindingEnergy = transitionManager->Shell(zTarget,0)->BindingEnergy();
250 if (verboseLevel>0) G4cout << " kBindingEnergy=" << kBindingEnergy/eV<< G4endl;
251
252 G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
253 if (verboseLevel>0) G4cout << " massTarget=" << massTarget<< G4endl;
254
255 G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2; //the mass of the system (projectile, target)
256 if (verboseLevel>0) G4cout << " systemMass=" << systemMass<< G4endl;
257
258 const G4double zkshell= 0.3;
259
260 G4double screenedzTarget = zTarget-zkshell; // screenedzTarget is the screened nuclear charge of the target
261
262 const G4double rydbergMeV= 13.6e-6;
263
264 G4double tetaK = kBindingEnergy/((screenedzTarget*screenedzTarget)*rydbergMeV); //tetaK denotes the reduced binding energy of the electron
265 if (verboseLevel>0) G4cout << " tetaK=" << tetaK<< G4endl;
266
267 G4double velocity =(2./(tetaK*screenedzTarget))*std::pow(((energyIncident*electron_mass_c2)/(massIncident*rydbergMeV)),0.5);
268 if (verboseLevel>0) G4cout << " velocity=" << velocity<< G4endl;
269
270 const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ;
271 if (verboseLevel>0) G4cout << " bohrPow2Barn=" << bohrPow2Barn<< G4endl;
272
273 G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.); //sigma0 is the initial cross section of K shell at stable state
274 if (verboseLevel>0) G4cout << " sigma0=" << sigma0<< G4endl;
275
276 const G4double kAnalyticalApproximation= 1.5;
277
278 G4double x = kAnalyticalApproximation/velocity;
279 if (verboseLevel>0) G4cout << " x=" << x<< G4endl;
280
281
282
283
284 G4double electrIonizationEnergy;
285
286 if ( x<0.035)
287 {
288 electrIonizationEnergy= 0.75*pi*(std::log(1./(x*x))-1.);
289 }
290 else
291 {
292 if ( x<3.)
293 {
294 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));
295 }
296
297 else
298 {
299 electrIonizationEnergy =2.*std::exp(-2.*x)/std::pow(x,1.6); }
300 }
301
302 if (verboseLevel>0) G4cout << " electrIonizationEnergy=" << electrIonizationEnergy<< G4endl;
303
304 G4double hFunction =(electrIonizationEnergy*2.)/(tetaK*std::pow(velocity,3)); //hFunction represents the correction for polarization effet
305
306 if (verboseLevel>0) G4cout << " hFunction=" << hFunction<< G4endl;
307
308 G4double gFunction = (1.+(9.*velocity)+(31.*velocity*velocity)+(98.*std::pow(velocity,3.))+(12.*std::pow(velocity,4.))+(25.*std::pow(velocity,5.))
309 +(4.2*std::pow(velocity,6.))+(0.515*std::pow(velocity,7.)))/std::pow(1.+velocity,9.); //gFunction represents the correction for binding effet
310 if (verboseLevel>0) G4cout << " gFunction=" << gFunction<< G4endl;
311
312 //-----------------------------------------------------------------------------------------------------------------------------
313
314 G4double sigmaPSS = 1.+(((2.*zIncident)/(screenedzTarget*tetaK))*(gFunction-hFunction)); //describes the perturbed stationnairy state of the affected atomic electon
315 if (verboseLevel>0) G4cout << " sigmaPSS=" << sigmaPSS<< G4endl;
316 if (verboseLevel>0) G4cout << " sigmaPSS*tetaK=" << sigmaPSS*tetaK<< G4endl;
317
318 //----------------------------------------------------------------------------------------------------------------------------
319
320 const G4double cNaturalUnit= 1/fine_structure_const; // it's the speed of light according to Atomic-Unit-System
321 if (verboseLevel>0) G4cout << " cNaturalUnit=" << cNaturalUnit<< G4endl;
322
323 G4double ykFormula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocity/sigmaPSS);
324 if (verboseLevel>0) G4cout << " ykFormula=" << ykFormula<< G4endl;
325
326 G4double relativityCorrection = std::pow((1.+(1.1*ykFormula*ykFormula)),0.5)+ykFormula;// the relativistic correction parameter
327 if (verboseLevel>0) G4cout << " relativityCorrection=" << relativityCorrection<< G4endl;
328
329 G4double reducedVelocity = velocity*std::pow(relativityCorrection,0.5); // presents the reduced collision velocity parameter
330 if (verboseLevel>0) G4cout << " reducedVelocity=" << reducedVelocity<< G4endl;
331
332 G4double etaOverTheta2 = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget)
333 /(sigmaPSS*tetaK)/(sigmaPSS*tetaK);
334 if (verboseLevel>0) G4cout << " etaOverTheta2=" << etaOverTheta2<< G4endl;
335
336 // low velocity formula
337
338 G4double 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
339 if (verboseLevel>0) G4cout << " universalFunction by Brandt 1981 =" << universalFunction<< G4endl;
340
341 // Alternative formula by Rice 1977, closer to tabulated data than the above function from Brandt 1981
342
343 G4double x_ = etaOverTheta2;
344 if (verboseLevel>0) G4cout << " x_=" << x_ << G4endl;
345
346 G4double b0 = pow(2.,17)/45.;
347 if (verboseLevel>0) G4cout << " b0=" << b0 << G4endl;
348
349 G4double b1 = 19*(sigmaPSS*tetaK) - 480./11.;
350 if (verboseLevel>0) G4cout << " b1=" << b1 << G4endl;
351
352 G4double b2 = (720./7.)*(23*(sigmaPSS*tetaK)*(sigmaPSS*tetaK)/11 - 97*(sigmaPSS*tetaK)/9. + 160./13);
353 if (verboseLevel>0) G4cout << " b2=" << b2 << G4endl;
354
355 universalFunction = b0*x_*x_*x_*x_*(1+b1*x_+b2*x_*x_);
356
357 if (verboseLevel>0) G4cout << " universalFunction by Rice 1977 =" << universalFunction<< G4endl;
358
359
360 if ( etaOverTheta2 < 0.01 )
361 {
362 if (verboseLevel>0) G4cout << " Notice : FK is computed from low velocity formula" << G4endl;
363
364 }
365
366 else
367
368 {
369
370 if ( etaOverTheta2 > 95 && (sigmaPSS*tetaK) > 0.4 && (sigmaPSS*tetaK) < 1.7 )
371 {
372 // From Rice 1977
373
374 if (verboseLevel>0) G4cout << " Notice : FK is computed from high velocity formula" << G4endl;
375
376 if (verboseLevel>0) G4cout << " sigmaPSS*tetaK=" << sigmaPSS*tetaK << G4endl;
377
378 G4double C1= tableC1->FindValue(sigmaPSS*tetaK);
379 G4double C2= tableC2->FindValue(sigmaPSS*tetaK);
380 G4double C3= tableC3->FindValue(sigmaPSS*tetaK);
381 if (verboseLevel>0) G4cout << " C1=" << C1 << G4endl;
382 if (verboseLevel>0) G4cout << " C2=" << C2 << G4endl;
383 if (verboseLevel>0) G4cout << " C3=" << C3 << G4endl;
384
385 G4double etaK = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
386 if (verboseLevel>0) G4cout << " etaK=" << etaK << G4endl;
387
388 G4double etaT = (sigmaPSS*tetaK)*(sigmaPSS*tetaK)*(95.); // at any theta, the largest tabulated etaOverTheta2 is 95
389 if (verboseLevel>0) G4cout << " etaT=" << etaT << G4endl;
390
391 G4double fKT = FunctionFK((sigmaPSS*tetaK),95.)*(etaT/(sigmaPSS*tetaK));
392 if (FunctionFK((sigmaPSS*tetaK),95.)<=0.) G4cout <<
393 "*** WARNING : G4ecpssrCrossSection::CalculateCrossSection is unable to interpolate FK function in high velocity region ! ***" << G4endl;
394 if (verboseLevel>0) G4cout << " FunctionFK=" << FunctionFK((sigmaPSS*tetaK),95.) << G4endl;
395 if (verboseLevel>0) G4cout << " fKT=" << fKT << G4endl;
396
397 G4double GK = C2/(4*etaK) + C3/(32*etaK*etaK);
398 if (verboseLevel>0) G4cout << " GK=" << GK << G4endl;
399 G4double GT = C2/(4*etaT) + C3/(32*etaT*etaT);
400 if (verboseLevel>0) G4cout << " GT=" << GT << G4endl;
401
402 G4double DT = fKT - C1*std::log(etaT) + GT;
403 if (verboseLevel>0) G4cout << " DT=" << DT << G4endl;
404
405 G4double fKK = C1*std::log(etaK) + DT - GK;
406 if (verboseLevel>0) G4cout << " fKK=" << fKK << G4endl;
407
408 G4double universalFunction3= fKK/(etaK/tetaK);
409 if (verboseLevel>0) G4cout << " universalFunction3=" << universalFunction3 << G4endl;
410
411 universalFunction=universalFunction3;
412
413 }
414 else
415 {
416 // From Rice 1977
417
418 if (verboseLevel>0) G4cout << " Notice : FK is computed from INTERPOLATED data" << G4endl;
419
420 G4double universalFunction2 = FunctionFK((sigmaPSS*tetaK),etaOverTheta2);
421 if (universalFunction2<=0) G4cout <<
422 "*** WARNING : G4ecpssrCrossSection::CalculateCrossSection is unable to interpolate FK function in medium velocity region ! ***" << G4endl;
423
424 if (verboseLevel>0) G4cout << " universalFunction2=" << universalFunction2 << " for theta=" << sigmaPSS*tetaK << " and etaOverTheta2=" << etaOverTheta2 << G4endl;
425
426 universalFunction=universalFunction2;
427 }
428
429 }
430
431 //----------------------------------------------------------------------------------------------------------------------
432
433 G4double sigmaPSSR = (sigma0/(sigmaPSS*tetaK))*universalFunction; //sigmaPSSR is the straight-line K-shell ionization cross section
434 if (verboseLevel>0) G4cout << " sigmaPSSR=" << sigmaPSSR<< G4endl;
435
436 //-----------------------------------------------------------------------------------------------------------------------
437
438 G4double pssDeltaK = (4./(systemMass*sigmaPSS*tetaK))*(sigmaPSS/velocity)*(sigmaPSS/velocity);
439 if (verboseLevel>0) G4cout << " pssDeltaK=" << pssDeltaK<< G4endl;
440
441 G4double energyLoss = std::pow(1-pssDeltaK,0.5); //energyLoss incorporates the straight-line energy-loss
442 if (verboseLevel>0) G4cout << " energyLoss=" << energyLoss<< G4endl;
443
444 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
445 if (verboseLevel>0) G4cout << " energyLossFunction=" << energyLossFunction<< G4endl;
446
447 //----------------------------------------------------------------------------------------------------------------------------------------------
448
449 G4double coulombDeflection = (4.*pi*zIncident/systemMass)*std::pow(tetaK*sigmaPSS,-2.)*std::pow(velocity/sigmaPSS,-3.)*(zTarget/screenedzTarget); //incorporates Coulomb deflection parameter
450
451 G4double cParameter = 2.*coulombDeflection/(energyLoss*(energyLoss+1.));
452
453 if (verboseLevel>0) G4cout << " cParameter=" << cParameter<< G4endl;
454
455 G4double coulombDeflectionFunction = 9.*ExpIntFunction(10,cParameter); //this function describes Coulomb-deflection effect
456 if (verboseLevel>0) G4cout << " ExpIntFunction(10,cParameter) =" << ExpIntFunction(10,cParameter) << G4endl;
457 if (verboseLevel>0) G4cout << " coulombDeflectionFunction =" << coulombDeflectionFunction << G4endl;
458
459 //--------------------------------------------------------------------------------------------------------------------------------------------------
460
461 G4double crossSection = 0;
462
463 crossSection = energyLossFunction* coulombDeflectionFunction*sigmaPSSR; //this ECPSSR cross section is estimated at perturbed-stationnairy-state(PSS)
464 //and it's reduced by the energy-loss(E),the Coulomb deflection(C),
465 //and the relativity(R) effects
466
467 //--------------------------------------------------------------------------------------------------------------------------------------------------
468
469 if (crossSection >= 0) {
470 return crossSection;
471 }
472 else {return 0;}
473}
474
475//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
476
477G4double G4ecpssrKCrossSection::FunctionFK(G4double k, G4double theta)
478{
479
480 G4double sigma = 0.;
481 G4double valueT1 = 0;
482 G4double valueT2 = 0;
483 G4double valueE21 = 0;
484 G4double valueE22 = 0;
485 G4double valueE12 = 0;
486 G4double valueE11 = 0;
487 G4double xs11 = 0;
488 G4double xs12 = 0;
489 G4double xs21 = 0;
490 G4double xs22 = 0;
491
492 // PROTECTION TO ALLOW INTERPOLATION AT MINIMUM AND MAXIMUM EtaK/Theta2 values
493 // (in particular for FK computation at 95 for high velocity formula)
494
495 if (
496 theta==9.5e-2 ||
497 theta==9.5e-1 ||
498 theta==9.5e+00 ||
499 theta==9.5e+01
500 ) theta=theta-1e-12;
501
502 if (
503 theta==1.e-2 ||
504 theta==1.e-1 ||
505 theta==1.e+00 ||
506 theta==1.e+01
507 ) theta=theta+1e-12;
508
509 // END PROTECTION
510
511 {
512 std::vector<double>::iterator t2 = std::upper_bound(dummyVec.begin(),dummyVec.end(), k);
513 std::vector<double>::iterator t1 = t2-1;
514
515 std::vector<double>::iterator e12 = std::upper_bound(aVecMap[(*t1)].begin(),aVecMap[(*t1)].end(), theta);
516 std::vector<double>::iterator e11 = e12-1;
517
518 std::vector<double>::iterator e22 = std::upper_bound(aVecMap[(*t2)].begin(),aVecMap[(*t2)].end(), theta);
519 std::vector<double>::iterator e21 = e22-1;
520
521 valueT1 =*t1;
522 valueT2 =*t2;
523 valueE21 =*e21;
524 valueE22 =*e22;
525 valueE12 =*e12;
526 valueE11 =*e11;
527
528 xs11 = FKData[valueT1][valueE11];
529 xs12 = FKData[valueT1][valueE12];
530 xs21 = FKData[valueT2][valueE21];
531 xs22 = FKData[valueT2][valueE22];
532
533/*
534verboseLevel=1;
535
536if (verboseLevel>0) G4cout << "x1= " << valueT1 << G4endl;
537if (verboseLevel>0) G4cout << " vector of y for x1" << G4endl;
538for_each (aVecMap[(*t1)].begin(),aVecMap[(*t1)].end(), print);
539
540if (verboseLevel>0) G4cout << G4endl;
541
542if (verboseLevel>0) G4cout << "x2= " << valueT2 << G4endl;
543if (verboseLevel>0) G4cout << " vector of y for x2" << G4endl;
544for_each (aVecMap[(*t2)].begin(),aVecMap[(*t2)].end(), print);
545
546if (verboseLevel>0) G4cout << G4endl;
547
548if (verboseLevel>0) G4cout
549<< " "
550<< valueT1 << " "
551<< valueT2 << " "
552<< valueE11 << " "
553<< valueE12 << " "
554<< valueE21<< " "
555<< valueE22 << " "
556<< xs11 << " "
557<< xs12 << " "
558<< xs21 << " "
559<< xs22 << " "
560<< G4endl;
561//verboseLevel=0;
562*/
563
564}
565
566 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
567
568 if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
569
570 if (xsProduct != 0.)
571 {
572 sigma = QuadInterpolator( valueE11, valueE12,
573 valueE21, valueE22,
574 xs11, xs12,
575 xs21, xs22,
576 valueT1, valueT2,
577 k, theta );
578 }
579
580 return sigma;
581}
582
583//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
584
585G4double G4ecpssrKCrossSection::LinLogInterpolate(G4double e1,
586 G4double e2,
587 G4double e,
588 G4double xs1,
589 G4double xs2)
590{
591 G4double d1 = std::log(xs1);
592 G4double d2 = std::log(xs2);
593 G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
594 return value;
595}
596
597//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
598
599G4double G4ecpssrKCrossSection::LogLogInterpolate(G4double e1,
600 G4double e2,
601 G4double e,
602 G4double xs1,
603 G4double xs2)
604{
605 G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
606 G4double b = std::log10(xs2) - a*std::log10(e2);
607 G4double sigma = a*std::log10(e) + b;
608 G4double value = (std::pow(10.,sigma));
609 return value;
610}
611
612//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
613
614G4double G4ecpssrKCrossSection::QuadInterpolator(G4double e11, G4double e12,
615 G4double e21, G4double e22,
616 G4double xs11, G4double xs12,
617 G4double xs21, G4double xs22,
618 G4double t1, G4double t2,
619 G4double t, G4double e)
620{
621// Log-Log
622/*
623 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
624 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
625 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
626*/
627
628// Lin-Log
629 G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
630 G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
631 G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
632 return value;
633}
634
635
636
637
638
Note: See TracBrowser for help on using the repository browser.