source: trunk/source/processes/electromagnetic/lowenergy/src/G4AnalyticalEcpssrKCrossSection.cc @ 1316

Last change on this file since 1316 was 1316, checked in by garnier, 14 years ago

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