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 "EvtGenBase/EvtVector4R.hh"
#include <vector>
#include <string>
 
using std::vector;
using std::string;
 
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