Hello,
I am trying to use the TSpline(3/5) class to parametrize some data points and then use the resulting TSpline object to retrieve the higher oder derivatives.
Reading the (somewhat poor) documentation of the TSpline(3/5) classes, I thought that the spline is build from data at knots requiring that the spline derivatives are continuous up to 2nd order for a TSpline3 and up to 4th order for a TSpline5. I am probably missing something or made a bug in the following simple macro which obviously give non continuous derivatives already at 1st order.
Thanks for any suggestion.
//_______________________________________________________________________________________________________
void GetTaylorTermsAtX(TSpline3 *spl, double x, double *SExpansion){
// SExpansion array of size 4 (value of function and 3 derivatives at point x
int i = spl->FindX(x);
double Xi;
double A0, A1, A2, A3;// An = d^n spl(x) / dx^n / n!
spl->GetCoeff(i, Xi, A0, A1, A2, A3);
double p = x - Xi;
SExpansion[0] = A0 + A1 * p + A2 * p * p + A3 * p * p * p;
SExpansion[1] = A1 + 2 * A2 * p + 3 * A3 * p * p;
SExpansion[2] = (2 * A2 + 6 * A3 * p) / 2.;
SExpansion[3] = A3;
return;
}// end of function GTaylorTermsAtX
//_______________________________________________________________________________________________________
//_______________________________________________________________________________________________________
void GetTaylorTermsAtX(TSpline5 *spl, double x, double *SExpansion){
// SExpansion array of size 6 (value of function and 5 derivatives at point x
int i = spl->FindX(x);
double Xi;
double A0, A1, A2, A3, A4, A5;// An = d^n spl(x) / dx^n / n!
spl->GetCoeff(i, Xi, A0, A1, A2, A3, A4, A5);
double p = x - Xi;
SExpansion[0] = A0 + A1 * p + A2 * p * p + A3 * p * p * p + A4 * p * p * p * p + A5 * p * p * p * p * p;
SExpansion[1] = A1 + 2 * A2 * p + 3 * A3 * p * p + 4 * A4 * p * p * p + 5 * A5 * p * p * p * p;
SExpansion[2] = (2 * A2 + 6 * A3 * p + 12 * A4 * p * p + 20 * A5 * p * p * p ) / 2.;
SExpansion[3] = (6 * A3 + 24 * A4 * p + 60 * A5 * p * p) / 6.;
SExpansion[4] = (24 * A4 + 120 * A5 * p) / 24.;
SExpansion[5] = A5;// (120 / 120)
return;
}// end of function GTaylorTermsAtX
//_______________________________________________________________________________________________________
//_______________________________________________________________________________________________________
void test_spline(){
const Int_t npts = 100;
Double_t theta[npts], GE2[npts];
theta[ 0] = 2.44346e+00; GE2[ 0] = 9.05542e-01;
theta[ 1] = 2.44364e+00; GE2[ 1] = 9.05537e-01;
theta[ 2] = 2.44381e+00; GE2[ 2] = 9.05532e-01;
theta[ 3] = 2.44398e+00; GE2[ 3] = 9.05527e-01;
theta[ 4] = 2.44416e+00; GE2[ 4] = 9.05522e-01;
theta[ 5] = 2.44433e+00; GE2[ 5] = 9.05517e-01;
theta[ 6] = 2.44451e+00; GE2[ 6] = 9.05512e-01;
theta[ 7] = 2.44468e+00; GE2[ 7] = 9.05507e-01;
theta[ 8] = 2.44486e+00; GE2[ 8] = 9.05503e-01;
theta[ 9] = 2.44503e+00; GE2[ 9] = 9.05498e-01;
theta[ 10] = 2.44521e+00; GE2[ 10] = 9.05493e-01;
theta[ 11] = 2.44538e+00; GE2[ 11] = 9.05488e-01;
theta[ 12] = 2.44556e+00; GE2[ 12] = 9.05483e-01;
theta[ 13] = 2.44573e+00; GE2[ 13] = 9.05478e-01;
theta[ 14] = 2.44590e+00; GE2[ 14] = 9.05474e-01;
theta[ 15] = 2.44608e+00; GE2[ 15] = 9.05469e-01;
theta[ 16] = 2.44625e+00; GE2[ 16] = 9.05464e-01;
theta[ 17] = 2.44643e+00; GE2[ 17] = 9.05459e-01;
theta[ 18] = 2.44660e+00; GE2[ 18] = 9.05454e-01;
theta[ 19] = 2.44678e+00; GE2[ 19] = 9.05449e-01;
theta[ 20] = 2.44695e+00; GE2[ 20] = 9.05444e-01;
theta[ 21] = 2.44713e+00; GE2[ 21] = 9.05440e-01;
theta[ 22] = 2.44730e+00; GE2[ 22] = 9.05435e-01;
theta[ 23] = 2.44748e+00; GE2[ 23] = 9.05430e-01;
theta[ 24] = 2.44765e+00; GE2[ 24] = 9.05425e-01;
theta[ 25] = 2.44782e+00; GE2[ 25] = 9.05420e-01;
theta[ 26] = 2.44800e+00; GE2[ 26] = 9.05415e-01;
theta[ 27] = 2.44817e+00; GE2[ 27] = 9.05411e-01;
theta[ 28] = 2.44835e+00; GE2[ 28] = 9.05406e-01;
theta[ 29] = 2.44852e+00; GE2[ 29] = 9.05401e-01;
theta[ 30] = 2.44870e+00; GE2[ 30] = 9.05396e-01;
theta[ 31] = 2.44887e+00; GE2[ 31] = 9.05391e-01;
theta[ 32] = 2.44905e+00; GE2[ 32] = 9.05386e-01;
theta[ 33] = 2.44922e+00; GE2[ 33] = 9.05382e-01;
theta[ 34] = 2.44940e+00; GE2[ 34] = 9.05377e-01;
theta[ 35] = 2.44957e+00; GE2[ 35] = 9.05372e-01;
theta[ 36] = 2.44974e+00; GE2[ 36] = 9.05367e-01;
theta[ 37] = 2.44992e+00; GE2[ 37] = 9.05362e-01;
theta[ 38] = 2.45009e+00; GE2[ 38] = 9.05357e-01;
theta[ 39] = 2.45027e+00; GE2[ 39] = 9.05353e-01;
theta[ 40] = 2.45044e+00; GE2[ 40] = 9.05348e-01;
theta[ 41] = 2.45062e+00; GE2[ 41] = 9.05343e-01;
theta[ 42] = 2.45079e+00; GE2[ 42] = 9.05338e-01;
theta[ 43] = 2.45097e+00; GE2[ 43] = 9.05333e-01;
theta[ 44] = 2.45114e+00; GE2[ 44] = 9.05328e-01;
theta[ 45] = 2.45131e+00; GE2[ 45] = 9.05324e-01;
theta[ 46] = 2.45149e+00; GE2[ 46] = 9.05319e-01;
theta[ 47] = 2.45166e+00; GE2[ 47] = 9.05314e-01;
theta[ 48] = 2.45184e+00; GE2[ 48] = 9.05309e-01;
theta[ 49] = 2.45201e+00; GE2[ 49] = 9.05304e-01;
theta[ 50] = 2.45219e+00; GE2[ 50] = 9.05300e-01;
theta[ 51] = 2.45236e+00; GE2[ 51] = 9.05295e-01;
theta[ 52] = 2.45254e+00; GE2[ 52] = 9.05290e-01;
theta[ 53] = 2.45271e+00; GE2[ 53] = 9.05285e-01;
theta[ 54] = 2.45289e+00; GE2[ 54] = 9.05280e-01;
theta[ 55] = 2.45306e+00; GE2[ 55] = 9.05276e-01;
theta[ 56] = 2.45323e+00; GE2[ 56] = 9.05271e-01;
theta[ 57] = 2.45341e+00; GE2[ 57] = 9.05266e-01;
theta[ 58] = 2.45358e+00; GE2[ 58] = 9.05261e-01;
theta[ 59] = 2.45376e+00; GE2[ 59] = 9.05256e-01;
theta[ 60] = 2.45393e+00; GE2[ 60] = 9.05252e-01;
theta[ 61] = 2.45411e+00; GE2[ 61] = 9.05247e-01;
theta[ 62] = 2.45428e+00; GE2[ 62] = 9.05242e-01;
theta[ 63] = 2.45446e+00; GE2[ 63] = 9.05237e-01;
theta[ 64] = 2.45463e+00; GE2[ 64] = 9.05232e-01;
theta[ 65] = 2.45481e+00; GE2[ 65] = 9.05228e-01;
theta[ 66] = 2.45498e+00; GE2[ 66] = 9.05223e-01;
theta[ 67] = 2.45515e+00; GE2[ 67] = 9.05218e-01;
theta[ 68] = 2.45533e+00; GE2[ 68] = 9.05213e-01;
theta[ 69] = 2.45550e+00; GE2[ 69] = 9.05208e-01;
theta[ 70] = 2.45568e+00; GE2[ 70] = 9.05204e-01;
theta[ 71] = 2.45585e+00; GE2[ 71] = 9.05199e-01;
theta[ 72] = 2.45603e+00; GE2[ 72] = 9.05194e-01;
theta[ 73] = 2.45620e+00; GE2[ 73] = 9.05189e-01;
theta[ 74] = 2.45638e+00; GE2[ 74] = 9.05184e-01;
theta[ 75] = 2.45655e+00; GE2[ 75] = 9.05180e-01;
theta[ 76] = 2.45673e+00; GE2[ 76] = 9.05175e-01;
theta[ 77] = 2.45690e+00; GE2[ 77] = 9.05170e-01;
theta[ 78] = 2.45707e+00; GE2[ 78] = 9.05165e-01;
theta[ 79] = 2.45725e+00; GE2[ 79] = 9.05161e-01;
theta[ 80] = 2.45742e+00; GE2[ 80] = 9.05156e-01;
theta[ 81] = 2.45760e+00; GE2[ 81] = 9.05151e-01;
theta[ 82] = 2.45777e+00; GE2[ 82] = 9.05146e-01;
theta[ 83] = 2.45795e+00; GE2[ 83] = 9.05141e-01;
theta[ 84] = 2.45812e+00; GE2[ 84] = 9.05137e-01;
theta[ 85] = 2.45830e+00; GE2[ 85] = 9.05132e-01;
theta[ 86] = 2.45847e+00; GE2[ 86] = 9.05127e-01;
theta[ 87] = 2.45865e+00; GE2[ 87] = 9.05122e-01;
theta[ 88] = 2.45882e+00; GE2[ 88] = 9.05118e-01;
theta[ 89] = 2.45899e+00; GE2[ 89] = 9.05113e-01;
theta[ 90] = 2.45917e+00; GE2[ 90] = 9.05108e-01;
theta[ 91] = 2.45934e+00; GE2[ 91] = 9.05103e-01;
theta[ 92] = 2.45952e+00; GE2[ 92] = 9.05099e-01;
theta[ 93] = 2.45969e+00; GE2[ 93] = 9.05094e-01;
theta[ 94] = 2.45987e+00; GE2[ 94] = 9.05089e-01;
theta[ 95] = 2.46004e+00; GE2[ 95] = 9.05084e-01;
theta[ 96] = 2.46022e+00; GE2[ 96] = 9.05080e-01;
theta[ 97] = 2.46039e+00; GE2[ 97] = 9.05075e-01;
theta[ 98] = 2.46057e+00; GE2[ 98] = 9.05070e-01;
theta[ 99] = 2.46074e+00; GE2[ 99] = 9.05065e-01;
const Int_t ngr = 4;
TSpline3 *GE2_spline = new TSpline3("GE2_spline", theta, GE2, npts);
GE2_spline->SetLineColor(kBlue);
GE2_spline->SetLineWidth(2);
// ---------------------------------------------------------------- //
// -------*-*-*-*-*------ Derivatives of GE2 -------*-*-*-*-*------ //
// ---------------------------------------------------------------- //
TCanvas *cnv_GE2 = new TCanvas("cnv_GE2", "cnv_GE2", 2500, 1500);
cnv_GE2->Divide(ngr / 2, 2);
Double_t leftmargin = 0.15, rightmargin = 0.15;
Double_t SExpansion[ngr];
TString gr_title;
TGraph *gr_GE2_S[ngr];
for (Int_t igr = 0; igr < ngr; ++igr){
gr_GE2_S[igr] = new TGraph(npts);
}
for (Int_t ipt = 0; ipt < npts; ++ipt){
GetTaylorTermsAtX(GE2_spline, theta[ipt], SExpansion);
for (Int_t igr = 0; igr < ngr; ++igr){
gr_GE2_S[igr]->SetPoint(ipt, theta[ipt], SExpansion[igr]);
}
}
for (Int_t igr = 0; igr < ngr; ++igr){
cnv_GE2->cd(igr+1);
gPad->SetLeftMargin(leftmargin);
gPad->SetRightMargin(rightmargin);
gr_GE2_S[igr]->Draw("apl");
gr_title.Form(";#theta [rd];G_{E}^{2(%d)}", igr);
gr_GE2_S[igr]->SetTitle(gr_title.Data());
if (igr == 0) GE2_spline->Draw("Lsame");
}
}// end of test_spline function
//_______________________________________________________________________________________________________
_ROOT Version: 5.34 to 6.14
Platform: Ubuntu 18.04 LTS
Compiler: gcc 7.3