#include <EvtCubicSpline.hh>

Public Member Functions
EvtCubicSpline &	operator= (const EvtCubicSpline &)

	EvtCubicSpline (const EvtVector4R &p4_p, const EvtVector4R &p4_d1, const EvtVector4R &p4_d2)

virtual	~EvtCubicSpline ()

const EvtVector4R &	p4_p ()

const EvtVector4R &	p4_d1 ()

const EvtVector4R &	p4_d2 ()

double	amplitude ()

double	theta ()

EvtComplex	resAmpl ()

Static Public Member Functions
static void	setParams (const vector< double > x, const vector< double > ym, const vector< double > yp)

static void	setParams (const int n, const double x, const double ym, const double *yp)

static void	fprime (const vector< double > &x, const vector< ControlPoint > &y, const int n, ControlPoint &yp2)

static void	fprime (const vector< double > &x, const vector< ControlPoint > &y, ControlPoint &yp2)

static int	find_bin (double mass1, const vector< double > &x, const int n)

static bool	Complex_Derivative (const vector< double > &x, const vector< ControlPoint > &y, const int n, vector< ControlPoint > &y2)

Static Public Attributes
static int	_nPoints = 0

static vector< double >	_mHHLimits

static vector< ControlPoint >	_yvalues

static vector< ControlPoint >	_y2values

Detailed Description

Definition at line 42 of file EvtCubicSpline.hh.

Constructor & Destructor Documentation

◆ EvtCubicSpline()

EvtCubicSpline::EvtCubicSpline	(	const EvtVector4R &	p4_p,
		const EvtVector4R &	p4_d1,
		const EvtVector4R &	p4_d2 )

Definition at line 61 of file EvtCubicSpline.cc.

                      : 
  _p4_p(p4_p),_p4_d1(p4_d1), _p4_d2(p4_d2)
  //  _m1a(m1a), _m1b(m1b), _g1(g1),
  //  _m2a(m2a), _m2b(m2b), _g2(g2)
{
//    _nPoints = 0;
}

◆ ~EvtCubicSpline()

EvtCubicSpline::~EvtCubicSpline ( )

virtual

Definition at line 38 of file EvtCubicSpline.cc.

38{}

Member Function Documentation

◆ amplitude()

double EvtCubicSpline::amplitude ( )

inline

Definition at line 71 of file EvtCubicSpline.hh.

71{ return _ampl; }

◆ Complex_Derivative()

bool EvtCubicSpline::Complex_Derivative	(	const vector< double > &	x,
		const vector< ControlPoint > &	y,
		const int	n,
		vector< ControlPoint > &	y2 )

static

Definition at line 117 of file EvtCubicSpline.cc.

                                                                                                                                    {
    int i,k;
    ControlPoint *u = new ControlPoint[n];
    double sig,p,qn,un;
    ControlPoint yp1, ypn; 
/*    yp1.real = 2.*(y[1].real - y[0].real) / (x[1] - x[0]);
    yp1.imag = 2.*(y[1].imag - y[0].imag) / (x[1] - x[0]);
    ypn.real = 2.*(y[n-1].real - y[n-2].real) / (x[n-1] - x[n-2]);
    ypn.imag = 2.*(y[n-1].imag - y[n-2].imag) / (x[n-1] - x[n-2]);*/
    fprime( x, y,yp1);
    fprime( x, y,n,ypn);
 
    /* The lower boundary condition is set either to be "natural" or else to have specified first derivative*/
    if(yp1.real > 0.99e30) {
        y2[0].real = 0.;
        u[0].real = 0.;
    }
    else{
        y2[0].real=-0.5;
        u[0].real=(3.0/(x[1]-x[0]))*((y[1].real-y[0].real)/(x[1]-x[0])-yp1.real);
    }
    if(yp1.imag > 0.99e30) {
        y2[0].imag = 0.;
        u[0].imag = 0.;
    }
    else{
        y2[0].imag=-0.5;
        u[0].imag=(3.0/(x[1]-x[0]))*((y[1].imag-y[0].imag)/(x[1]-x[0])-yp1.imag);
    }
    
    
/* This is the decomposition loop of the tridiagonal algorithm. y2 and u are used for temporary storage of the decomposed factors*/
    
    for(i=1;i<n-1;i++){
        sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
        p=sig*y2[i-1].real+2.0;
        y2[i].real=(sig-1.0)/p;
        u[i].real=(y[i+1].real-y[i].real)/(x[i+1]-x[i]) - (y[i].real-y[i-1].real)/(x[i]-x[i-1]);
        u[i].real=(6.0*u[i].real/(x[i+1]-x[i-1])-sig*u[i-1].real)/p;
        p=sig*y2[i-1].imag+2.0;
        y2[i].imag=(sig-1.0)/p;
        u[i].imag=(y[i+1].imag-y[i].imag)/(x[i+1]-x[i]) - (y[i].imag-y[i-1].imag)/(x[i]-x[i-1]);
        u[i].imag=(6.0*u[i].imag/(x[i+1]-x[i-1])-sig*u[i-1].imag)/p;
    }
    
    /* The upper boundary condition is set either to be "natural" or else to have specified first derivative*/
    
    if(ypn.real > 0.99e30) {
        qn = 0.;
        un = 0.;
    }
    else{
        qn=0.5;
        un=(3.0/(x[n-1]-x[n-2]))*(ypn.real-(y[n-1].real-y[n-2].real)/(x[n-1]-x[n-2]));
    }
    y2[n-1].real=(un-qn*u[n-2].real)/(qn*y2[n-2].real+1.0);
    if(ypn.imag > 0.99e30) {
        qn = 0.;
        un = 0.;
    }
    else{
        qn=0.5;
        un=(3.0/(x[n-1]-x[n-2]))*(ypn.imag-(y[n-1].imag-y[n-2].imag)/(x[n-1]-x[n-2]));
    }
    y2[n-1].imag=(un-qn*u[n-2].imag)/(qn*y2[n-2].imag+1.0);
    
    /* This is the backsubstitution loop of the tridiagonal algorithm */
    
    for(k=n-2;k>=0;k--) {
        y2[k].real=y2[k].real*y2[k+1].real+u[k].real;
        y2[k].imag=y2[k].imag*y2[k+1].imag+u[k].imag;
    }
    delete [] u;
    return true;
}

Referenced by setParams().

◆ find_bin()

static int EvtCubicSpline::find_bin	(	double	mass1,
		const vector< double > &	x,
		const int	n )

inlinestatic

Definition at line 106 of file EvtCubicSpline.hh.

                                                                         {
      int mhi = 0;
      for (;mhi<n;){
          if (mass1<x[mhi]) break;
          mhi++;
      }
      return mhi==0?1:(mhi==n?n-1:mhi);
  }

Referenced by resAmpl().

◆ fprime() [1/2]

static void EvtCubicSpline::fprime	(	const vector< double > &	x,
		const vector< ControlPoint > &	y,
		const int	n,
		ControlPoint &	yp2 )

inlinestatic

Definition at line 81 of file EvtCubicSpline.hh.

                                                                                                                   {
      double s1 = x[n-1] - x[n-2],
       s2 = x[n-1] - x[n-3],
       s3 = x[n-1] - x[n-4],
      s12 = s1-s2,
      s13 = s1-s3,
      s23 = s2-s3;
      yp2.real = -( s1*s2/( s13*s23*s3 ) )*y[n-4].real+( s1*s3/( s12*s2*s23 ) )*y[n-3].real
           -( s2*s3/( s1*s12*s13 ) )*y[n-2].real+( 1./s1+1./s2+1./s3 )*y[n-1].real;
      yp2.imag = -( s1*s2/( s13*s23*s3 ) )*y[n-4].imag+( s1*s3/( s12*s2*s23 ) )*y[n-3].imag
           -( s2*s3/( s1*s12*s13 ) )*y[n-2].imag+( 1./s1+1./s2+1./s3 )*y[n-1].imag;
  }

Referenced by Complex_Derivative().

◆ fprime() [2/2]

static void EvtCubicSpline::fprime	(	const vector< double > &	x,
		const vector< ControlPoint > &	y,
		ControlPoint &	yp2 )

inlinestatic

Definition at line 93 of file EvtCubicSpline.hh.

                                                                                                      {
      double s1 = x[0]-x[1],
     s2 = x[0]-x[2],
     s3 = x[0]-x[3],
     s12 = s1-s2,
     s13 = s1-s3,
     s23 = s2-s3;
 
     yp2.real = -( s1*s2/( s13*s23*s3 ) )*y[3].real+( s1*s3/( s12*s2*s23 ) )*y[2].real
           -( s2*s3/( s1*s12*s13 ) )*y[1].real+( 1./s1+1./s2+1./s3 )*y[0].real;
     yp2.imag = -( s1*s2/( s13*s23*s3 ) )*y[3].imag+( s1*s3/( s12*s2*s23 ) )*y[2].imag
           -( s2*s3/( s1*s12*s13 ) )*y[1].imag+( 1./s1+1./s2+1./s3 )*y[0].imag;
  }

◆ operator=()

EvtCubicSpline & EvtCubicSpline::operator= ( const EvtCubicSpline & n )

Definition at line 42 of file EvtCubicSpline.cc.

{
  if ( &n == this ) return *this;
  _p4_p = n._p4_p;
  _p4_d1 = n._p4_d1;
  _p4_d2 = n._p4_d2;
  _ampl = n._ampl;
  _theta = n._theta;
/*  _nPoints = n._nPoints;
  for (int i=0;i<_nPoints;i++){
      _mHHLimits.push_back(n._mHHLimits[i]);
      _yvalues.push_back(n._yvalues[i]);
      _y2values.push_back(n._y2values[i]);
  }*/
  return  *this;
}

◆ p4_d1()

const EvtVector4R & EvtCubicSpline::p4_d1 ( )

inline

Definition at line 66 of file EvtCubicSpline.hh.

66{ return _p4_d1; }

◆ p4_d2()

const EvtVector4R & EvtCubicSpline::p4_d2 ( )

inline

Definition at line 67 of file EvtCubicSpline.hh.

67{ return _p4_d2; }

◆ p4_p()

const EvtVector4R & EvtCubicSpline::p4_p ( )

inline

Definition at line 65 of file EvtCubicSpline.hh.

65{ return _p4_p; }

◆ resAmpl()

EvtComplex EvtCubicSpline::resAmpl ( )

Definition at line 195 of file EvtCubicSpline.cc.

                                   {
 
  //   EvtComplex ampl(cos(_theta*pi180inv), sin(_theta*pi180inv));
  //   ampl *= _ampl;
  if (_nPoints ==0) return EvtComplex(0,0);
 
  // SCALARS ONLY
  double mAB = (_p4_d1+_p4_d2).mass();
  int khiAB = find_bin(mAB, _mHHLimits, _nPoints);
  int kloAB = khiAB -1 ;
  double dmHH, aa, bb, aa3, bb3;
  double pwa_coefs_real_kloAB = _yvalues[kloAB].real;
  double pwa_coefs_imag_kloAB = _yvalues[kloAB].imag;
  double pwa_coefs_real_khiAB = _yvalues[khiAB].real;
  double pwa_coefs_imag_khiAB = _yvalues[khiAB].imag;
  double pwa_coefs_prime_real_kloAB = _y2values[kloAB].real;
  double pwa_coefs_prime_imag_kloAB = _y2values[kloAB].imag;
  double pwa_coefs_prime_real_khiAB = _y2values[khiAB].real;
  double pwa_coefs_prime_imag_khiAB = _y2values[khiAB].imag;
  dmHH = _mHHLimits[khiAB] - _mHHLimits[kloAB];
  aa = ( _mHHLimits[khiAB] - mAB )/dmHH;
  bb = 1 - aa;
  aa3 = aa * aa * aa; bb3 = bb * bb * bb;
  double ret_real = aa * pwa_coefs_real_kloAB + bb * pwa_coefs_real_khiAB + ((aa3 - aa)*pwa_coefs_prime_real_kloAB + (bb3 - bb) * pwa_coefs_prime_real_khiAB) * (dmHH*dmHH)/6.0;
  double ret_imag = aa * pwa_coefs_imag_kloAB + bb * pwa_coefs_imag_khiAB + ((aa3 - aa)*pwa_coefs_prime_imag_kloAB + (bb3 - bb) * pwa_coefs_prime_imag_khiAB) * (dmHH*dmHH)/6.0;
  return EvtComplex(ret_real,ret_imag);
}

Referenced by EvtDDalitz::decay().

◆ setParams() [1/2]

void EvtCubicSpline::setParams	(	const int	n,
		const double *	x,
		const double *	ym,
		const double *	yp )

static

Definition at line 84 of file EvtCubicSpline.cc.

                                                                                              {
    vector<double> vx, vym, vyp;
    for (int i=0;i<n;i++){
        vx.push_back(x[i]);
        vym.push_back(ym[i]);
        vyp.push_back(yp[i]);
    }
    setParams(vx, vym, vyp);
}

◆ setParams() [2/2]

void EvtCubicSpline::setParams	(	const vector< double >	x,
		const vector< double >	ym,
		const vector< double >	yp )

static

Definition at line 71 of file EvtCubicSpline.cc.

                                                                                                      {
    if (_nPoints) return;
    double e1, e2, e3;
    for (int i=0;i<x.size();i++){
        e1 = x[i]; e2 = ym[i]; e3 = yp[i];
        _mHHLimits.push_back(e1);
        _yvalues.push_back(ControlPoint(e2*cos(e3),e2*sin(e3)));
        _y2values.push_back(ControlPoint(0,0));
        _nPoints ++; 
    }
    Complex_Derivative(_mHHLimits, _yvalues, _nPoints, _y2values);
}

Referenced by EvtDDalitz::decay(), and setParams().

◆ theta()

double EvtCubicSpline::theta ( )

inline

Definition at line 74 of file EvtCubicSpline.hh.

74{ return _theta; }

Member Data Documentation

◆ _mHHLimits

vector< double > EvtCubicSpline::_mHHLimits

static

Definition at line 46 of file EvtCubicSpline.hh.

Referenced by resAmpl(), and setParams().

◆ _nPoints

int EvtCubicSpline::_nPoints = 0

static

Definition at line 45 of file EvtCubicSpline.hh.

Referenced by resAmpl(), and setParams().

◆ _y2values

vector< ControlPoint > EvtCubicSpline::_y2values

static

Definition at line 48 of file EvtCubicSpline.hh.

Referenced by resAmpl(), and setParams().

◆ _yvalues

vector< ControlPoint > EvtCubicSpline::_yvalues

static

Definition at line 47 of file EvtCubicSpline.hh.

Referenced by resAmpl(), and setParams().

The documentation for this class was generated from the following files:

Public Member Functions

Static Public Member Functions

Static Public Attributes

Detailed Description

Constructor & Destructor Documentation

◆ EvtCubicSpline()

◆ ~EvtCubicSpline()

Member Function Documentation

◆ amplitude()

◆ Complex_Derivative()

◆ find_bin()

◆ fprime() [1/2]

◆ fprime() [2/2]

◆ operator=()

◆ p4_d1()

◆ p4_d2()

◆ p4_p()

◆ resAmpl()

◆ setParams() [1/2]

◆ setParams() [2/2]

◆ theta()

Member Data Documentation

◆ _mHHLimits

◆ _nPoints

◆ _y2values

◆ _yvalues