/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Generator/BesEvtGen/BesEvtGen-00-03-58/src/EvtGen/EvtGenBase/EvtPolInt.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of models developed at BES collaboration
//      based on the EvtGen framework.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/BesCopyright
//      Copyright (A) 2006      Ping Rong-Gang, Pang Cai-Ying@IHEP
//
// Module: EvtEvtPloInt.cc
//
// Description: Routine to deal with polynomial interpolation 
//
// Modification history:
//
//    Ping R.-G.    December, 2010   Module created
//
//------------------------------------------------------------------------

#include "EvtGenBase/EvtPolInt.hh"

void EvtPolInt::polynomial(){
  double y,dy;
  int i,m,ns=0;
  double den,dif,dift,ho,hp,w;

  int n=vx.size();
  vector <double> c(n),d(n);
  dif = fabs(xx-vx[0]);
  for (i=0;i<n;i++){
    if ((dift = fabs(xx-vx[i]))<dif){
      ns=i;
      dif = dift;
    }
        c[i] = vy[i];
        d[i] = vy[i];
        }
    y = vy[ns--];
    for (m=1;m<n;m++){
      for (i=0;i<n-m;i++){
        ho = vx[i] -xx;
        hp = vx[i+m] -xx;
        w  = c[i+1] - d[i];
        if(( den = ho-hp) == 0.0) std::cout<<"Error in routine polint"<<std::endl;
        den = w/den;
        d[i] = hp*den;
        c[i] = ho*den;
      }
      y += (dy=(2*(ns+1)< (n-m) ? c[ns+1]: d[ns--]));
    }
  value = y;
  error = dy;
  // std::cout<<"value= "<<value<<std::endl;
  if(value <0) value = 0;
  return;
}

void EvtPolInt::ratint()
{
        double y,dy;
        const double TINY=1.0e-25;
        int m,i,ns=0;
        double w,t,hh,h,dd;

        int n=vx.size();
        vector <double> c(n),d(n);
        hh=fabs(xx-vx[0]);
        for (i=0;i<n;i++) {
                h=fabs(xx-vx[i]);
                if (h == 0.0) {
                        y=vy[i];
                        dy=0.0;
                        return;
                } else if (h < hh) {
                        ns=i;
                        hh=h;
                }
                c[i]=vy[i];
                d[i]=vy[i]+TINY;
        }
        y=vy[ns--];
        for (m=1;m<n;m++) {
                for (i=0;i<n-m;i++) {
                        w=c[i+1]-d[i];
                        h=vx[i+m]-xx;
                        t=(vx[i]-xx)*d[i]/h;
                        dd=t-c[i+1];
                        if (dd == 0.0) std::cout<<"Error in routine ratint"<<std::endl;
                        dd=w/dd;
                        d[i]=c[i+1]*dd;
                        c[i]=t*dd;
                }
                y += (dy=(2*(ns+1) < (n-m) ? c[ns+1] : d[ns--]));
        }
  value = y;
  error = dy;
  if(value <0) value = 0;
  return;
}

vector <double> EvtPolInt::spline()
{
  vector <double> y2;
  y2.clear();
  for(int i=0;i<size;i++){y2.push_back(0.0);}
  double yp1=0;
  double ypn=0;
        int i,k;
        double p,qn,sig,un;

        int n=size;
        vector<double> u; //(n-1);
        if (yp1 > 0.99e30){
        u[0]=0.0;
        y2.push_back(0.0);}
        else {
        y2.push_back(-0.5);     
        u[0]=(3.0/(vx[1]-vx[0]))*((vy[1]-vy[0])/(vx[1]-vx[0])-yp1);
        }
        for (i=1;i<n-1;i++) {
                sig=(vx[i]-vx[i-1])/(vx[i+1]-vx[i-1]);
                p=sig*y2[i-1]+2.0;
                double yy2=(sig-1.0)/p;
                y2.push_back(yy2);
                u[i]=(vy[i+1]-vy[i])/(vx[i+1]-vx[i]) - (vy[i]-vy[i-1])/(vx[i]-vx[i-1]);
                u[i]=(6.0*u[i]/(vx[i+1]-vx[i-1])-sig*u[i-1])/p;
        }
        if (ypn > 0.99e30)
                qn=un=0.0;
        else {
                qn=0.5;
                un=(3.0/(vx[n-1]-vx[n-2]))*(ypn-(vy[n-1]-vy[n-2])/(vx[n-1]-vx[n-2]));
        }
        y2[n-1]=(un-qn*u[n-2])/(qn*y2[n-2]+1.0);
        for (k=n-2;k>=0;k--)
        y2[k]=y2[k]*y2[k+1]+u[k];
        return y2;  
    }

void EvtPolInt::splint(){
  vector <double> y2a = spline();
  double y;
        int k;
        double h,b,a;

        int n=vx.size();
        int klo=0;
        int khi=n-1;
        while (khi-klo > 1) {
                k=(khi+klo) >> 1;
                if (vx[k] > xx) khi=k;
                else klo=k;
        }
        h=vx[khi]-vx[klo];
        if (h == 0.0) std::cout<<"Bad xa input to routine splint"<<std::endl;
        a=(vx[khi]-xx)/h;
        b=(xx-vx[klo])/h;
        y=a*vy[klo]+b*vy[khi]+((a*a*a-a)*y2a[klo]
                +(b*b*b-b)*y2a[khi])*(h*h)/6.0;
  value = y;
  return;
}

double EvtPolInt::getvalue(){
  polynomial();
  //  ratint();
  //  splint();
return value;
}

double EvtPolInt::geterror(){
  polynomial();
  //ratint();
return error;
}

---