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 "EvtGenBase/EvtPatches.hh"
#include <math.h>
#include "EvtGenBase/EvtVector4R.hh"
#include "EvtGenBase/EvtKine.hh"
#include "EvtGenBase/EvtComplex.hh"
#include "EvtGenBase/EvtCubicSpline.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtConst.hh"
#include <fstream>
#include <iostream>
#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);
}