// https://root-forum.cern.ch/t/using-tspline-for-fit/18255/2

#include <stdio.h>
#include <iostream>
#include <fstream>
#include <string>
#include <TLegend.h>
#include <TMath.h>


using namespace std;

#include <TCanvas.h>
#include <TH1.h>
#include <TF1.h>
#include <TFile.h>
#include <TTree.h>
#include <TChain.h>
#include <TNtuple.h>
#include <TLegend.h>

#include <Math/Interpolator.h>


Double_t spline_4nodes(Double_t *x, Double_t *par)
{
   /*Fit parameters:
   par[0-3]=X of nodes (to be fixed in the fit!)
   par[4-7]=Y of nodes
   par[8-9]=first derivative at begin and end (to be fixed in the fit!)
   */
   
   
    Double_t xx = x[0];

   Double_t xn[5] = { par[0], par[1], par[2], par[3], par[4]  };
   
   Double_t yn[5] = { par[5], par[6], par[7], par[8], par[9] };


//   TSpline3 sp3("sp3", xn, yn, 5, "b1e1", b1, e1);

 //  return sp3.Eval(xx);

    ROOT::Math::Interpolator inter (5, ROOT::Math::Interpolation:: kPOLYNOMIAL );
/*
  Interpolator (
    unsigned int ndata,                   : how many data points
    ROOT::Math::Interpolation::Type type  : kLINEAR, kPOLYNOMIAL, kCSPLINE, kAKIMA, kCSPLINE_PERIODIC, kAKIMA_PERIODIC
  )

  Setting up the interpolation tool. Choose among the following methods:
   
  o kLINEAR           : straight line between points
  o kPOLYNOMIAL       : to be used for small number of points since introduces large oscillations
  o kCSPLINE          : cubic spline with natural  boundary conditions
  o kCSPLINE_PERIODIC : cubic spline with periodic boundary conditions
  o kAKIMA            : Akima spline with natural  boundary conditions ( requires min. 5 points)
  o kAKIMA_PERIODIC   : Akima spline with periodic boundaries          ( requires min. 5 points)

  See: root.cern.ch/root/html/ROOT__Math__Interpolator.html#ROOT__Math__Interpolator:Interpolator
*/

 
  inter.SetData (5,  xn,  yn);
  
  return inter.Eval (xx);


}



void macro_interpolation ()
{
	
  // gROOT->SetStyle("Plain");   //set plain TStyle
   gStyle->SetOptStat(0); //draw statistics on plots, 0 for no output
 // gStyle->SetOptStat(111111); //draw statistics on plots, 0 for no output
   gStyle->SetOptFit(1);    //draw fit results on plot, 0 for no output
   gStyle->SetPalette(57);     //set color map
   gStyle->SetOptTitle(0);    //suppress title box



  


  float xmin = -1.0, xmax = 1.0; 
  
  /*
  int Nprob = 100;
  double* Xprob  = new double [Nprob];
  double* Yinter = new double [Nprob];
  
  for ( int i = 0; i < Nprob; ++i ) {
    Xprob [i] = (double) i*(xmax-xmin)/(Nprob-1) + xmin;
    Yinter[i] = inter.Eval ( Xprob[i] );
  //  cout << Yinter[i]  << endl ; 
  }
 */
  TCanvas* can1 = new TCanvas ("c1", "Interpolation", 600, 400); 
 
  TGraphErrors *gf = new TGraphErrors("C:/root_v5.34.36/potential/beta/C6.txt", "%lf,%lf,%lf,%lf\n");
  gf->GetXaxis()->SetRangeUser(-1, 1);
  gf->SetLineColor(14);
  gf->SetLineWidth(3);
  gf->SetMarkerStyle(22);
  gf->SetTitle ("Triangles: input.  Curve: probed interpolator response");
  gf->Draw ("AP");
  
  
 int fit_pts_x_y_dy_dx = 10 ;
 
  double x0 = gf->GetX()[0];  // phonon sigma for k = 1
 double x1 = gf->GetX()[1];  // Intensity for k =1  
 double x2 = gf->GetX()[2];  // zero-phonon line for k =1  
 double x3 = gf->GetX()[3];  //mode frequency for k = 1
 double x4 = gf->GetX()[4];  // phonon sigma for k = 1
 double x5 = gf->GetX()[5];  // Intensity for k =1  
 double x6 = gf->GetX()[6];  // zero-phonon line for k =1  
 double x7 = gf->GetX()[7];  //mode frequency for k = 1
 double x8 = gf->GetX()[8];  // phonon sigma for k = 1
 double x9 = gf->GetX()[9];  // Intensity for k =1  
 double x10 = gf->GetX()[10];  // zero-phonon line for k =1  
 double x11 = gf->GetX()[11];  //mode frequency for k = 1
 double x12 = gf->GetX()[12];  // phonon sigma for k = 1
 double x13 = gf->GetX()[13];  // Intensity for k =1  
 double x14 = gf->GetX()[14];  // zero-phonon line for k =1  
 double x15 = gf->GetX()[15];  //mode frequency for k = 1
 double x16 = gf->GetX()[16];  // phonon sigma for k = 1
 double x17 = gf->GetX()[17];  // Intensity for k =1  
 
 double y0 = gf->GetY()[0];  // phonon sigma for k = 1
 double y1 = gf->GetY()[1];  // Intensity for k =1  
 double y2 = gf->GetY()[2];  // zero-phonon line for k =1  
 double y3 = gf->GetY()[3];  //mode frequency for k = 1
 double y4 = gf->GetY()[4];  // phonon sigma for k = 1
 double y5 = gf->GetY()[5];  // Intensity for k =1  
 double y6 = gf->GetY()[6];  // zero-phonon line for k =1  
 double y7 = gf->GetY()[7];  //mode frequency for k = 1
 double y8 = gf->GetY()[8];  // phonon sigma for k = 1
 double y9 = gf->GetY()[9];  // Intensity for k =1  
 double y10 = gf->GetY()[10];  // zero-phonon line for k =1  
 double y11 = gf->GetY()[11];  //mode frequency for k = 1
 double y12 = gf->GetY()[12];  // phonon sigma for k = 1
 double y13 = gf->GetY()[13];  // Intensity for k =1  
 double y14 = gf->GetY()[14];  // zero-phonon line for k =1  
 double y15 = gf->GetY()[15];  //mode frequency for k = 1
 double y16 = gf->GetY()[16];  // phonon sigma for k = 1
 double y17 = gf->GetY()[17];  // Intensity for k =1  
 
 
 double x_0 = x0;  // phonon sigma for k = 1
 double x_1 = x4;  // Intensity for k =1  
 double x_2 = x8;  // Intensity for k =1  
 double x_3 = x13;  // zero-phonon line for k =1  
 double x_4 = x17;  //mode frequency for k = 1
 
 double y_0 = y0;  //mode frequency for k = 1
 double y_1 = y4;  //mode frequency for k = 1
 double y_2 = y8;  //mode frequency for k = 1
 double y_3 = y13;  //mode frequency for k = 1
 double y_4 = y17;  //mode frequency for k = 1
 
 double b_1 = 0E-5;  //mode frequency for k = 1
 double e_1 = 0E-5;  //mode frequency for k = 1
 
  TF1 *f_spline4 = new TF1("f_spline4", spline_4nodes, xmin, xmax, fit_pts_x_y_dy_dx); // npars = 2*nodes+2
  
     f_spline4->SetParameter(0, x_0);
     f_spline4->SetParameter(1, x_1);
	 f_spline4->SetParameter(2, x_2);
	 f_spline4->SetParameter(3, x_3);
	 f_spline4->SetParameter(4, x_4);
	 
	 
	 f_spline4->SetParameter(5, y_0);
	 f_spline4->SetParameter(6, y_1);
	 f_spline4->SetParameter(7, y_2);
	 f_spline4->SetParameter(8, y_3);
	 f_spline4->SetParameter(9, y_4);

     f_spline4->SetParName(0, "x_0");
     f_spline4->SetParName(1, "x_1");
	 f_spline4->SetParName(2,"x_2");
	 f_spline4->SetParName(3,"x_3");
	 f_spline4->SetParName(4,"x_4");
	 	
	 f_spline4->SetParName(5, "y_0");
     f_spline4->SetParName(6, "y_1");
	 f_spline4->SetParName(7,"y_2");
	 f_spline4->SetParName(8,"y_3");
	 f_spline4->SetParName(9,"y_4");
	 
	 

  gf->Fit(f_spline4); 
  gf->SetLineColor (28);
  gf->SetLineWidth (1);
  gf->SetTitle ("Integral");
  gf->Draw ("SAME L"); 

 /*  const int eval_pts = 5;
   Double_t x_eval[eval_pts] = {  f_spline4->GetParameter(0),  f_spline4->GetParameter(1),  f_spline4->GetParameter(2),  f_spline4->GetParameter(3),  f_spline4->GetParameter(4) };
   Double_t y_eval[eval_pts] = {  f_spline4->GetParameter(5),  f_spline4->GetParameter(6),  f_spline4->GetParameter(7),  f_spline4->GetParameter(8),  f_spline4->GetParameter(9) };

   Double_t b1_eval =  f_spline4->GetParameter(10);
   Double_t e1_eval =  f_spline4->GetParameter(11);
   
   TSpline3 sp_eval("sp_eval", x_eval, y_eval, eval_pts, "b1e1_eval", b1_eval, e1_eval);

   cout << sp_eval.Eval(-0.75) <<endl;
*/
}