EvtPolInt Class Reference

#include <EvtPolInt.hh>

List of all members.

Public Member Functions

 EvtPolInt (vector< double > xi, vector< double > yi, double x)
virtual ~EvtPolInt ()
void polynomial ()
void ratint ()
vector< double > spline ()
void splint ()
double getvalue ()
double geterror ()

Private Attributes

vector< double > vx
vector< double > vy
double xx
double value
double error
int size

Detailed Description

Definition at line 31 of file EvtPolInt.hh.

Constructor & Destructor Documentation

EvtPolInt::EvtPolInt ( vector< double >  xi,
vector< double >  yi,
double  x 
) [inline]

Definition at line 33 of file EvtPolInt.hh.

References genRecEmupikp::i, size, vx, vy, and xx.

00033                                                            {
00034     xx = x;
00035     size = xi.size();
00036     vx.clear();
00037     vy.clear();
00038     for(int i=0;i<size;i++){
00039       vx.push_back(xi[i]);
00040       vy.push_back(yi[i]);
00041     }
00042   }

virtual EvtPolInt::~EvtPolInt (  )  [inline, virtual]

Definition at line 43 of file EvtPolInt.hh.

00043 {}

Member Function Documentation

double EvtPolInt::geterror (  ) 

Definition at line 172 of file EvtPolInt.cc.

References error, and polynomial().

00172                           {
00173   polynomial();
00174   //ratint();
00175 return error;
00176 }

double EvtPolInt::getvalue (  ) 

Definition at line 165 of file EvtPolInt.cc.

References polynomial(), and value.

00165                           {
00166   polynomial();
00167   //  ratint();
00168   //  splint();
00169 return value;
00170 }

void EvtPolInt::polynomial (  ) 

Definition at line 23 of file EvtPolInt.cc.

References error, genRecEmupikp::i, ns, value, vx, vy, w, and xx.

Referenced by geterror(), and getvalue().

00023                           {
00024   double y,dy;
00025   int i,m,ns=0;
00026   double den,dif,dift,ho,hp,w;
00027 
00028   int n=vx.size();
00029   vector <double> c(n),d(n);
00030   dif = fabs(xx-vx[0]);
00031   for (i=0;i<n;i++){
00032     if ((dift = fabs(xx-vx[i]))<dif){
00033       ns=i;
00034       dif = dift;
00035     }
00036         c[i] = vy[i];
00037         d[i] = vy[i];
00038         }
00039     y = vy[ns--];
00040     for (m=1;m<n;m++){
00041       for (i=0;i<n-m;i++){
00042         ho = vx[i] -xx;
00043         hp = vx[i+m] -xx;
00044         w  = c[i+1] - d[i];
00045         if(( den = ho-hp) == 0.0) std::cout<<"Error in routine polint"<<std::endl;
00046         den = w/den;
00047         d[i] = hp*den;
00048         c[i] = ho*den;
00049       }
00050       y += (dy=(2*(ns+1)< (n-m) ? c[ns+1]: d[ns--]));
00051     }
00052   value = y;
00053   error = dy;
00054   // std::cout<<"value= "<<value<<std::endl;
00055   if(value <0) value = 0;
00056   return;
00057 }

void EvtPolInt::ratint (  ) 

Definition at line 59 of file EvtPolInt.cc.

References error, genRecEmupikp::i, ns, t(), TINY, value, vx, vy, w, and xx.

00060 {
00061         double y,dy;
00062         const double TINY=1.0e-25;
00063         int m,i,ns=0;
00064         double w,t,hh,h,dd;
00065 
00066         int n=vx.size();
00067         vector <double> c(n),d(n);
00068         hh=fabs(xx-vx[0]);
00069         for (i=0;i<n;i++) {
00070                 h=fabs(xx-vx[i]);
00071                 if (h == 0.0) {
00072                         y=vy[i];
00073                         dy=0.0;
00074                         return;
00075                 } else if (h < hh) {
00076                         ns=i;
00077                         hh=h;
00078                 }
00079                 c[i]=vy[i];
00080                 d[i]=vy[i]+TINY;
00081         }
00082         y=vy[ns--];
00083         for (m=1;m<n;m++) {
00084                 for (i=0;i<n-m;i++) {
00085                         w=c[i+1]-d[i];
00086                         h=vx[i+m]-xx;
00087                         t=(vx[i]-xx)*d[i]/h;
00088                         dd=t-c[i+1];
00089                         if (dd == 0.0) std::cout<<"Error in routine ratint"<<std::endl;
00090                         dd=w/dd;
00091                         d[i]=c[i+1]*dd;
00092                         c[i]=t*dd;
00093                 }
00094                 y += (dy=(2*(ns+1) < (n-m) ? c[ns+1] : d[ns--]));
00095         }
00096   value = y;
00097   error = dy;
00098   if(value <0) value = 0;
00099   return;
00100 }

vector< double > EvtPolInt::spline (  ) 

Definition at line 102 of file EvtPolInt.cc.

References genRecEmupikp::i, size, vx, and vy.

Referenced by splint().

00103 {
00104   vector <double> y2;
00105   y2.clear();
00106   for(int i=0;i<size;i++){y2.push_back(0.0);}
00107   double yp1=0;
00108   double ypn=0;
00109         int i,k;
00110         double p,qn,sig,un;
00111 
00112         int n=size;
00113         vector<double> u; //(n-1);
00114         if (yp1 > 0.99e30){
00115         u[0]=0.0;
00116         y2.push_back(0.0);}
00117         else {
00118         y2.push_back(-0.5);     
00119         u[0]=(3.0/(vx[1]-vx[0]))*((vy[1]-vy[0])/(vx[1]-vx[0])-yp1);
00120         }
00121         for (i=1;i<n-1;i++) {
00122                 sig=(vx[i]-vx[i-1])/(vx[i+1]-vx[i-1]);
00123                 p=sig*y2[i-1]+2.0;
00124                 double yy2=(sig-1.0)/p;
00125                 y2.push_back(yy2);
00126                 u[i]=(vy[i+1]-vy[i])/(vx[i+1]-vx[i]) - (vy[i]-vy[i-1])/(vx[i]-vx[i-1]);
00127                 u[i]=(6.0*u[i]/(vx[i+1]-vx[i-1])-sig*u[i-1])/p;
00128         }
00129         if (ypn > 0.99e30)
00130                 qn=un=0.0;
00131         else {
00132                 qn=0.5;
00133                 un=(3.0/(vx[n-1]-vx[n-2]))*(ypn-(vy[n-1]-vy[n-2])/(vx[n-1]-vx[n-2]));
00134         }
00135         y2[n-1]=(un-qn*u[n-2])/(qn*y2[n-2]+1.0);
00136         for (k=n-2;k>=0;k--)
00137         y2[k]=y2[k]*y2[k+1]+u[k];
00138         return y2;  
00139     }

void EvtPolInt::splint (  ) 

Definition at line 141 of file EvtPolInt.cc.

References spline(), value, vx, vy, and xx.

00141                       {
00142   vector <double> y2a = spline();
00143   double y;
00144         int k;
00145         double h,b,a;
00146 
00147         int n=vx.size();
00148         int klo=0;
00149         int khi=n-1;
00150         while (khi-klo > 1) {
00151                 k=(khi+klo) >> 1;
00152                 if (vx[k] > xx) khi=k;
00153                 else klo=k;
00154         }
00155         h=vx[khi]-vx[klo];
00156         if (h == 0.0) std::cout<<"Bad xa input to routine splint"<<std::endl;
00157         a=(vx[khi]-xx)/h;
00158         b=(xx-vx[klo])/h;
00159         y=a*vy[klo]+b*vy[khi]+((a*a*a-a)*y2a[klo]
00160                 +(b*b*b-b)*y2a[khi])*(h*h)/6.0;
00161   value = y;
00162   return;
00163 }

Member Data Documentation

double EvtPolInt::error [private]

Definition at line 54 of file EvtPolInt.hh.

Referenced by geterror(), polynomial(), and ratint().

int EvtPolInt::size [private]

Definition at line 55 of file EvtPolInt.hh.

Referenced by EvtPolInt(), and spline().

double EvtPolInt::value [private]

Definition at line 54 of file EvtPolInt.hh.

Referenced by getvalue(), polynomial(), ratint(), and splint().

vector<double> EvtPolInt::vx [private]

Definition at line 53 of file EvtPolInt.hh.

Referenced by EvtPolInt(), polynomial(), ratint(), spline(), and splint().

vector<double> EvtPolInt::vy [private]

Definition at line 53 of file EvtPolInt.hh.

Referenced by EvtPolInt(), polynomial(), ratint(), spline(), and splint().

double EvtPolInt::xx [private]

Definition at line 54 of file EvtPolInt.hh.

Referenced by EvtPolInt(), polynomial(), ratint(), and splint().

Generated on Tue Nov 29 23:19:13 2016 for BOSS_7.0.2 by  doxygen 1.4.7