EvtCubicSpline.hh

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: EvtGen/EvtCubicSpline.hh
//
// Description:resonance-defining class
//
// Modification history:
//
//    ponyisi  18 Feb 2008  created
//
//------------------------------------------------------------------------
 
#ifndef EVTCUBICSPLINE_HH
#define EVTCUBICSPLINE_HH
 
#include "EvtVector4R.hh"
#include <string>
#include <vector>
 
using std::string;
using std::vector;
 
class EvtComplex;
 
// Helper class
struct ControlPoint {
  double real;
  double imag;
  ControlPoint( double r = 0, double i = 0. ) {
    real = r;
    imag = i;
  }
};
 
// class declaration
 
class EvtCubicSpline {
public:
  static int                  _nPoints;
  static vector<double>       _mHHLimits;
  static vector<ControlPoint> _yvalues;
  static vector<ControlPoint> _y2values;
  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 setParamsFromFile(const string swaveparfile, const double scale = 1.);
 
  // operator
  EvtCubicSpline& operator=( const EvtCubicSpline& );
 
  // constructor with all information about the resonance
  EvtCubicSpline( const EvtVector4R& p4_p, const EvtVector4R& p4_d1,
                  const EvtVector4R& p4_d2 );
 
  // destructor
  virtual ~EvtCubicSpline();
 
  // accessors
  // return 4-momenta of the particles involved
  inline const EvtVector4R& p4_p() { return _p4_p; }
  inline const EvtVector4R& p4_d1() { return _p4_d1; }
  inline const EvtVector4R& p4_d2() { return _p4_d2; }
 
  // return amplitude
  inline double amplitude() { return _ampl; }
 
  // return theta
  inline double theta() { return _theta; }
 
  // functions
 
  // calculate amplitude for this resonance
  EvtComplex resAmpl();
 
  inline static void fprime( const vector<double>& x, const vector<ControlPoint>& y,
                             const int n, ControlPoint& yp2 ) {
    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;
  }
  inline static void fprime( const vector<double>& x, const vector<ControlPoint>& y,
                             ControlPoint& yp2 ) {
    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;
  }
  static int find_bin( double mass1, const vector<double>& x, const int n ) {
    int mhi = 0;
    for ( ; mhi < n; )
    {
      if ( mass1 < x[mhi] ) break;
      mhi++;
    }
    return mhi == 0 ? 1 : ( mhi == n ? n - 1 : mhi );
  }
  static bool Complex_Derivative( const vector<double>& x, const vector<ControlPoint>& y,
                                  const int n, vector<ControlPoint>& y2 );
 
private:
  EvtVector4R _p4_p, _p4_d1, _p4_d2;
  double      _ampl, _theta;
 
  //      double _m1a, _m1b, _g1;
  //      double _m2a, _m2b, _g2;
};
 
#endif