SplineStructs.h

Go to the documentation of this file.

#pragma once
 
// MaCh3  includes
#include "Parameters/ParameterHandlerGeneric.h"
#include "Parameters/ParameterStructs.h"
 
#include <cmath>
 
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wuseless-cast"
#pragma GCC diagnostic ignored "-Wfloat-conversion"
#pragma GCC diagnostic ignored "-Wnull-dereference"
 
 
// *******************
struct FastSplineInfo {
// *******************
  FastSplineInfo() {
    nPts = -999;
    CurrSegment = 0;
    splineParsPointer = nullptr;
  }
 
  virtual ~FastSplineInfo() {
 
  }
 
  M3::int_t nPts;
 
  std::vector<M3::float_t> xPts;
 
  M3::int_t CurrSegment;
 
  const double* splineParsPointer;
};
 
 
// ***************************************************************************
inline void ApplyKnotWeightCap(TGraph* xsecgraph, const int splineParsIndex, ParameterHandlerGeneric* XsecCov) {
// ***************************************************************************
  if(xsecgraph == nullptr){
    MACH3LOG_ERROR("hmmm looks like you're trying to apply capping for spline parameter {} but it hasn't been set in xsecgraph yet", XsecCov->GetParFancyName(splineParsIndex));
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  // EM: cap the weights of the knots if specified in the config
  if(
    (XsecCov->GetParSplineKnotUpperBound(splineParsIndex) != M3::DefSplineKnotUpBound)
    || (XsecCov->GetParSplineKnotLowerBound(splineParsIndex) != M3::DefSplineKnotLowBound))
  {
    for(int knotId = 0; knotId < xsecgraph->GetN(); knotId++){
      double x,y;
 
      // EM: get the x and y of the point. Also double check that the requested point was valid just to be super safe
      if( xsecgraph->GetPoint(knotId, x, y) == -1) {
        MACH3LOG_ERROR("Invalid knot requested: {}", knotId);
        throw MaCh3Exception(__FILE__ , __LINE__ );
      }
 
      y = std::min(y, XsecCov->GetParSplineKnotUpperBound(splineParsIndex));
      y = std::max(y, XsecCov->GetParSplineKnotLowerBound(splineParsIndex));
 
      xsecgraph->SetPoint(knotId, x, y);
    }
 
    // EM: check if our cap made the spline flat, if so we set to 1 knot to avoid problems later on
    bool isFlat = true;
    for(int knotId = 0; knotId < xsecgraph->GetN(); knotId++){
      double x,y;
 
      // EM: get the x and y of the point. Also double check that the requested point was valid just to be super safe
      if( xsecgraph->GetPoint(knotId, x, y) == -1) {
        MACH3LOG_ERROR("Invalid knot requested: {}", knotId);
        throw MaCh3Exception(__FILE__ , __LINE__ );
      }
      if(std::abs(y - 1.0) > 1e-5) isFlat = false;
    }
 
    if(isFlat){
      xsecgraph->Set(1);
      xsecgraph->SetPoint(0, 0.0, 1.0);
    }
  }
}
 
// ***************************************************************************
inline void ApplyKnotWeightCapTSpline3(TSpline3* &Spline, const int splineParsIndex, ParameterHandlerGeneric* XsecCov) {
// ***************************************************************************
  if(Spline == nullptr) {
    MACH3LOG_ERROR("hmmm looks like you're trying to apply capping for spline parameter {} but it hasn't been set in Spline yet", XsecCov->GetParFancyName(splineParsIndex));
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  std::string oldName = Spline->GetName();
  // EM: cap the weights of the knots if specified in the config
  if((XsecCov->GetParSplineKnotUpperBound(splineParsIndex) != M3::DefSplineKnotUpBound) || (XsecCov->GetParSplineKnotLowerBound(splineParsIndex) != M3::DefSplineKnotLowBound))
  {
    const int NValues = Spline->GetNp();
    std::vector<double> XVals(NValues);
    std::vector<double> YVals(NValues);
    for(int knotId = 0; knotId < NValues; knotId++){
      double x,y;
      // Extract X and Y
      Spline->GetKnot(knotId, x, y);
 
      y = std::min(y, XsecCov->GetParSplineKnotUpperBound(splineParsIndex));
      y = std::max(y, XsecCov->GetParSplineKnotLowerBound(splineParsIndex));
 
      XVals[knotId] = x;
      YVals[knotId] = y;
    }
    delete Spline;
    // If we capped we have to make new TSpline3 to recalculate coefficients
    Spline = new TSpline3(oldName.c_str(), XVals.data(), YVals.data(), NValues);
  }
}
 
// ************************
class TResponseFunction_red {
// ************************
public:
  TResponseFunction_red() { }
  virtual ~TResponseFunction_red() { }
  virtual double Eval(const double var)=0;
  virtual void Print()=0;
  virtual M3::int_t GetNp()=0;
};
 
// ************************
class TF1_red: public TResponseFunction_red {
// ************************
public:
  TF1_red() : TResponseFunction_red() {
    length = 0;
    Par = nullptr;
  }
 
  virtual ~TF1_red() {
    if (Par != NULL) {
      delete[] Par;
      Par = nullptr;
    }
  }
 
  TF1_red(M3::int_t nSize, M3::float_t* Array) : TResponseFunction_red() {
    length = nSize;
    for (int i = 0; i < length; ++i) {
      Par[i] = Array[i];
    }
  }
 
  TF1_red(TF1* &Function) : TResponseFunction_red() {
    Par = nullptr;
    SetFunc(Function);
  }
 
  inline void SetFunc(TF1* &Func) {
    length = M3::int_t(Func->GetNpar());
    if (Par != nullptr) delete[] Par;
    Par = new M3::float_t[length];
    for (int i = 0; i < length; ++i) {
      Par[i] = M3::float_t(Func->GetParameter(i));
    }
    delete Func;
    Func = nullptr;
  }
 
  inline double Eval(const double var) override {
    return Par[1]+Par[0]*var;
 
    /* FIXME in future we might introduce more TF1
    //If we have 5 parameters we're using a fifth order polynomial
    if (Type == kFifthOrderPolynomial) {
      return 1+Par[0]*var+Par[1]*var*var+Par[2]*var*var*var+Par[3]*var*var*var*var+Par[4]*var*var*var*var*var;
    } else if (Type == kTwoLinears) {
      return (var <= 0)*(Par[2]+Par[0]*var)+(var > 0)*(Par[2]+Par[1]*var);
    } else if (Type == kLinear) {
      return (Par[1]+Par[0]*var);
    } else if (Type == kPseudoHeaviside) {
      return (var <= 0)*(1+Par[0]*var) + (1 >= var)*(var > 0)*(1+Par[1]*var) + (var > 1)*(Par[3]+Par[2]*var);
    }else {
      MACH3LOG_ERROR("    Class only knows about 5th order polynomial, two superposed linear functions, linear function, or pseudo Heaviside.");
      MACH3LOG_ERROR("    You have tried something else than this, which remains unimplemented.");
      MACH3LOG_ERROR("{}: {}", __FILE__, __LINE__);
      throw MaCh3Exception(__FILE__ , __LINE__ );
    }
    */
  }
 
  inline void SetParameter(M3::int_t Parameter, M3::float_t Value) {
    Par[Parameter] = Value;
  }
 
  double GetParameter(M3::int_t Parameter) {
    if (Parameter > length) {
      MACH3LOG_ERROR("You requested parameter number {} but length is {} parameters", Parameter, length);
      throw MaCh3Exception(__FILE__ , __LINE__ );
      return -999.999;
    }
    return Par[Parameter];
  }
 
  inline void SetSize(M3::int_t nSpline) {
    length = nSpline;
    Par = new M3::float_t[length];
  }
  inline int GetSize() { return length; }
  inline void Print() override {
    MACH3LOG_INFO("Printing TF1_red:");
    MACH3LOG_INFO("  Length  = {}", length);
    for (int i = 0; i < length; i++) {
      MACH3LOG_INFO("  Coeff {}  = {}", i, Par[i]);
    }
  }
 
  inline TF1* ConstructTF1(const std::string& function, const int xmin, const int xmax) {
    TF1 *func = new TF1("TF1", function.c_str(), xmin, xmax);
    for(int i = 0; i < length; ++i) {
      func->SetParameter(i, Par[i]);
    }
    return func;
  }
 
  inline M3::int_t GetNp() override { return length; }
 
private:
  M3::float_t* Par;
  M3::int_t length;
};
 
// ************************
class TSpline3_red: public TResponseFunction_red {
// ************************
public:
  TSpline3_red() : TResponseFunction_red() {
    nPoints = 0;
    Par = nullptr;
    XPos = nullptr;
    YResp = nullptr;
  }
 
  TSpline3_red(TSpline3* &spline, SplineInterpolation InterPolation = kTSpline3) : TResponseFunction_red() {
    Par = nullptr;
    XPos = nullptr;
    YResp = nullptr;
    SetFunc(spline, InterPolation);
  }
 
  TSpline3_red(M3::float_t *X, M3::float_t *Y, M3::int_t N, M3::float_t **P) : 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("*********************************************************************************************");
      }
    }
  }
  inline void SetFunc(TSpline3* &spline, SplineInterpolation InterPolation = kTSpline3) {
    nPoints = M3::int_t(spline->GetNp());
    if (Par != nullptr) {
      for (int i = 0; i < nPoints; ++i) {
        delete[] Par[i];
        Par[i] = nullptr;
      }
      delete[] Par;
      Par = nullptr;
    }
    if (XPos != nullptr) delete[] XPos;
    if (YResp != nullptr) 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) {
        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);
 
        double tempb = 0;
        if (k == nPoints - 1) {
          tempb = Par[k-1][0];
        } else {
          spline->GetCoeff(k + 1, x2, y2, b2, c2, d2);
          tempb = (y2-y1)/(x2-x1);
        }
        XPos[k]   = M3::float_t(x1);
        YResp[k]  = M3::float_t(y1);
        Par[k][0] = M3::float_t(tempb);  // linear slope
        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 if (InterPolation == kKochanekBartels)
    {
      // Allocate memory for tangents and coefficients
      M3::float_t* Tangents = new M3::float_t[nPoints];
      for (int i = 0; i < nPoints; ++i)
      {
        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;
      }
 
      // KS: Setting all parameters to 0 gives a Catmull-Rom spline.
      // Tension (T):
      //   - T =  0.0: Default smooth interpolation (Catmull-Rom).
      //   - T =  1.0: Maximum tension; curve becomes linear between knots (no curvature).
      //   - T = -1.0: Overshooting; creates loops or "bouncy" effects around knots.
      constexpr M3::float_t T = 0.0;
 
      // Continuity (C):
      //   - C =  0.0: Default smooth transition (continuous first derivative).
      //   - C =  1.0: Sharp corner/crease at knot (discontinuous first derivative).
      //   - C = -1.0: Inverted curve; creates loops or kinks at knot.
      constexpr M3::float_t C = 0.0;
 
      // Bias (B):
      //   - B =  0.0: Default symmetric curve around knot.
      //   - B =  1.0: Curve is "pulled" toward the next knot (leading effect).
      //   - B = -1.0: Curve is "pulled" toward the previous knot (lagging effect).
      constexpr M3::float_t B = 0.0;
 
      // Calculate tangents for internal knots
      for (int i = 1; i < nPoints - 1; ++i)
      {
        M3::float_t d0 = (YResp[i]   - YResp[i-1]) / (XPos[i]   - XPos[i-1]);
        M3::float_t d1 = (YResp[i+1] - YResp[i])   / (XPos[i+1] - XPos[i]);
 
        M3::float_t term1 = (1.0 - T) * (1.0 + B) * (1.0 + C) * 0.5 * d0;
        M3::float_t term2 = (1.0 - T) * (1.0 - B) * (1.0 - C) * 0.5 * d1;
 
        Tangents[i] = term1 + term2;
      }
 
      // Boundary conditions (simple choice: secant slopes)
      Tangents[0]        = (YResp[1]       - YResp[0])       / (XPos[1]       - XPos[0]);
      Tangents[nPoints-1]= (YResp[nPoints-1] - YResp[nPoints-2]) / (XPos[nPoints-1] - XPos[nPoints-2]);
 
      // Compute cubic coefficients for each segment
      for (int i = 0; i < nPoints - 1; ++i)
      {
        M3::float_t dx = XPos[i+1] - XPos[i];
        M3::float_t dy = YResp[i+1] - YResp[i];
        M3::float_t m0 = Tangents[i];
        M3::float_t m1 = Tangents[i+1];
 
        Par[i][0] = m0;                                  // b
        Par[i][1] = (3*dy/(dx*dx)) - (2*m0 + m1)/dx;     // c
        Par[i][2] = (m0 + m1 - 2*dy/dx) / (dx*dx);       // d
 
        MACH3LOG_TRACE("KB segment {}: dx={}, dy={}, b={}, c={}, d={}",
                       i, dx, dy, Par[i][0], Par[i][1], Par[i][2]);
      }
      // Last "virtual" segment
      Par[nPoints-1][0] = 0.0;
      Par[nPoints-1][1] = 0.0;
      Par[nPoints-1][2] = 0.0;
 
      delete[] Tangents;
    }
    else
    {
      MACH3LOG_ERROR("Unsupported interpolation type {}", static_cast<int>(InterPolation));
      throw MaCh3Exception(__FILE__ , __LINE__ );
    }
 
    delete spline;
    spline = nullptr;
  }
 
  virtual ~TSpline3_red() {
    if(Par != nullptr) {
      for (int i = 0; i < nPoints; ++i) {
        if (Par[i] != nullptr) {
          delete[] Par[i];
        }
      }
      delete[] Par;
    }
    if(XPos != nullptr) delete[] XPos;
    if(YResp != nullptr) delete[] YResp;
    Par = nullptr;
    XPos = YResp = nullptr;
  }
 
  inline int FindX(double x) {
    // 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;
  }
 
  inline double Eval(double var) override {
    // 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;
  }
 
  inline M3::int_t GetNp() override { return nPoints; }
  // Get the ith knot's x and y position
  inline void GetKnot(int i, M3::float_t &xtmp, M3::float_t &ytmp) {
    xtmp = XPos[i];
    ytmp = YResp[i];
  }
 
  inline void GetCoeff(int segment, M3::float_t &x, M3::float_t &y, M3::float_t &b, M3::float_t &c, M3::float_t &d) {
    b = Par[segment][0];
    c = Par[segment][1];
    d = Par[segment][2];
    x = XPos[segment];
    y = YResp[segment];
  }
 
  inline TSpline3* ConstructTSpline3() {
    // 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
    for (Int_t i = 0; i < nPoints; ++i) {
      spline->SetPointCoeff(i, Par[i][0], Par[i][1], Par[i][2]);
    }
 
    return spline;
  }
 
  inline void Print() override {
    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]);
    }
  }
 
  protected: //changed to protected from private so can be accessed by derived classes
    M3::int_t nPoints;
    M3::float_t **Par;
    M3::float_t *XPos;
    M3::float_t *YResp;
};
 
// *****************************************
inline bool isFlat(TSpline3_red* &spl) {
// *****************************************
  int Np = spl->GetNp();
  M3::float_t x, y, b, c, d;
  // Go through spline segment parameters,
  // Get y values for each spline knot,
  // Every knot must evaluate to 1.0 to create a flat spline
  for(int i = 0; i < Np; i++) {
    spl->GetCoeff(i, x, y, b, c, d);
    if (y != 1) {
      return false;
    }
  }
  return true;
}
 
// *********************************
inline std::vector<std::vector<TSpline3_red*> > ReduceTSpline3(std::vector<std::vector<TSpline3*> > &MasterSpline) {
// *********************************
  std::vector<std::vector<TSpline3*> >::iterator OuterIt;
  std::vector<TSpline3*>::iterator InnerIt;
 
  // The return vector
  std::vector<std::vector<TSpline3_red*> > ReducedVector;
  ReducedVector.reserve(MasterSpline.size());
 
  // Loop over each parameter
  int OuterCounter = 0;
  for (OuterIt = MasterSpline.begin(); OuterIt != MasterSpline.end(); ++OuterIt, ++OuterCounter) {
    // Make the temp vector
    std::vector<TSpline3_red*> TempVector;
    TempVector.reserve(OuterIt->size());
    int InnerCounter = 0;
    // Loop over each TSpline3 pointer
    for (InnerIt = OuterIt->begin(); InnerIt != OuterIt->end(); ++InnerIt, ++InnerCounter) {
      // Here's our delicious TSpline3 object
      TSpline3 *spline = (*InnerIt);
      // Now make the reduced TSpline3 pointer
      TSpline3_red *red = nullptr;
      if (spline != NULL) {
        red = new TSpline3_red(spline);
        (*InnerIt) = spline;
      }
      // Push back onto new vector
      TempVector.push_back(red);
    } // End inner for loop
    ReducedVector.push_back(TempVector);
  } // End outer for loop
  // Now have the reduced vector
  return ReducedVector;
}
 
// *********************************
inline std::vector<std::vector<TF1_red*> > ReduceTF1(std::vector<std::vector<TF1*> > &MasterSpline) {
// *********************************
  std::vector<std::vector<TF1*> >::iterator OuterIt;
  std::vector<TF1*>::iterator InnerIt;
 
  // The return vector
  std::vector<std::vector<TF1_red*> > ReducedVector;
  ReducedVector.reserve(MasterSpline.size());
 
  // Loop over each parameter
  int OuterCounter = 0;
  for (OuterIt = MasterSpline.begin(); OuterIt != MasterSpline.end(); ++OuterIt, ++OuterCounter) {
    // Make the temp vector
    std::vector<TF1_red*> TempVector;
    TempVector.reserve(OuterIt->size());
    int InnerCounter = 0;
    // Loop over each TSpline3 pointer
    for (InnerIt = OuterIt->begin(); InnerIt != OuterIt->end(); ++InnerIt, ++InnerCounter) {
      // Here's our delicious TSpline3 object
      TF1* spline = (*InnerIt);
      // Now make the reduced TSpline3 pointer (which deleted TSpline3)
      TF1_red* red = nullptr;
      if (spline != NULL) {
        red = new TF1_red(spline);
        (*InnerIt) = spline;
      }
      // Push back onto new vector
      TempVector.push_back(red);
    } // End inner for loop
    ReducedVector.push_back(TempVector);
  } // End outer for loop
  // Now have the reduced vector
  return ReducedVector;
}
 
// *********************************
inline TResponseFunction_red* CreateResponseFunction(TGraph* &graph,
                                                     const RespFuncType SplineRespFuncType,
                                                     const SplineInterpolation SplineInterpolationType,
                                                     const std::string& Title) {
// *********************************
  TResponseFunction_red* RespFunc = nullptr;
 
  if (graph && graph->GetN() > 1)
  {
    if(SplineRespFuncType == kTSpline3_red)
    {
      // Here's the TSpline3
      TSpline3* spline = nullptr;
      TSpline3_red *spline_red = nullptr;
 
      // Create the TSpline3* from the TGraph* and build the coefficients
      spline = new TSpline3(Title.c_str(), graph);
      spline->SetNameTitle(Title.c_str(), Title.c_str());
 
      // Make the reduced TSpline3 format and delete the old spline
      spline_red = new TSpline3_red(spline, SplineInterpolationType);
 
      RespFunc = spline_red;
    }
    else if(SplineRespFuncType == kTF1_red)
    {
      // The TF1 object we build from fitting the TGraph
      TF1 *Fitter = nullptr;
      TF1_red *tf1_red = nullptr;
 
      if(graph->GetN() != 2) {
        MACH3LOG_ERROR("Trying to make TF1 from more than 2 knots.  Knots = {}", graph->GetN());
        MACH3LOG_ERROR("Currently support only linear with 2 knots :(");
        throw MaCh3Exception(__FILE__ , __LINE__ );
      }
 
      // Try simple linear function
      Fitter = new TF1(Title.c_str(), "([1]+[0]*x)", graph->GetX()[0], graph->GetX()[graph->GetN()-1]);
      //CW: For 2p2h shape C and O we can't fit a polynomial: try a linear combination of two linear functions around 0
      //Fitter = new TF1(Title.c_str(), "(x<=0)*(1+[0]*x)+(x>0)*([1]*x+1)", graph->GetX()[0], graph->GetX()[graph->GetN()-1]);
      // Fit 5hd order polynomial for all other parameters
      //Fitter = new TF1(Title.c_str(), "1+[0]*x+[1]*x*x+[2]*x*x*x+[3]*x*x*x*x+[4]*x*x*x*x*x", graph->GetX()[0], graph->GetX()[graph->GetN()-1]);
      //Pseudo Heaviside for Pauli Blocking
      //Fitter = new TF1(Title.c_str(), "(x <= 0)*(1+[0]*x) + (1 >= x)*(x > 0)*(1+[1]*x) + (x > 1)*([3]+[2]*x)", graph->GetX()[0], graph->GetX()[graph->GetN()-1]);
 
      // Fit the TF1 to the graph
      graph->Fit(Fitter, "Q0");
      // Make the reduced TF1 if we want
      tf1_red = new TF1_red(Fitter);
 
      RespFunc = tf1_red;
    }
    else
    {
      MACH3LOG_ERROR("Unsupported response function type");
      throw MaCh3Exception(__FILE__ , __LINE__ );
    }
  }
  else
  {
    RespFunc = nullptr;
  }
  return RespFunc;
}
 
#pragma GCC diagnostic pop