EvtCubicSpline.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
//      Copyright (C) 1998      Caltech, UCSB
//
// Module: EvtCubicSpline.cc
//
// Description: resonance-defining class
//
// Modification history:
//
//    ponyisi  18 Feb 2008  created
//
//------------------------------------------------------------------------
//
#include "EvtCubicSpline.hh"
#include "EvtComplex.hh"
#include "EvtConst.hh"
#include "EvtKine.hh"
#include "EvtPatches.hh"
#include "EvtReport.hh"
#include "EvtVector4R.hh"
#include <fstream>
#include <iostream>
#include <math.h>
#include <stdlib.h>
 
int                  EvtCubicSpline::_nPoints = 0;
vector<double>       EvtCubicSpline::_mHHLimits;
vector<ControlPoint> EvtCubicSpline::_yvalues;
vector<ControlPoint> EvtCubicSpline::_y2values;
 
EvtCubicSpline::~EvtCubicSpline() {}
 
// operator
 
EvtCubicSpline& EvtCubicSpline::operator=( const EvtCubicSpline& n ) {
  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;
}
 
// constructor
 
EvtCubicSpline::EvtCubicSpline( const EvtVector4R& p4_p, const EvtVector4R& p4_d1,
                                const EvtVector4R& p4_d2 )
    : _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;
}
 
void EvtCubicSpline::setParams( const vector<double> x, const vector<double> ym,
                                const vector<double> yp ) {
  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 );
}
 
void EvtCubicSpline::setParams( const int n, const double* x, const double* ym,
                                const double* yp ) {
  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 );
}
 
/*
void EvtCubicSpline::setParamsFromFile(const string swaveparfile, const double scale){
    if (_nPoints) return;
    string location = getenv("BESEVTGENROOT");
    location += "/share/"; location += swaveparfile;
    std::ifstream file(location.c_str());
    std::cout<<"Reading swave parameters from file "<<location<<std::endl;
    if(file==0){std::cout<<" EvtCubicSpline::EvtCubicSpline: swave parameter file not
available"<<std::endl;abort();} double e1, e2, e3; while(!file.eof()){ file>>e1>>e2>>e3; if
(e1<0) break; _mHHLimits.push_back(e1); e2 *= scale;
        _yvalues.push_back(ControlPoint(e2*cos(e3),e2*sin(e3)));
        _y2values.push_back(ControlPoint(0,0));
        _nPoints ++;
        e1 = -1;
    }
    file.close();
    Complex_Derivative(_mHHLimits, _yvalues, _nPoints, _y2values);
}*/
 
bool EvtCubicSpline::Complex_Derivative( const vector<double>&       x,
                                         const vector<ControlPoint>& y, const int n,
                                         vector<ControlPoint>& y2 ) {
  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;
}
 
// amplitude function
 
EvtComplex EvtCubicSpline::resAmpl() {
 
  //   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 );
}