#ifndef MySpline2_C
#define MySpline2_C

#include "Geom.C"

#include <TF1.h>
#include <TString.h>
#include <TSystemDirectory.h>
#include <TList.h>
#include <TMath.h>

#include <iostream>
#include <vector>
#include <fstream>
#include <string>
#include <cassert>

namespace MySpline2
{
  class MySpline2
  {
  public:
    
    // The following function implements recursively
    // a Horner's-scheme-evaluated polynomial.
    // NOTE: The N parameter here is the NUMBER OF TERMS in the
    // polynomial, it is NOT the order.  A quadratic has N = 3.
    // The higher-order term is c[N], the constant term is c[0].
    static double Poly(unsigned int N, const double x, const double * const c)
    {
      if( N > 1 )
        {
          return  c[0] + x*Poly(N-1,x,c+1);
        }
      else if( N == 1 )
        {
          return c[0];
        }
      else
        {
          return 0;
        }
    }

    static double Poly(unsigned int N, const double * const x, const double * const c)
    {
      return Poly(N,x[0],c);
    }
    
    // This is a function for the nth derivative of an N-Polynomial
    // as evaluated by the other functions above.  Remember that N is the
    // number of coefficients in the polynomial, NOT the order.  
    static double PolyDeriv(unsigned int N, unsigned int n, 
                            const double x, const double * const c)
    { // nth derivative of Nth order polynomial.
      // f_N^{(n)} = \sum_{j = 0}^{N-n} (j+n)!/j! c_{j+n} x^j
      // The result is just an N-mth order polynomial with transformed 
      // coefficients.  The jth order coefficient is (j+n)!/j!*c_{j+n}.
      // TODO: Find a way to evaluate the coefficient multipliers (the
      // factorials) at compile-time.
      if(n == 0) return Poly(N,x,c);
      //double cprime[N-n];
      std::vector<double> cprime(N-n,0); // Fill with N-n zeros.
      for(unsigned int j = 0; j < N-n; j++)
        {
          cprime[j] = c[j+n];
          for(unsigned int k = j+n; k > j; k--)
            {
              cprime[j] *= k;
            }
        }
      //return Poly(N-n,x,cprime);
      return Poly(N-n,x,&(cprime[0]));
    }

    static double PolyDeriv(unsigned int N, unsigned int n,
                            const double * const x, const double * const c)
    {
      return PolyDeriv(N,n,x[0],c);
    }
  };



  // This class represents a piecewise continuous or smooth polynomial
  // with a fixed number & position of junctions (knots).  The number
  // of coefficients of the segments N is also fixed and uniform.  The
  // knot x values are fixed, but not the y values.
  // When a PieceWise object is instantiated, you can ask its
  // n_c_total member how many parameters it needs.  The first N
  // parameters are the coefficients of the first segment.  Subsequent
  // parameters are the high-order coefficients of the other segments.
  // Continuity (smoothness order of 0) is enforced.
  class PieceWise
  {
  public:
    PieceWise() : N_knots(0), N(0), smoothness(0), n_c_total(0), n_c_seg(0) {};
    PieceWise(unsigned int N_knots_, unsigned int N_, unsigned int smoothness_,
              const double * const knots_) : 
      N_knots(N_knots_), N(N_), smoothness(smoothness_),
      n_c_total((N-(smoothness+1))*N_knots + N), n_c_seg(N-(smoothness+1)),
      knots_x(knots_, knots_+N_knots) { };
    //PieceWise(TString fname);
    
    double operator()(const double * const x, const double * const c) const 
    { return Deriv(0,x[0],c); };
    double operator()(const double * const x) const { return Deriv(0,x[0],coeffs.data()); };

    double Eval(const double x, const double * const c) const { return Deriv(0,x,c); };
    double Eval(const double x) const { return Deriv(0,x,coeffs.data()); };
    
    double Deriv(unsigned int d, const double x) const { return Deriv(d,x,coeffs.data()); };
    double Deriv(unsigned int d, const double x, const double * c) const
    { // Number of coefficients in c is n_c_total.

      // The first N are used for the first segment, right here:
      if(N_knots < 1 or x <= knots_x[0]) return MySpline2::PolyDeriv(N,d,x,c);
      
      // The following blocks of N-(s+1) elements of c are for each
      // following segment.  They are doing to be the HIGH-ORDER
      // elements of the segments, because the first s coefficients of the
      // later segments are determined by previous segments.

      // This algorithm finds the index i of the segment that we are on.
      // So if i == 3, then knots_x[i-1] < x <= knots_x[i]
      // if i == N_knots then x is past the last knot.
      const auto i = std::lower_bound(knots_x.begin(),
                                      knots_x.begin()+N_knots,x) - knots_x.begin();
      // x_knot is the x value at the knot on the left of this segment.
      double x_knot = (knots_x[i-1]); 
      
      // We need to construct an array of coefficients for this segment.
      double * cprime = new double[N];
      
      // We are in segment i > 0.  We have to copy over the
      // parts of "c" that we care about.  We first skip N
      // elements from the first segment, then n_c_seg*i-1 for
      // the other segments before this one.  We copy n_c_seg
      // coefficients into the high-order part of cprime.
      std::copy(c+N+n_c_seg*(i-1), c+N+n_c_seg*i, cprime+smoothness+1);
      // Fill the rest with zeros.
      std::fill(cprime,cprime+smoothness+1,0);

      // Now compute each of the other cprimes based on the
      // derivative values of the previous segment at the knot.
      // The weird loop is because it's unsigned and I want to iterate
      // over [smoothness..0] inclusive.  Uses the famous "goes to" operator.
      for(unsigned int s = smoothness+1; s --> 0 ;) // "s goes to 0"
        {
          cprime[s] = this->Deriv(s,x_knot,c) - 
            MySpline2::PolyDeriv(N,s,x_knot,cprime);
        }

      auto value = MySpline2::PolyDeriv(N,d,x,cprime);
      delete[] cprime;
      return value;
    };

    void SetCoeffs(const double * const newcoeffs)
    {
      //std::copy(newcoeffs,newcoeffs+n_c_total,coeffs);
      coeffs.clear();
      coeffs.assign(newcoeffs,newcoeffs+n_c_total);
    };
    
    unsigned int N_knots;
    unsigned int N;
    unsigned int smoothness;

    // The total number of coefficients is given by 
    // N + (N - (smoothness+1) )*N_knots.
    // The first N is for the first (unconstrained) segment.  Then
    // there are N-(s+1) coefficients per knot.  It's s+1 because s=0
    // still imposes a constraint: continuity.
    unsigned int n_c_total; 
    
    // The number of coefficients per segment n_c_seg is given by
    // N-(smoothness+1)
    unsigned int n_c_seg;
    std::vector<double> knots_x;
    std::vector<double> coeffs;
  }; // End of class PieceWise
  

  class Poly
  { // This class allows for efficient evaluation of mostly-fixed coefficient
    // polynomials.  It calculates the derived coefficients for derivatives upon
    // creation of updating of the coefficients and stores them in an internal vector.
  public:
    Poly(unsigned int N_, const double * const c_) : N(N_)
    { // N_: number of coefficients.
      // c_: array of N_ coefficients.  c_[0] is the constant term, 
      //     c_[N_-1] is the high-order term.
      SetCoeffs(N,c_);
    };
    Poly() : N(0) {};

    // Evaluate the polynomial at position x using the current coefficients.
    double Eval(double x) const { return Eval(N,x,coeffs[0].data()); };
    // Evaluate the nth derivative of the polynomial at the position x.
    double Deriv(unsigned int n, double x) const { return Eval(N-n,x,coeffs[n].data()); };


    void SetCoeffs(const double * const c)
    { // Set new coefficients for the polynomial.
      // This SetCoeffs does not change the number of coefficients.
      // Assumes enough coefficients are passed for the current N.
      coeffs[0].assign(c,c+N);
      for(unsigned int n = 1; n < N; n++)
      {
        for(unsigned int j = 0; j < N-n;j++)
        {
          coeffs[n][j] = c[j+n];
          for(unsigned int k = j+n; k > j; k--)
          {
            coeffs[n][j] *= k;
          }
        }
      }
    };

    void SetCoeffs(unsigned int N_, const double * const c)
    { // Set new coefficients for the polynomial, and change the number.
      // This overload is for when you want to change the number
      // of coefficients.  It's a bit faster than making a whole new Poly
      // object because it re-uses the coeff vector.
      N = N_;
      coeffs.clear();
      coeffs.emplace_back(c,c+N);
      for(unsigned int n = 1; n < N; n++)
      {
        std::vector<double> cprime(N-n,0);
        for(unsigned int j = 0; j < N-n;j++)
        {
          cprime[j] = c[j+n];
          for(unsigned int k = j+n; k > j; k--)
          {
            cprime[j] *= k;
          }
        }
        coeffs.push_back(std::move(cprime));
      }
    };

    // Retrieve a const pointer to the current coefficients.
    double const * GetCoeffs() const { return coeffs[0].data(); };

    // Print all the current coefficients & derived ones.
    void PrintCoeffs() const
    {
      for(unsigned int n = 0; n < N; n++)
      {
        std::cout << n << ": ";
        for(double val : coeffs[n])
        {
          std::cout << val << ", ";
        }
        std::cout << std::endl;
      }
    };

  private:
    inline static double Eval(unsigned int Nc, const double x, const double * const c)
    { // Recursive Horner's scheme polynomial evaluation.
      if( Nc > 1 ) { return c[0] + x*Eval(Nc-1,x,c+1); }
      else if( Nc == 1 ) { return c[0]; }
      else { return 0; }
    };

    unsigned int N; // I could just use coeffs[0].size()...
    std::vector<std::vector<double> > coeffs; 

  }; // End of class Poly

  class Spline
  { // This class represents a piecewise-defined polynomial with a certain smoothness
    // or continuity enforced.  The number of knots (the junctions between piecewise
    // segments) and X position of the knots must be specified ahead of time, but the 
    // coefficients of the polynomials may be specified later.  The spline then allows
    // for convenient calculation of the spline Y values and derivatives.
    // Depending on the smoothness required, each non-leading segment has a number of 
    // dependent coefficients.  These are automatically calculated in the constructor,
    // in SetCoeffs, and in the operator() and Eval overloads that take a cnew parameter.
  public:
    // The user is responsible for ensuring that "knots" has N_knots values and that they
    // are sorted.  The user is also responsible for ensuring that "coeffs_" has enough 
    // values for the requested spline.  The number of coefficients can be obtained from the
    // static method GetNCoeffsTotal.
    Spline() : N_knots(0), N(0), smoothness(0), n_c_total(0), n_c_seg(0) {};
    Spline(unsigned int N_knots_, unsigned int N_, unsigned int smoothness_,
           const double * const knots, const double * const coeffs_ = nullptr) : 
      N_knots(N_knots_), N(N_), smoothness(smoothness_),
      n_c_total((N-(smoothness+1))*N_knots + N), n_c_seg(N-(smoothness+1)), 
      knots_x(knots,knots+N_knots),
      //coeffs(coeffs_,coeffs_+n_c_total),
      scratch(n_c_total)
    {
      assert(N > smoothness && "smoothness argument cannot exceed poly coefficient count.");
      // fill coeffs with zeros if coefficients not provided.
      if(!coeffs_) coeffs.assign(n_c_total,0); 
      else coeffs.assign(coeffs_,coeffs_+n_c_total); 
      Poly junk(N,coeffs.data()); // Temporary Poly object with coefficients.
      segments.assign(N_knots+1,junk);  // Fill segments array with copies of junk.
      SetCoeffs(coeffs.data()); // Now calculate the correct coeffs for all segments.
    };

    // Convenience overloads.
    double operator()(const double * const x, const double * const cnew) 
    { return Eval(x[0],cnew);};
    double operator()(const double * const x) const { return Eval(x[0]); };
    double Eval(const double * const x, const double * const cnew)  
    { return Eval(x[0],cnew);};
    
    // This Eval looks up the correct segment and returns the spline Y-value there.
    double Eval(double x) const { return segments[GetSegment(x)].Eval(x); };

    // This Eval takes a cnew argument to change the coefficients.  It only updates
    // the derived coefficients if necessary.
    double Eval(double x, const double * const cnew)
    {
      TryNewCoeffs(cnew);
      return Eval(x);
    };

    // The deriv methods do the same thing as the Eval, but for derivatives.
    double Deriv(unsigned int n, const double * const x, const double * const cnew) 
    { return Deriv(n,x[0],cnew);};
    double Deriv(unsigned int n, double x) const
    {
      return segments[GetSegment(x)].Deriv(n,x);
    };
    double Deriv(unsigned int n, double x, const double * const cnew)
    {
      TryNewCoeffs(cnew);
      return Deriv(n,x);
    }

    void SetCoeffs(const double * const newcoeffs)
    { // Calculate new dependent coefficients for non-leading segments.

      // First we just copy the first N coefficients for the leading segment.
      segments[0].SetCoeffs(newcoeffs);
      // This loop is over each of the other dependent segments.
      for(unsigned int i_seg = 1; i_seg < N_knots+1; i_seg++)
      { // x_knot is the x value at the knot on the left of this segment.
        double x_knot = knots_x[i_seg-1]; 
        
        // copy over the independent coefficients for this segment.  This depends
        // on the order and smoothness.  The copying is done into the high-order
        // terms of the scratch array.
        std::copy(newcoeffs+N+n_c_seg*(i_seg-1),
                  newcoeffs+N+n_c_seg*i_seg,
                  scratch.data()+smoothness+1);
        // Temporarily fill the low-order terms of the scratch array with zeros.
        std::fill(scratch.data(),scratch.data()+smoothness+1,0);

        Poly & poly = segments[i_seg]; // Get a reference to Poly for this segment.
        // Set the segment poly to use the partially-filled coeffs.
        poly.SetCoeffs(scratch.data());
        for(unsigned int s = smoothness+1; s --> 0;) // "s goes to zero"
        {
          // This mysterious calculation & loop enforces the smoothness by setting
          // the low-order terms of the segment (previously set to zero) to the values
          // that make it match the previous segment.
          scratch[s] = segments[i_seg-1].Deriv(s,x_knot) - poly.Deriv(s,x_knot);
          poly.SetCoeffs(scratch.data());
        }
      }
      // I don't actually use the coeffs vector for anything, but it's useful for
      // recalling the current coefficients without needing to extract them from
      // the Poly segments that use the real coefficients.
      coeffs.assign(newcoeffs,newcoeffs+n_c_seg);
    };

    unsigned int GetSegment(double x) const
    { // Return index of segment for given x value.
      return std::lower_bound(knots_x.begin(), knots_x.begin()+N_knots,x) - knots_x.begin();
    }

    // Return pointer to current coeffs.
    double const * GetCoeffs() const { return coeffs.data(); };

    // Print values of all coefficients for all segments including derivatives.
    void PrintCoeffs() const
    {
      for(unsigned int n_seg = 0; n_seg < segments.size(); n_seg++)
      {
        std::cout << "Segment " << n_seg << ":" << std::endl;
        segments[n_seg].PrintCoeffs();
      }
    };

    // Boring getters for various private values.
    // I would have made the values public & const, but that makes the class
    // difficult to copy/assign and stuff, so this is easier.
    unsigned int GetN_knots() const {return N_knots;};
    unsigned int GetN() const {return N;};
    unsigned int GetSmoothness() const {return smoothness;};
    unsigned int GetNCoeffsTotal() const {return n_c_total;};
    unsigned int GetNCoeffsPerSegment() const {return n_c_seg;};
    std::vector<double> GetKnots() const {return knots_x;};

    // Static overload to give the number of coefficients required to construct
    // a Spline with the given number of knots, order, and smoothness.
    static unsigned int GetNCoeffsTotal(unsigned int N_knots, 
                                        unsigned int N, 
                                        unsigned int smoothness)
    { return (N-(smoothness+1))*N_knots + N ;};


    TF1 makeTF1()
    {
      double low,high;
      if(N_knots > 1)
      {
        low = knots_x[0] - (knots_x[1] - knots_x[0]);
        high= knots_x[N_knots-1] + (knots_x[N_knots-1] - knots_x[N_knots-2]);
      }
      else
      { // Deal with N_knots == 0 or 1 in the same case.
        low = N_knots == 1 ? knots_x[0] - 1 : -1;
        high = low + 2;
      }

      TF1 fspline("fspline",this,low,high,n_c_total);
      fspline.SetParameters(coeffs.data());
      return fspline;
    };

private:
    unsigned int N_knots;
    unsigned int N;
    unsigned int smoothness;

    unsigned int n_c_total;
    unsigned int n_c_seg;
    std::vector<double> knots_x;
    std::vector<Poly> segments;
    std::vector<double> coeffs;
    std::vector<double> scratch;



    void TryNewCoeffs(const double * const cnew)
    {
      const double * const cold = coeffs.data();
      for(unsigned int i = 0; i < n_c_total; i++)
      {
        if(cold[i] != cnew[i])
        {
          SetCoeffs(cnew);
          return;
        }
      }
    }

  }; // End of class Spline

  ///////////////////////////////////////////////////////////////////////////////////
  // Utility functions for MySpline2 ////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////////////
  bool CheckHeader(const TObjArray * const toa,unsigned int index, TString expected)
  {
    TString present = TString(toa->At(index)->GetName()).Strip(TString::kBoth,' ');
    bool result = (present == expected);
    if(not result)
      {
        std::cout << present << "," << expected << std::endl;
      }
    return present == expected;
  }
  
  //PieceWise::PieceWise(TString fname)
  PieceWise GetFromFile(TString fname)
  {
    std::ifstream myfile(fname);
    if(!myfile) 
      {
        std::cout << "Error, file " << fname << " not found." << std::endl;
        return PieceWise();
      }
    std::string stdline;
    std::getline(myfile,stdline); // Skip past the first header line.
    std::getline(myfile,stdline); // This header line defines the columns.
    TString line(stdline.c_str());
    line = line.Strip(TString::kLeading,'#');
    line = line.Strip(TString::kBoth,' ');
    auto header_toa = line.Tokenize(',');
    
    std::getline(myfile,stdline);
    line = stdline.c_str();
    auto values_toa = line.Tokenize(",");
    const auto nentries = header_toa->GetEntries();
    assert(nentries == values_toa->GetEntries());
    
    
    assert(CheckHeader(header_toa,0,"order/i"));
    const unsigned int order = TString(values_toa->At(0)->GetName()).Atoi();
    assert(CheckHeader(header_toa,1,"smoothness/i"));
    const unsigned int smoothness = TString(values_toa->At(1)->GetName()).Atoi();
    assert(CheckHeader(header_toa,2,"n_knots/i"));
    const unsigned int n_knots = TString(values_toa->At(2)->GetName()).Atoi();
    double * knots = new double[n_knots];
    
    for(unsigned int i = 0; i < n_knots; i++)
      {
        assert(CheckHeader(header_toa,3+i,TString::Format("knot%d/D",i)));
        knots[i] = (TString(values_toa->At(3+i)->GetName()).Atof());
      }
    PieceWise pw(n_knots,order,smoothness,knots);
    unsigned int n_coeffs = pw.n_c_total;
    
    double * coeffs = new double[n_coeffs];
    for(unsigned int i = 0; i < n_coeffs; i++)
      {
        assert(CheckHeader(header_toa,3+n_knots+i,TString::Format("c%d/D",i)));
        coeffs[i] = (TString(values_toa->At(3+n_knots+i)->GetName()).Atof());
      }
    
    myfile.close();

    pw.SetCoeffs(coeffs);
    delete[] knots;
    delete[] coeffs;
    return pw;
  }

  //enum Side {LEFT, RIGHT, BOTH};
  //const TString SideStrings[] = {"LEFT","RIGHT","BOTH"};
  
  TString GetFile(int runnum, int order, Geom::Side side = Geom::BOTH,
                        const TString init = "Spline_Correction_")
  {
    // This part of the function locates the correct spline correction calib file.
    const TSystemDirectory tsd("tsd","calib/Spline_Correction");
    const TList * cal_list = tsd.GetListOfFiles();
    const int nfiles = cal_list->GetEntries();
    //const TString init("Spline_Correction");
    std::vector<int> filetimes; // This vector will store the desired file timestamps.
    for(int i = 0;i< nfiles;i++)
      {
        const TString fname(cal_list->At(i)->GetName());

        // Reject any files that don't start with the right form.
        if(not fname.BeginsWith(init)) continue;

        // This splits the string at each of the tokens listed.
        const TObjArray * toa = fname.Tokenize("_.");
        const Int_t filerun = TString(toa->At(2)->GetName()).Atoi();
        const Int_t fileorder = TString(toa->At(3)->GetName()).Atoi();
        const TString fileside = TString(toa->At(4)->GetName());
        const Int_t filetime = TString(toa->At(5)->GetName()).Atoi();

        // Reject any files that don't have the desired run number or correction order.
        TString filesideUP = fileside;
        filesideUP.ToUpper();
        if((filerun != runnum) or 
           (fileorder != order) or 
           (filesideUP != Geom::SideStrings[side])) continue;
        

        // Get the timestamp from the filename (1 second precision).
        filetimes.push_back(filetime);
      }
    // Here I should delete the temporary TLists and things from the file searching above.
    
    // Now filetimes contains the list of timestamps for the desired files, 
    // but there might be 0.
    if(filetimes.size() <= 0)
      {
        std::cout << "Could not find a " << init << " file for run " << runnum 
                  << " and order " << order << std::endl;
        return "";
      }
    
    // Now of the desired files, I take the latest one (max timestamp)
    const Int_t latest = *std::max_element(filetimes.begin(),filetimes.end());

    // Reconstruct the filename, it's easier than actually holding a file handle.
    const TString cal_fname = TString::Format("%s%05d_%d_%s_%d.csv",
                                              init.Data(),runnum,order,
                                              Geom::SideStrings[side].Data(),latest);

    assert(bool(cal_list->FindObject(cal_fname)));
    std::cout << "Found file " << cal_fname << " for " << init << std::endl;
    return TString::Format("calib/Spline_Correction/%s",cal_fname.Data());
  }

  Spline GetSpline(TString fname)
  {
    std::ifstream myfile(fname);
    if(!myfile) 
      {
        std::cout << "Error, file " << fname << " not found." << std::endl;
        return Spline();
      }
    std::string stdline;
    std::getline(myfile,stdline); // Skip past the first header line.
    std::getline(myfile,stdline); // This header line defines the columns.
    TString line(stdline.c_str());
    line = line.Strip(TString::kLeading,'#');
    line = line.Strip(TString::kBoth,' ');
    auto header_toa = line.Tokenize(',');
    
    std::getline(myfile,stdline);
    line = stdline.c_str();
    auto values_toa = line.Tokenize(",");
    const unsigned int nentries = header_toa->GetEntries();
    assert(nentries == static_cast<unsigned int>(values_toa->GetEntries()));
    
    
    assert(CheckHeader(header_toa,0,"order/i"));
    const unsigned int order = TString(values_toa->At(0)->GetName()).Atoi();
    assert(CheckHeader(header_toa,1,"smoothness/i"));
    const unsigned int smoothness = TString(values_toa->At(1)->GetName()).Atoi();
    assert(CheckHeader(header_toa,2,"n_knots/i"));
    const unsigned int n_knots = TString(values_toa->At(2)->GetName()).Atoi();
    //double * knots = new double[n_knots];
    std::vector<double> knots;//(n_knots);

    for(unsigned int i = 0; i < n_knots; i++)
      {
        assert(CheckHeader(header_toa,3+i,TString::Format("knot%d/D",i)));
        //knots[i] = (TString(values_toa->At(3+i)->GetName()).Atof());
        knots.push_back((TString(values_toa->At(3+i)->GetName()).Atof()));
      }

    std::vector<double> coeffs;
    unsigned int i = 0;
    while(3+n_knots+i < nentries)
      {
        assert(CheckHeader(header_toa,3+n_knots+i,TString::Format("c%d/D",i)));
        coeffs.push_back(TString(values_toa->At(3+n_knots+i)->GetName()).Atof());
        i++;
      }
    
    myfile.close();

    Spline spline(n_knots,order,smoothness,knots.data(),coeffs.data());
    //delete[] knots;
    
    return spline;
  }

  Spline GetSpline(int runnum, int order, Geom::Side side = Geom::BOTH,
                   const TString init = "Spline_Correction_")
  {
    return GetSpline(GetFile(runnum,order,side,init));
  }

  TF1 makeTF1(PieceWise & pw)
  {
    double min,max;
    if(pw.N_knots == 0)
    {
      min = -1;
      max = 1;
    }
    else if(pw.N_knots == 1)
    {
      min = pw.knots_x[0]-1;
      max = min + 2;
    }
    else
    {
      double deltalow = pw.knots_x[1] - pw.knots_x[0];
      double deltahigh = pw.knots_x[pw.N_knots-1] - pw.knots_x[pw.N_knots-2];
      min = pw.knots_x[0] - deltalow;
      max = pw.knots_x[pw.N_knots-1] + deltahigh;
    }
    TF1 f("f_spline",pw,min,max,pw.n_c_total);
    f.SetParameters(pw.coeffs.data());
    return f;
  }

  class CorrectionVerifier
  {
  public:
    CorrectionVerifier(){slope = 0.0;};
    CorrectionVerifier(Spline & corr)
    {
      slope = corr.Eval(minimum_time)/minimum_time;
    };

    double Verify(double arr_time, double proposed) const
    {
      if(TMath::Abs(arr_time) > minimum_time)
      {
        return proposed;
      }
      else // This also catches NaN proposals.
      {
        return slope*arr_time;
      }
    };

    // The initialization of this static const data member cannot be in the class
    // definition for some reason, so look for it below this class.
    static const double minimum_time;
  private:
    double slope;
  };

  // For some reason the initialization of this static const data member
  // cannot be in the class definition, so I wrote it here.
  const double CorrectionVerifier::minimum_time = 0.01; // microseconds
  
} // End of namespace MySpline2

#if defined(__MAKECINT__)
#pragma link C++ defined_in MySpline2;
#pragma link C++ nestedclass;
//#pragma link C++ class MySpline2::MySpline2;

#pragma link C++ class MySpline2::PieceWise;
#pragma link C++ class MySpline2::Poly;
#pragma link C++ class MySpline2::Spline;

#pragma link C++ class MySpline2::CorrectionVerifier;
#endif

#endif // include guard