TSpline3_red Class Reference

CW: Reduced TSpline3 class. More...

#include <Splines/SplineStructs.h>

Inheritance diagram for TSpline3_red:

Collaboration diagram for TSpline3_red:

Public Member Functions
	TSpline3_red ()
	Empty constructor.
	TSpline3_red (TSpline3 *&spline, SplineInterpolation InterPolation=kTSpline3)
	The constructor that takes a TSpline3 pointer and copies in to memory.
	TSpline3_red (M3::float_t *X, M3::float_t *Y, M3::int_t N, M3::float_t **P)
	constructor taking parameters
void	SetFunc (TSpline3 *&spline, SplineInterpolation InterPolation=kTSpline3)
	Set the function.
virtual	~TSpline3_red ()
	Empty destructor.
int	FindX (double x)
	Find the segment relevant to this variation in x See root/hist/hist/src/TSpline3FindX(double) or samplePDFND....::FindSplineSegment.
double	Eval (double var) override
	CW: Evaluate the weight from a variation.
M3::int_t	GetNp () override
	CW: Get the number of points.
void	GetKnot (int i, M3::float_t &xtmp, M3::float_t &ytmp)
void	GetCoeff (int segment, M3::float_t &x, M3::float_t &y, M3::float_t &b, M3::float_t &c, M3::float_t &d)
	CW: Get the coefficient of a given segment.
TSpline3 *	ConstructTSpline3 ()
	CW: Make a TSpline3 from the reduced splines.
void	Print () override
	Print detailed info.
Public Member Functions inherited from TResponseFunction_red
	TResponseFunction_red ()
	Empty constructor.
virtual	~TResponseFunction_red ()
	Empty destructor.
virtual double	Eval (const double var)=0
	Evaluate a variation.
virtual void	Print ()=0
	KS: Printer.
virtual M3::int_t	GetNp ()=0
	DL: Get number of points.

Protected Attributes
M3::int_t	nPoints
	Number of points/knot in TSpline3.
M3::float_t **	Par
	Always uses a third order polynomial, so hard-code the number of coefficients in implementation.
M3::float_t *	XPos
	Positions of each x for each knot.
M3::float_t *	YResp
	y-value for each knot

Detailed Description

CW: Reduced TSpline3 class.

Definition at line 267 of file SplineStructs.h.

Constructor & Destructor Documentation

TSpline3_red() [1/3]

TSpline3_red::TSpline3_red ( )

inline

Empty constructor.

Definition at line 271 of file SplineStructs.h.

                 : TResponseFunction_red() {
    nPoints = 0;
    Par = NULL;
    XPos = NULL;
    YResp = NULL;
  }

TSpline3_red() [2/3]

TSpline3_red::TSpline3_red ( TSpline3 *&  spline, SplineInterpolation  InterPolation = kTSpline3  )

inline

The constructor that takes a TSpline3 pointer and copies in to memory.

Definition at line 279 of file SplineStructs.h.

                                                                                 : TResponseFunction_red() {
    Par = NULL;
    XPos = NULL;
    YResp = NULL;
    SetFunc(spline, InterPolation);
  }

TSpline3_red() [3/3]

TSpline3_red::TSpline3_red ( M3::float_t *  X, M3::float_t *  Y, M3::int_t  N, M3::float_t **  P  )

inline

constructor taking parameters

Definition at line 287 of file SplineStructs.h.

                                                                     : TResponseFunction_red() {
    nPoints = N;
    // Save the parameters for each knot
    Par = new M3::float_t*[nPoints];
    // Save the positions of the knots
    XPos = new M3::float_t[nPoints];
    // Save the y response at each knot
    YResp = new M3::float_t[nPoints];
    for(int j = 0; j < N; ++j){
      Par[j] = new M3::float_t[3];
      Par[j][0] = P[j][0];
      Par[j][1] = P[j][1];
      Par[j][2] = P[j][2];
      XPos[j] = X[j];
      YResp[j] = Y[j];
 
      if((Par[j][0] == -999) | (Par[j][1] ==-999) | (Par[j][2] ==-999) | (XPos[j] ==-999) | (YResp[j] ==-999)){
        MACH3LOG_ERROR("******************* Bad parameter values when constructing TSpline3_red *********************");
        MACH3LOG_ERROR("Passed values (i, x, y, b, c, d): {}, {}, {}, {}, {}, {}", j, X[j], Y[j], P[j][0], P[j][1], P[j][2]);
        MACH3LOG_ERROR("Set values (i, x, y, b, c, d): {}, {}, {}, {}, {}, {}", j, XPos[j], YResp[j], Par[j][0], Par[j][1], Par[j][2]);
        MACH3LOG_ERROR("*********************************************************************************************");
      }
    }
  }

~TSpline3_red()

virtual TSpline3_red::~TSpline3_red ( )

inlinevirtual

Empty destructor.

Definition at line 580 of file SplineStructs.h.

                          {
    if(Par != NULL) {
      for (int i = 0; i < nPoints; ++i) {
        if (Par[i] != NULL) {
          delete[] Par[i];
        }
      }
      delete[] Par;
    }
    if(XPos != NULL) delete[] XPos;
    if(YResp != NULL) delete[] YResp;
    Par = NULL;
    XPos = YResp = NULL;
  }

Member Function Documentation

ConstructTSpline3()

TSpline3 * TSpline3_red::ConstructTSpline3 ( )

inline

CW: Make a TSpline3 from the reduced splines.

Definition at line 661 of file SplineStructs.h.

                                       {
    // KS: Sadly ROOT only accepts double...
    #ifdef _LOW_MEMORY_STRUCTS_
    std::vector<Double_t> xPosDoubles(nPoints);
    std::vector<Double_t> yPosDoubles(nPoints);
    for (Int_t i = 0; i < nPoints; ++i) {
      xPosDoubles[i] = static_cast<Double_t>(XPos[i]); // Convert float to double
      yPosDoubles[i] = static_cast<Double_t>(YResp[i]); // Convert float to double
    }
    TSpline3 *spline = new TSpline3("Spline", xPosDoubles.data(), yPosDoubles.data(), static_cast<int>(nPoints));
    #else
    TSpline3 *spline = new TSpline3("Spline", XPos, YResp, nPoints);
    #endif
    return spline;
  }

Eval()

double TSpline3_red::Eval ( double  var )

inlineoverridevirtual

CW: Evaluate the weight from a variation.

Implements TResponseFunction_red.

Definition at line 631 of file SplineStructs.h.

                                          {
    // Get the segment for this variation
    int segment = FindX(var);
    // The get the coefficients for this variation
    M3::float_t x = M3::float_t(-999.99), y = M3::float_t(-999.99), b = M3::float_t(-999.99), c = M3::float_t(-999.99), d = M3::float_t(-999.99);
    GetCoeff(segment, x, y, b, c, d);
    double dx = var - x;
    // Evaluate the third order polynomial
    double weight = y+dx*(b+dx*(c+d*dx));
    return weight;
  }

FindX()

int TSpline3_red::FindX ( double  x )

inline

Find the segment relevant to this variation in x See root/hist/hist/src/TSpline3FindX(double) or samplePDFND....::FindSplineSegment.

Definition at line 597 of file SplineStructs.h.

                             {
    // The segment we're interested in (klow in ROOT code)
    int segment = 0;
    int kHigh = nPoints-1;
    // If the variation is below the lowest saved spline point
    if (x <= XPos[0]){
      segment = 0;
      // If the variation is above the highest saved spline point
    } else if (x >= XPos[nPoints-1]) {
      // Yes, the -2 is indeed correct, see TSpline.cxx:814 and //see: https://savannah.cern.ch/bugs/?71651
      segment = kHigh;
      // If the variation is between the maximum and minimum, perform a binary search
    } else {
      // The top point we've got
      int kHalf = 0;
      // While there is still a difference in the points (we haven't yet found the segment)
      // This is a binary search, incrementing segment and decrementing kHalf until we've found the segment
      while (kHigh - segment > 1) {
        // Increment the half-step
        kHalf = (segment + kHigh)/2;
        // If our variation is above the kHalf, set the segment to kHalf
        if (x > XPos[kHalf]) {
          segment = kHalf;
          // Else move kHigh down
        } else {
          kHigh = kHalf;
        }
      } // End the while: we've now done our binary search
    } // End the else: we've now found our point
    if (segment >= nPoints-1 && nPoints > 1) segment = nPoints-2;
    return segment;
  }

GetCoeff()

void TSpline3_red::GetCoeff ( int  segment, M3::float_t &  x, M3::float_t &  y, M3::float_t &  b, M3::float_t &  c, M3::float_t &  d  )

inline

CW: Get the coefficient of a given segment.

Definition at line 652 of file SplineStructs.h.

                                                                                                          {
    b = Par[segment][0];
    c = Par[segment][1];
    d = Par[segment][2];
    x = XPos[segment];
    y = YResp[segment];
  }

GetKnot()

void TSpline3_red::GetKnot ( int  i, M3::float_t &  xtmp, M3::float_t &  ytmp  )

inline

Definition at line 646 of file SplineStructs.h.

                                                               {
    xtmp = XPos[i];
    ytmp = YResp[i];
  }

GetNp()

M3::int_t TSpline3_red::GetNp ( )

inlineoverridevirtual

CW: Get the number of points.

Implements TResponseFunction_red.

Definition at line 644 of file SplineStructs.h.

{ return nPoints; }

Print()

void TSpline3_red::Print ( )

inlineoverridevirtual

Print detailed info.

Implements TResponseFunction_red.

Definition at line 678 of file SplineStructs.h.

                               {
    MACH3LOG_INFO("Printing TSpline_red:");
    MACH3LOG_INFO(" Nknots = {}", nPoints);
    for (int i = 0; i < nPoints; ++i) {
      MACH3LOG_INFO("  i = {} x = {} y = {} b = {} c = {} d = {}",
                    i, XPos[i], YResp[i], Par[i][0], Par[i][1], Par[i][2]);
    }
  }

SetFunc()

void TSpline3_red::SetFunc ( TSpline3 *&  spline, SplineInterpolation  InterPolation = kTSpline3  )

inline

Set the function.

Definition at line 312 of file SplineStructs.h.

                                                                                        {
    nPoints = M3::int_t(spline->GetNp());
    if (Par != NULL) {
      for (int i = 0; i < nPoints; ++i) {
        delete[] Par[i];
        Par[i] = NULL;
      }
      delete[] Par;
      Par = NULL;
    }
    if (XPos != NULL) delete[] XPos;
    if (YResp != NULL) delete[] YResp;
    // Save the parameters for each knot
    Par = new M3::float_t*[nPoints];
    // Save the positions of the knots
    XPos = new M3::float_t[nPoints];
    // Save the y response at each knot
    YResp = new M3::float_t[nPoints];
 
    //KS: Default TSpline3 ROOT implementation
    if(InterPolation == kTSpline3)
    {
      for (int i = 0; i < nPoints; ++i) {
        // 3 is the size of the TSpline3 coefficients
        Par[i] = new M3::float_t[3];
        double x = -999.99, y = -999.99, b = -999.99, c = -999.99, d = -999.99;
        spline->GetCoeff(i, x, y, b, c, d);
        XPos[i]   = M3::float_t(x);
        YResp[i]  = M3::float_t(y);
        Par[i][0] = M3::float_t(b);
        Par[i][1] = M3::float_t(c);
        Par[i][2] = M3::float_t(d);
      }
    }
    //CW: Reduce to use linear spline interpolation for certain parameters
    // Not the most elegant way: use TSpline3 object but set coefficients to zero and recalculate spline points; the smart way (but more human intensive) would be to save memory here and simply not store the zeros at all
    // Get which parameters should be linear from the fit manager
    // Convert the spline number to global xsec parameter
    // Loop over the splines points
    // KS: kLinearFunc should be used with TF1, this is just as safety
    else if(InterPolation == kLinear || InterPolation == kLinearFunc)
    {
      for (int k = 0; k < nPoints; ++k) {
        // 3 is the size of the TSpline3 coefficients
        Par[k] = new M3::float_t[3];
        Double_t x1, y1, b1, c1, d1, x2, y2, b2, c2, d2 = 0;
        spline->GetCoeff(k, x1, y1, b1, c1, d1);
        spline->GetCoeff(k+1, x2, y2, b2, c2, d2);
        double tempb = (y2-y1)/(x2-x1);
 
        XPos[k]   = M3::float_t(x1);
        YResp[k]  = M3::float_t(y1);
        Par[k][0] = M3::float_t(tempb);
        Par[k][1] = M3::float_t(0);
        Par[k][2] = M3::float_t(0);
      }
    }
    //EM: Akima spline is similar to regular cubic spline but is allowed to be discontinuous in 2nd derivative and coefficients in any segment
    // only depend on th 2 nearest points on either side
    else if(InterPolation == kAkima)
    {
      // get the knot values for the spline
      for (int i = 0; i < nPoints; ++i) {
        // 3 is the size of the TSpline3 coefficients
        Par[i] = new M3::float_t[3];
 
        double x = -999.99, y = -999.99;
        spline->GetKnot(i, x, y);
 
        XPos[i]   = M3::float_t(x);
        YResp[i]  = M3::float_t(y);
      }
 
      M3::float_t* mvals = new M3::float_t[nPoints + 3];
      M3::float_t* svals = new M3::float_t[nPoints + 1];
 
      for (int i = -2; i <= nPoints; ++i) {
        // if segment is first or last or 2nd to first or last, needs to be dealt with slightly differently;
        // need to estimate the values for additinal points which would lie outside of the spline
        if(i ==-2){
          mvals[i+2] = M3::float_t(3.0 * (YResp[1] - YResp[0]) / (XPos[1] - XPos[0]) - 2.0*(YResp[2] - YResp[1]) / (XPos[2] - XPos[1]));
        }
        else if(i==-1){
          mvals[i+2] = M3::float_t(2.0 * (YResp[1] - YResp[0]) / (XPos[1] - XPos[0]) - (YResp[2] - YResp[1]) / (XPos[2] - XPos[1]));
        }
        else if(i==nPoints){
          mvals[i+2] = M3::float_t(3.0 * (YResp[nPoints-1] - YResp[nPoints-2]) / (XPos[nPoints-1] - XPos[nPoints-2]) - 2.0*(YResp[nPoints-2] - YResp[nPoints-3]) / (XPos[nPoints-2] - XPos[nPoints-3]));
        }
        else if(i == nPoints - 1){
          mvals[i+2] = M3::float_t(2.0 * (YResp[nPoints-1] - YResp[nPoints-2]) / (XPos[nPoints-1] - XPos[nPoints-2]) - (YResp[nPoints-2] - YResp[nPoints-3]) / (XPos[nPoints-2] - XPos[nPoints-3]));
        }
        //standard internal segment
        else{
          mvals[i+2] = (YResp[i+1] - YResp[i])/ (XPos[i+1] - XPos[i]);
        }
      }
 
      for(int i =2; i<=nPoints+2; i++){
        if (std::abs(mvals[i+1] - mvals[i]) + std::abs(mvals[i-1] - mvals[i-2]) != 0.0){
          svals[i-2] = (std::abs(mvals[i+1] - mvals[i]) * mvals[i-1] + std::abs(mvals[i-1] - mvals[i-2]) *mvals[i]) / (std::abs(mvals[i+1] - mvals[i]) + std::abs(mvals[i-1] - mvals[i-2]));
        }
        else{svals[i-2] = mvals[i];}
      }
 
      // calculate the coefficients for the spline
      for(int i = 0; i <nPoints; i++){
        M3::float_t b, c, d = M3::float_t(-999.999);
 
        b = svals[i];
        c = M3::float_t(3.0* (YResp[i+1] - YResp[i]) / (XPos[i+1] - XPos[i]) -2.0 *svals[i] - svals[i +1]) /(XPos[i+1] - XPos[i]);
        d = M3::float_t((svals[i + 1] +svals[i]) - 2.0*(YResp[i+1] - YResp[i]) / (XPos[i+1] - XPos[i])) / ((XPos[i+1] - XPos[i]) * (XPos[i+1] - XPos[i]));
 
        Par[i][0] = b;
        Par[i][1] = c;
        Par[i][2] = d;
      }
 
      // check the input spline for linear segments, if there are any then overwrite the calculated coefficients
      // this will pretty much only ever be the case if they are set to be linear in samplePDFND i.e. the user wants it to be linear
      for(int i = 0; i <nPoints-1; i++){
        double x = -999.99, y = -999.99, b = -999.99, c = -999.99, d = -999.99;
        spline->GetCoeff(i, x, y, b, c, d);
 
        if((c == 0.0 && d == 0.0)){
          Par[i][0] = M3::float_t(b);
          Par[i][1] = M3::float_t(0.0);
          Par[i][2] = M3::float_t(0.0);
        }
      }
      delete[] mvals;
      delete[] svals;
    }
    //EM: Monotone spline is similar to regular cubic spline but enforce the condition that the interpolated value at any point
    // must be between its two nearest knots, DOES NOT make the entire spline monotonic, only the segments
    else if(InterPolation == kMonotonic)
    {
      // values of the secants at each point (for calculating monotone spline)
      M3::float_t * Secants = new M3::float_t[nPoints -1];
      // values of the tangens at each point (for calculating monotone spline)
      M3::float_t *  Tangents = new M3::float_t[nPoints];
 
      // get the knot values for the spline
      for (int i = 0; i < nPoints; ++i) {
        // 3 is the size of the TSpline3 coefficients
        Par[i] = new M3::float_t[3];
 
        double x = -999.99, y = -999.99;
        spline->GetKnot(i, x, y);
 
        XPos[i]   = M3::float_t(x);
        YResp[i]  = M3::float_t(y);
 
        Tangents[i] = 0.0;
      }
 
      // deal with the case of two points (just do linear interpolation between them)
      if (nPoints == 2){
        Par[0][0] = (YResp[1] - YResp[0]) / (XPos[1] - XPos[0]);
        Par[0][1] = 0.0;
        Par[0][2] = 0.0;
        // extra "virtual" segment at end to make Par array shape fit with knot arrays shapes
        Par[1][1] = 0.0;
        Par[1][2] = 0.0;
 
        return;
      } // if nPoints !=2 do full monotonic spline treatment:
      else
      {
        // first pass over knots to calculate the secants
        for (int i = 0; i < nPoints-1; ++i) {
          Secants[i] = (YResp[i+1] - YResp[i]) / (XPos[i+1] - XPos[i]);
          MACH3LOG_TRACE("Secant {}: {}", i, Secants[i]);
        }
 
        Tangents[0] = Secants[0];
        Tangents[nPoints-1] = Secants[nPoints -2];
 
        M3::float_t alpha;
        M3::float_t beta;
 
        // second pass over knots to calculate tangents
        for (int i = 1; i < nPoints-1; ++i) {
          if ((Secants[i-1] >= 0.0 && Secants[i] >= 0.0) | (Secants[i-1] < 0.0 && Secants[i] < 0.0)){ //check for same sign
            Tangents[i] = M3::float_t((Secants[i-1] + Secants[i]) /2.0);
          }
        }
 
        // third pass over knots to rescale tangents
        for (int i = 0; i < nPoints-1; ++i) {
          if (Secants[i] == 0.0){
            Tangents[i] = 0.0;
            Tangents[i+1] = 0.0;
          }
 
          else{
            alpha = Tangents[i]  / Secants[i];
            beta = Tangents[i+1] / Secants[i];
 
            if (alpha <0.0){
              Tangents[i] = 0.0;
            }
            if (beta < 0.0){
              Tangents[i+1] = 0.0;
            }
 
            if (alpha * alpha + beta * beta >9.0){
              M3::float_t tau = M3::float_t(3.0 / std::sqrt(alpha * alpha + beta * beta));
              Tangents[i]   = tau * alpha * Secants[i];
              Tangents[i+1] = tau * beta  * Secants[i];
            }
          }
        } // finished rescaling tangents
        // fourth pass over knots to calculate the coefficients for the spline
        M3::float_t dx;
        for(int i = 0; i <nPoints-1; i++){
          M3::float_t b, c, d = M3::float_t(-999.999);
          dx = XPos[i+1] - XPos[i];
 
          b = Tangents[i] * dx;
          c = M3::float_t(3.0* (YResp[i+1] - YResp[i]) -2.0 *dx * Tangents[i] - dx * Tangents[i +1]);
          d = M3::float_t(2.0* (YResp[i] - YResp[i+1]) + dx * (Tangents[i] + Tangents[i+1]));
 
          Par[i][0] = b /  dx;
          Par[i][1] = c / (dx * dx);
          Par[i][2] = d / (dx * dx * dx);
 
          if((Par[i][0] == -999) | (Par[i][1] == -999) | (Par[i][2] ==-999) | (Par[i][0] == -999.999) | (Par[i][1] == -999.999) | (Par[i][2] ==-999.999)){
            MACH3LOG_INFO("Bad spline parameters for segment {}: (b, c, d) = {}, {}, {}. This will cause problems with GPU.",
                          i, Par[i][0], Par[i][1], Par[i][2]);
          }
          MACH3LOG_TRACE("b: {}", b);
          MACH3LOG_TRACE("dx: {}, x_0: {}, x_1: {}", dx, XPos[i], XPos[i+1]);
          MACH3LOG_TRACE("    y_0: {}, y_1: {}", YResp[i], YResp[i+1]);
        }
 
        // include params for final "segment" outside of the spline so that par array fits with x and y arrays,
        // should never actually get used but if not set then the GPU code gets very angry
        Par[nPoints-1][0] = 0.0;
        Par[nPoints-1][1] = 0.0;
        Par[nPoints-1][2] = 0.0;
 
        // check the input spline for linear segments, if there are any then overwrite the calculated coefficients
        // this will pretty much only ever be the case if they are set to be linear in samplePDFND i.e. the user wants it to be linear
        for(int i = 0; i <nPoints-1; i++){
          double x = -999.99, y = -999.99, b = -999.99, c = -999.99, d = -999.99;
          spline->GetCoeff(i, x, y, b, c, d);
 
          if((c == 0.0 && d == 0.0)){
            Par[i][0] = M3::float_t(b);
            Par[i][1] = 0.0;
            Par[i][2] = 0.0;
          }
        }
        delete[] Secants;
        delete[] Tangents;
      } // end of if(nPoints !=2)
    }
    else
    {
      MACH3LOG_ERROR("Unsupported interpolation type {}", static_cast<int>(InterPolation));
      throw MaCh3Exception(__FILE__ , __LINE__ );
    }
 
    delete spline;
    spline = NULL;
  }

Member Data Documentation

nPoints

M3::int_t TSpline3_red::nPoints

protected

Number of points/knot in TSpline3.

Definition at line 689 of file SplineStructs.h.

Par

M3::float_t** TSpline3_red::Par

protected

Always uses a third order polynomial, so hard-code the number of coefficients in implementation.

Definition at line 691 of file SplineStructs.h.

XPos

M3::float_t* TSpline3_red::XPos

protected

Positions of each x for each knot.

Definition at line 693 of file SplineStructs.h.

YResp

M3::float_t* TSpline3_red::YResp

protected

y-value for each knot

Definition at line 695 of file SplineStructs.h.

---

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

• Splines/SplineStructs.h