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

Last change on this file since 1197 was 1197, checked in by garnier, 15 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.