#include "TSpline.h"
#include "TVirtualPad.h"
#include "TF1.h"
#include "TSystem.h"
#include "Riostream.h"
#include "TClass.h"
ClassImp(TSplinePoly)
ClassImp(TSplinePoly3)
ClassImp(TSplinePoly5)
ClassImp(TSpline3)
ClassImp(TSpline5)
ClassImp(TSpline)
TSpline::TSpline(const TSpline &sp) :
TNamed(sp),
TAttLine(sp),
TAttFill(sp),
TAttMarker(sp),
fDelta(sp.fDelta),
fXmin(sp.fXmin),
fXmax(sp.fXmax),
fNp(sp.fNp),
fKstep(sp.fKstep),
fHistogram(sp.fHistogram),
fGraph(sp.fGraph),
fNpx(sp.fNpx)
{
}
TSpline& TSpline::operator=(const TSpline &sp)
{
if(this!=&sp) {
TNamed::operator=(sp);
TAttLine::operator=(sp);
TAttFill::operator=(sp);
TAttMarker::operator=(sp);
fDelta=sp.fDelta;
fXmin=sp.fXmin;
fXmax=sp.fXmax;
fNp=sp.fNp;
fKstep=sp.fKstep;
fHistogram=sp.fHistogram;
fGraph=sp.fGraph;
fNpx=sp.fNpx;
}
return *this;
}
void TSpline::Draw(Option_t *option)
{
TString opt = option;
opt.ToLower();
if (gPad && !opt.Contains("same")) gPad->Clear();
AppendPad(option);
}
void TSpline::Paint(Option_t *option)
{
Int_t i;
Double_t xv;
TString opt = option;
opt.ToLower();
Double_t xmin, xmax, pmin, pmax;
pmin = gPad->PadtoX(gPad->GetUxmin());
pmax = gPad->PadtoX(gPad->GetUxmax());
xmin = fXmin;
xmax = fXmax;
if (opt.Contains("same")) {
if (xmax < pmin) return;
if (xmin > pmax) return;
if (xmin < pmin) xmin = pmin;
if (xmax > pmax) xmax = pmax;
} else {
gPad->Clear();
}
if (fHistogram)
if ((!gPad->GetLogx() && fHistogram->TestBit(TH1::kLogX)) ||
(gPad->GetLogx() && !fHistogram->TestBit(TH1::kLogX)))
{ delete fHistogram; fHistogram = 0;}
if (fHistogram) {
fHistogram->GetXaxis()->SetLimits(xmin,xmax);
} else {
if (xmin > 0 && gPad->GetLogx()) {
Double_t *xbins = new Double_t[fNpx+1];
Double_t xlogmin = TMath::Log10(xmin);
Double_t xlogmax = TMath::Log10(xmax);
Double_t dlogx = (xlogmax-xlogmin)/((Double_t)fNpx);
for (i=0;i<=fNpx;i++) {
xbins[i] = gPad->PadtoX(xlogmin+ i*dlogx);
}
fHistogram = new TH1F("Spline",GetTitle(),fNpx,xbins);
fHistogram->SetBit(TH1::kLogX);
delete [] xbins;
} else {
fHistogram = new TH1F("Spline",GetTitle(),fNpx,xmin,xmax);
}
if (!fHistogram) return;
fHistogram->SetDirectory(0);
}
for (i=1;i<=fNpx;i++) {
xv = fHistogram->GetBinCenter(i);
fHistogram->SetBinContent(i,this->Eval(xv));
}
fHistogram->SetBit(TH1::kNoStats);
fHistogram->SetLineColor(GetLineColor());
fHistogram->SetLineStyle(GetLineStyle());
fHistogram->SetLineWidth(GetLineWidth());
fHistogram->SetFillColor(GetFillColor());
fHistogram->SetFillStyle(GetFillStyle());
fHistogram->SetMarkerColor(GetMarkerColor());
fHistogram->SetMarkerStyle(GetMarkerStyle());
fHistogram->SetMarkerSize(GetMarkerSize());
char *o = (char *) opt.Data();
Int_t j=0;
i=0;
Bool_t graph=kFALSE;
do
if(o[i]=='p') graph=kTRUE ; else o[j++]=o[i];
while(o[i++]);
if (opt.Length() == 0 ) fHistogram->Paint("lf");
else if (opt == "same") fHistogram->Paint("lfsame");
else fHistogram->Paint(opt.Data());
if(graph) {
if(!fGraph) {
Double_t *xx = new Double_t[fNp];
Double_t *yy = new Double_t[fNp];
for(i=0; i<fNp; ++i)
GetKnot(i,xx[i],yy[i]);
fGraph=new TGraph(fNp,xx,yy);
delete [] xx;
delete [] yy;
}
fGraph->SetMarkerColor(GetMarkerColor());
fGraph->SetMarkerStyle(GetMarkerStyle());
fGraph->SetMarkerSize(GetMarkerSize());
fGraph->Paint("p");
}
}
void TSpline::Streamer(TBuffer &R__b)
{
if (R__b.IsReading()) {
UInt_t R__s, R__c;
Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
if (R__v > 1) {
TSpline::Class()->ReadBuffer(R__b, this, R__v, R__s, R__c);
return;
}
TNamed::Streamer(R__b);
TAttLine::Streamer(R__b);
TAttFill::Streamer(R__b);
TAttMarker::Streamer(R__b);
fNp = 0;
R__b.CheckByteCount(R__s, R__c, TSpline::IsA());
} else {
TSpline::Class()->WriteBuffer(R__b,this);
}
}
TSpline3::TSpline3(const char *title,
Double_t x[], Double_t y[], Int_t n, const char *opt,
Double_t valbeg, Double_t valend) :
TSpline(title,-1,x[0],x[n-1],n,kFALSE),
fValBeg(valbeg), fValEnd(valend), fBegCond(0), fEndCond(0)
{
fName="Spline3";
if(opt) SetCond(opt);
fPoly = new TSplinePoly3[n];
for (Int_t i=0; i<n; ++i) {
fPoly[i].X() = x[i];
fPoly[i].Y() = y[i];
}
BuildCoeff();
}
TSpline3::TSpline3(const char *title,
Double_t xmin, Double_t xmax,
Double_t y[], Int_t n, const char *opt,
Double_t valbeg, Double_t valend) :
TSpline(title,(xmax-xmin)/(n-1), xmin, xmax, n, kTRUE),
fValBeg(valbeg), fValEnd(valend),
fBegCond(0), fEndCond(0)
{
fName="Spline3";
if(opt) SetCond(opt);
fPoly = new TSplinePoly3[n];
for (Int_t i=0; i<n; ++i) {
fPoly[i].X() = fXmin+i*fDelta;
fPoly[i].Y() = y[i];
}
BuildCoeff();
}
TSpline3::TSpline3(const char *title,
Double_t x[], const TF1 *func, Int_t n, const char *opt,
Double_t valbeg, Double_t valend) :
TSpline(title,-1, x[0], x[n-1], n, kFALSE),
fValBeg(valbeg), fValEnd(valend),
fBegCond(0), fEndCond(0)
{
fName="Spline3";
if(opt) SetCond(opt);
fPoly = new TSplinePoly3[n];
for (Int_t i=0; i<n; ++i) {
fPoly[i].X() = x[i];
fPoly[i].Y() = ((TF1*)func)->Eval(x[i]);
}
BuildCoeff();
}
TSpline3::TSpline3(const char *title,
Double_t xmin, Double_t xmax,
const TF1 *func, Int_t n, const char *opt,
Double_t valbeg, Double_t valend) :
TSpline(title,(xmax-xmin)/(n-1), xmin, xmax, n, kTRUE),
fValBeg(valbeg), fValEnd(valend),
fBegCond(0), fEndCond(0)
{
fName="Spline3";
if(opt) SetCond(opt);
fPoly = new TSplinePoly3[n];
for (Int_t i=0; i<n; ++i) {
Double_t x=fXmin+i*fDelta;
fPoly[i].X() = x;
fPoly[i].Y() = ((TF1*)func)->Eval(x);
}
BuildCoeff();
}
TSpline3::TSpline3(const char *title,
const TGraph *g, const char *opt,
Double_t valbeg, Double_t valend) :
TSpline(title,-1,0,0,g->GetN(),kFALSE),
fValBeg(valbeg), fValEnd(valend),
fBegCond(0), fEndCond(0)
{
fName="Spline3";
if(opt) SetCond(opt);
fPoly = new TSplinePoly3[fNp];
for (Int_t i=0; i<fNp; ++i) {
Double_t xx, yy;
g->GetPoint(i,xx,yy);
fPoly[i].X()=xx;
fPoly[i].Y()=yy;
}
fXmin = fPoly[0].X();
fXmax = fPoly[fNp-1].X();
BuildCoeff();
}
TSpline3::TSpline3(const TSpline3& sp3) :
TSpline(sp3),
fPoly(sp3.fPoly),
fValBeg(sp3.fValBeg),
fValEnd(sp3.fValEnd),
fBegCond(sp3.fBegCond),
fEndCond(sp3.fEndCond)
{
}
TSpline3& TSpline3::operator=(const TSpline3& sp3)
{
if(this!=&sp3) {
TSpline::operator=(sp3);
fPoly=sp3.fPoly;
fValBeg=sp3.fValBeg;
fValEnd=sp3.fValEnd;
fBegCond=sp3.fBegCond;
fEndCond=sp3.fEndCond;
}
return *this;
}
void TSpline3::SetCond(const char *opt)
{
const char *b1 = strstr(opt,"b1");
const char *e1 = strstr(opt,"e1");
const char *b2 = strstr(opt,"b2");
const char *e2 = strstr(opt,"e2");
if (b1 && b2)
Error("SetCond","Cannot specify first and second derivative at first point");
if (e1 && e2)
Error("SetCond","Cannot specify first and second derivative at last point");
if (b1) fBegCond=1;
else if (b2) fBegCond=2;
if (e1) fEndCond=1;
else if (e2) fEndCond=2;
}
void TSpline3::Test()
{
Double_t hx;
Double_t diff[3];
Double_t a[800], c[4];
Int_t i, j, k, m, n;
Double_t x[200], y[200], z;
Int_t jj, mm;
Int_t mm1, nm1;
Double_t com[3];
printf("1 TEST OF TSpline3 WITH NONEQUIDISTANT KNOTS\n");
n = 5;
x[0] = -3;
x[1] = -1;
x[2] = 0;
x[3] = 3;
x[4] = 4;
y[0] = 7;
y[1] = 11;
y[2] = 26;
y[3] = 56;
y[4] = 29;
m = 2;
mm = m << 1;
mm1 = mm-1;
printf("\n-N = %3d M =%2d\n",n,m);
TSpline3 *spline = new TSpline3("Test",x,y,n);
for (i = 0; i < n; ++i)
spline->GetCoeff(i,hx, a[i],a[i+200],a[i+400],a[i+600]);
delete spline;
for (i = 0; i < mm1; ++i) diff[i] = com[i] = 0;
for (k = 0; k < n; ++k) {
for (i = 0; i < mm; ++i) c[i] = a[k+i*200];
printf(" ---------------------------------------%3d --------------------------------------------\n",k+1);
printf("%12.8f\n",x[k]);
if (k == n-1) {
printf("%16.8f\n",c[0]);
} else {
for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
printf("\n");
for (i = 0; i < mm1; ++i)
if ((z=TMath::Abs(a[k+i*200])) > com[i]) com[i] = z;
z = x[k+1]-x[k];
for (i = 1; i < mm; ++i)
for (jj = i; jj < mm; ++jj) {
j = mm+i-jj;
c[j-2] = c[j-1]*z+c[j-2];
}
for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
printf("\n");
for (i = 0; i < mm1; ++i)
if (!(k >= n-2 && i != 0))
if((z = TMath::Abs(c[i]-a[k+1+i*200]))
> diff[i]) diff[i] = z;
}
}
printf(" MAXIMUM ABSOLUTE VALUES OF DIFFERENCES \n");
for (i = 0; i < mm1; ++i) printf("%18.9E",diff[i]);
printf("\n");
printf(" MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS \n");
if (TMath::Abs(c[0]) > com[0])
com[0] = TMath::Abs(c[0]);
for (i = 0; i < mm1; ++i) printf("%16.8f",com[i]);
printf("\n");
m = 2;
for (n = 10; n <= 100; n += 10) {
mm = m << 1;
mm1 = mm-1;
nm1 = n-1;
for (i = 0; i < nm1; i += 2) {
x[i] = i+1;
x[i+1] = i+2;
y[i] = 1;
y[i+1] = 0;
}
if (n % 2 != 0) {
x[n-1] = n;
y[n-1] = 1;
}
printf("\n-N = %3d M =%2d\n",n,m);
spline = new TSpline3("Test",x,y,n);
for (i = 0; i < n; ++i)
spline->GetCoeff(i,hx,a[i],a[i+200],a[i+400],a[i+600]);
delete spline;
for (i = 0; i < mm1; ++i)
diff[i] = com[i] = 0;
for (k = 0; k < n; ++k) {
for (i = 0; i < mm; ++i)
c[i] = a[k+i*200];
if (n < 11) {
printf(" ---------------------------------------%3d --------------------------------------------\n",k+1);
printf("%12.8f\n",x[k]);
if (k == n-1) printf("%16.8f\n",c[0]);
}
if (k == n-1) break;
if (n <= 10) {
for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
printf("\n");
}
for (i = 0; i < mm1; ++i)
if ((z=TMath::Abs(a[k+i*200])) > com[i])
com[i] = z;
z = x[k+1]-x[k];
for (i = 1; i < mm; ++i)
for (jj = i; jj < mm; ++jj) {
j = mm+i-jj;
c[j-2] = c[j-1]*z+c[j-2];
}
if (n <= 10) {
for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
printf("\n");
}
for (i = 0; i < mm1; ++i)
if (!(k >= n-2 && i != 0))
if ((z = TMath::Abs(c[i]-a[k+1+i*200]))
> diff[i]) diff[i] = z;
}
printf(" MAXIMUM ABSOLUTE VALUES OF DIFFERENCES \n");
for (i = 0; i < mm1; ++i) printf("%18.9E",diff[i]);
printf("\n");
printf(" MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS \n");
if (TMath::Abs(c[0]) > com[0])
com[0] = TMath::Abs(c[0]);
for (i = 0; i < mm1; ++i) printf("%16.8E",com[i]);
printf("\n");
}
}
Int_t TSpline3::FindX(Double_t x) const
{
Int_t klow=0;
if(x<=fXmin) klow=0;
else if(x>=fXmax) klow=fNp-1;
else {
if(fKstep) {
klow = TMath::Min(Int_t((x-fXmin)/fDelta),fNp-1);
} else {
Int_t khig=fNp-1, khalf;
while(khig-klow>1)
if(x>fPoly[khalf=(klow+khig)/2].X())
klow=khalf;
else
khig=khalf;
}
if(!(fPoly[klow].X()<=x && x<=fPoly[klow+1].X()))
Error("Eval",
"Binary search failed x(%d) = %f < %f < x(%d) = %f\n",
klow,fPoly[klow].X(),x,fPoly[klow+1].X());
}
return klow;
}
Double_t TSpline3::Eval(Double_t x) const
{
Int_t klow=FindX(x);
return fPoly[klow].Eval(x);
}
Double_t TSpline3::Derivative(Double_t x) const
{
Int_t klow=FindX(x);
return fPoly[klow].Derivative(x);
}
void TSpline3::SaveAs(const char *filename, Option_t * ) const
{
ofstream *f = new ofstream(filename,ios::out);
if (f == 0 || gSystem->AccessPathName(filename,kWritePermission)) {
Error("SaveAs","Cannot open file:%s\n",filename);
return;
}
char buffer[512];
Int_t nch = strlen(filename);
sprintf(buffer,"double %s",filename);
char *dot = strstr(buffer,".");
if (dot) *dot = 0;
strcat(buffer,"(double x) {\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," const int fNp = %d, fKstep = %d;\n",fNp,fKstep);
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," const double fDelta = %g, fXmin = %g, fXmax = %g;\n",fDelta,fXmin,fXmax);
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," const double fX[%d] = {",fNp);
nch = strlen(buffer); f->write(buffer,nch);
buffer[0] = 0;
Int_t i;
char numb[20];
for (i=0;i<fNp;i++) {
sprintf(numb," %g,",fPoly[i].X());
nch = strlen(numb);
if (i == fNp-1) numb[nch-1]=0;
strcat(buffer,numb);
if (i%5 == 4 || i == fNp-1) {
nch = strlen(buffer); f->write(buffer,nch);
if (i != fNp-1) sprintf(buffer,"\n ");
}
}
sprintf(buffer," };\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," const double fY[%d] = {",fNp);
nch = strlen(buffer); f->write(buffer,nch);
buffer[0] = 0;
for (i=0;i<fNp;i++) {
sprintf(numb," %g,",fPoly[i].Y());
nch = strlen(numb);
if (i == fNp-1) numb[nch-1]=0;
strcat(buffer,numb);
if (i%5 == 4 || i == fNp-1) {
nch = strlen(buffer); f->write(buffer,nch);
if (i != fNp-1) sprintf(buffer,"\n ");
}
}
sprintf(buffer," };\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," const double fB[%d] = {",fNp);
nch = strlen(buffer); f->write(buffer,nch);
buffer[0] = 0;
for (i=0;i<fNp;i++) {
sprintf(numb," %g,",fPoly[i].B());
nch = strlen(numb);
if (i == fNp-1) numb[nch-1]=0;
strcat(buffer,numb);
if (i%5 == 4 || i == fNp-1) {
nch = strlen(buffer); f->write(buffer,nch);
if (i != fNp-1) sprintf(buffer,"\n ");
}
}
sprintf(buffer," };\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," const double fC[%d] = {",fNp);
nch = strlen(buffer); f->write(buffer,nch);
buffer[0] = 0;
for (i=0;i<fNp;i++) {
sprintf(numb," %g,",fPoly[i].C());
nch = strlen(numb);
if (i == fNp-1) numb[nch-1]=0;
strcat(buffer,numb);
if (i%5 == 4 || i == fNp-1) {
nch = strlen(buffer); f->write(buffer,nch);
if (i != fNp-1) sprintf(buffer,"\n ");
}
}
sprintf(buffer," };\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," const double fD[%d] = {",fNp);
nch = strlen(buffer); f->write(buffer,nch);
buffer[0] = 0;
for (i=0;i<fNp;i++) {
sprintf(numb," %g,",fPoly[i].D());
nch = strlen(numb);
if (i == fNp-1) numb[nch-1]=0;
strcat(buffer,numb);
if (i%5 == 4 || i == fNp-1) {
nch = strlen(buffer); f->write(buffer,nch);
if (i != fNp-1) sprintf(buffer,"\n ");
}
}
sprintf(buffer," };\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," int klow=0;\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," // If out of boundaries, extrapolate. It may be badly wrong\n");
sprintf(buffer," if(x<=fXmin) klow=0;\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," else if(x>=fXmax) klow=fNp-1;\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," else {\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," if(fKstep) {\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," // Equidistant knots, use histogramming\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," klow = int((x-fXmin)/fDelta);\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," if (klow < fNp-1) klow = fNp-1;\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," } else {\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," int khig=fNp-1, khalf;\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," // Non equidistant knots, binary search\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," while(khig-klow>1)\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," if(x>fX[khalf=(klow+khig)/2]) klow=khalf;\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," else khig=khalf;\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," }\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," }\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," // Evaluate now\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," double dx=x-fX[klow];\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," return (fY[klow]+dx*(fB[klow]+dx*(fC[klow]+dx*fD[klow])));\n");
nch = strlen(buffer); f->write(buffer,nch);
f->write("}\n",2);
if (f) { f->close(); delete f;}
}
void TSpline3::BuildCoeff()
{
Int_t i, j, l, m;
Double_t divdf1,divdf3,dtau,g=0;
l = fNp-1;
for (m=1; m<fNp ; ++m) {
fPoly[m].C() = fPoly[m].X() - fPoly[m-1].X();
fPoly[m].D() = (fPoly[m].Y() - fPoly[m-1].Y())/fPoly[m].C();
}
if(fBegCond==0) {
if(fNp == 2) {
fPoly[0].D() = 1.;
fPoly[0].C() = 1.;
fPoly[0].B() = 2.*fPoly[1].D();
} else {
fPoly[0].D() = fPoly[2].C();
fPoly[0].C() = fPoly[1].C() + fPoly[2].C();
fPoly[0].B() =((fPoly[1].C()+2.*fPoly[0].C())*fPoly[1].D()*fPoly[2].C()+fPoly[1].C()*fPoly[1].C()*fPoly[2].D())/fPoly[0].C();
}
} else if (fBegCond==1) {
fPoly[0].B() = fValBeg;
fPoly[0].D() = 1.;
fPoly[0].C() = 0.;
} else if (fBegCond==2) {
fPoly[0].D() = 2.;
fPoly[0].C() = 1.;
fPoly[0].B() = 3.*fPoly[1].D() - fPoly[1].C()/2.*fValBeg;
}
if(fNp > 2) {
for (m=1; m<l; ++m) {
g = -fPoly[m+1].C()/fPoly[m-1].D();
fPoly[m].B() = g*fPoly[m-1].B() + 3.*(fPoly[m].C()*fPoly[m+1].D()+fPoly[m+1].C()*fPoly[m].D());
fPoly[m].D() = g*fPoly[m-1].C() + 2.*(fPoly[m].C() + fPoly[m+1].C());
}
if(fEndCond == 0) {
if (fNp > 3 || fBegCond != 0) {
g = fPoly[fNp-2].C() + fPoly[fNp-1].C();
fPoly[fNp-1].B() = ((fPoly[fNp-1].C()+2.*g)*fPoly[fNp-1].D()*fPoly[fNp-2].C()
+ fPoly[fNp-1].C()*fPoly[fNp-1].C()*(fPoly[fNp-2].Y()-fPoly[fNp-3].Y())/fPoly[fNp-2].C())/g;
g = -g/fPoly[fNp-2].D();
fPoly[fNp-1].D() = fPoly[fNp-2].C();
} else {
fPoly[fNp-1].B() = 2.*fPoly[fNp-1].D();
fPoly[fNp-1].D() = 1.;
g = -1./fPoly[fNp-2].D();
}
} else if (fEndCond == 1) {
fPoly[fNp-1].B() = fValEnd;
goto L30;
} else if (fEndCond == 2) {
fPoly[fNp-1].B() = 3.*fPoly[fNp-1].D() + fPoly[fNp-1].C()/2.*fValEnd;
fPoly[fNp-1].D() = 2.;
g = -1./fPoly[fNp-2].D();
}
} else {
if(fEndCond == 0) {
if (fBegCond > 0) {
fPoly[fNp-1].B() = 2.*fPoly[fNp-1].D();
fPoly[fNp-1].D() = 1.;
g = -1./fPoly[fNp-2].D();
} else {
fPoly[fNp-1].B() = fPoly[fNp-1].D();
goto L30;
}
} else if(fEndCond == 1) {
fPoly[fNp-1].B() = fValEnd;
goto L30;
} else if(fEndCond == 2) {
fPoly[fNp-1].B() = 3.*fPoly[fNp-1].D() + fPoly[fNp-1].C()/2.*fValEnd;
fPoly[fNp-1].D() = 2.;
g = -1./fPoly[fNp-2].D();
}
}
fPoly[fNp-1].D() = g*fPoly[fNp-2].C() + fPoly[fNp-1].D();
fPoly[fNp-1].B() = (g*fPoly[fNp-2].B() + fPoly[fNp-1].B())/fPoly[fNp-1].D();
L30: j = l-1;
do {
fPoly[j].B() = (fPoly[j].B() - fPoly[j].C()*fPoly[j+1].B())/fPoly[j].D();
--j;
} while (j>=0);
for (i=1; i<fNp; ++i) {
dtau = fPoly[i].C();
divdf1 = (fPoly[i].Y() - fPoly[i-1].Y())/dtau;
divdf3 = fPoly[i-1].B() + fPoly[i].B() - 2.*divdf1;
fPoly[i-1].C() = (divdf1 - fPoly[i-1].B() - divdf3)/dtau;
fPoly[i-1].D() = (divdf3/dtau)/dtau;
}
}
void TSpline3::Streamer(TBuffer &R__b)
{
if (R__b.IsReading()) {
UInt_t R__s, R__c;
Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
if (R__v > 1) {
TSpline3::Class()->ReadBuffer(R__b, this, R__v, R__s, R__c);
return;
}
TSpline::Streamer(R__b);
if (fNp > 0) {
fPoly = new TSplinePoly3[fNp];
for(Int_t i=0; i<fNp; ++i) {
fPoly[i].Streamer(R__b);
}
}
R__b >> fValBeg;
R__b >> fValEnd;
R__b >> fBegCond;
R__b >> fEndCond;
} else {
TSpline3::Class()->WriteBuffer(R__b,this);
}
}
TSpline5::TSpline5(const char *title,
Double_t x[], Double_t y[], Int_t n,
const char *opt, Double_t b1, Double_t e1,
Double_t b2, Double_t e2) :
TSpline(title,-1, x[0], x[n-1], n, kFALSE)
{
Int_t beg, end;
const char *cb1, *ce1, *cb2, *ce2;
fName="Spline5";
BoundaryConditions(opt,beg,end,cb1,ce1,cb2,ce2);
fPoly = new TSplinePoly5[fNp];
for (Int_t i=0; i<n; ++i) {
fPoly[i+beg].X() = x[i];
fPoly[i+beg].Y() = y[i];
}
SetBoundaries(b1,e1,b2,e2,cb1,ce1,cb2,ce2);
BuildCoeff();
}
TSpline5::TSpline5(const char *title,
Double_t xmin, Double_t xmax,
Double_t y[], Int_t n,
const char *opt, Double_t b1, Double_t e1,
Double_t b2, Double_t e2) :
TSpline(title,(xmax-xmin)/(n-1), xmin, xmax, n, kTRUE)
{
Int_t beg, end;
const char *cb1, *ce1, *cb2, *ce2;
fName="Spline5";
BoundaryConditions(opt,beg,end,cb1,ce1,cb2,ce2);
fPoly = new TSplinePoly5[fNp];
for (Int_t i=0; i<n; ++i) {
fPoly[i+beg].X() = fXmin+i*fDelta;
fPoly[i+beg].Y() = y[i];
}
SetBoundaries(b1,e1,b2,e2,cb1,ce1,cb2,ce2);
BuildCoeff();
}
TSpline5::TSpline5(const char *title,
Double_t x[], const TF1 *func, Int_t n,
const char *opt, Double_t b1, Double_t e1,
Double_t b2, Double_t e2) :
TSpline(title,-1, x[0], x[n-1], n, kFALSE)
{
Int_t beg, end;
const char *cb1, *ce1, *cb2, *ce2;
fName="Spline5";
BoundaryConditions(opt,beg,end,cb1,ce1,cb2,ce2);
fPoly = new TSplinePoly5[fNp];
for (Int_t i=0; i<n; i++) {
fPoly[i+beg].X() = x[i];
fPoly[i+beg].Y() = ((TF1*)func)->Eval(x[i]);
}
SetBoundaries(b1,e1,b2,e2,cb1,ce1,cb2,ce2);
BuildCoeff();
}
TSpline5::TSpline5(const char *title,
Double_t xmin, Double_t xmax,
const TF1 *func, Int_t n,
const char *opt, Double_t b1, Double_t e1,
Double_t b2, Double_t e2) :
TSpline(title,(xmax-xmin)/(n-1), xmin, xmax, n, kTRUE)
{
Int_t beg, end;
const char *cb1, *ce1, *cb2, *ce2;
fName="Spline5";
BoundaryConditions(opt,beg,end,cb1,ce1,cb2,ce2);
fPoly = new TSplinePoly5[fNp];
for (Int_t i=0; i<n; ++i) {
Double_t x=fXmin+i*fDelta;
fPoly[i+beg].X() = x;
fPoly[i+beg].Y() = ((TF1*)func)->Eval(x);
}
SetBoundaries(b1,e1,b2,e2,cb1,ce1,cb2,ce2);
BuildCoeff();
}
TSpline5::TSpline5(const char *title,
const TGraph *g,
const char *opt, Double_t b1, Double_t e1,
Double_t b2, Double_t e2) :
TSpline(title,-1,0,0,g->GetN(),kFALSE)
{
Int_t beg, end;
const char *cb1, *ce1, *cb2, *ce2;
fName="Spline5";
BoundaryConditions(opt,beg,end,cb1,ce1,cb2,ce2);
fPoly = new TSplinePoly5[fNp];
for (Int_t i=0; i<fNp-beg; ++i) {
Double_t xx, yy;
g->GetPoint(i,xx,yy);
fPoly[i+beg].X()=xx;
fPoly[i+beg].Y()=yy;
}
SetBoundaries(b1,e1,b2,e2,cb1,ce1,cb2,ce2);
fXmin = fPoly[0].X();
fXmax = fPoly[fNp-1].X();
BuildCoeff();
}
TSpline5::TSpline5(const TSpline5& sp5) :
TSpline(sp5),
fPoly(sp5.fPoly)
{
}
TSpline5& TSpline5::operator=(const TSpline5& sp5)
{
if(this!=&sp5) {
TSpline::operator=(sp5);
fPoly=sp5.fPoly;
}
return *this;
}
void TSpline5::BoundaryConditions(const char *opt,Int_t &beg,Int_t &end,
const char *&cb1,const char *&ce1,
const char *&cb2,const char *&ce2)
{
cb1=ce1=cb2=ce2=0;
beg=end=0;
if(opt) {
cb1 = strstr(opt,"b1");
ce1 = strstr(opt,"e1");
cb2 = strstr(opt,"b2");
ce2 = strstr(opt,"e2");
if(cb2) {
fNp=fNp+2;
beg=2;
} else if(cb1) {
fNp=fNp+1;
beg=1;
}
if(ce2) {
fNp=fNp+2;
end=2;
} else if(ce1) {
fNp=fNp+1;
end=1;
}
}
}
void TSpline5::SetBoundaries(Double_t b1, Double_t e1, Double_t b2, Double_t e2,
const char *cb1, const char *ce1, const char *cb2,
const char *ce2)
{
if(cb2) {
fPoly[0].X() = fPoly[1].X() = fPoly[2].X();
fPoly[0].Y() = fPoly[2].Y();
fPoly[2].Y()=b2;
if(cb1)
fPoly[1].Y()=b1;
else
fPoly[1].Y()=(fPoly[3].Y()-fPoly[0].Y())/(fPoly[3].X()-fPoly[2].X());
} else if(cb1) {
fPoly[0].X() = fPoly[1].X();
fPoly[0].Y() = fPoly[1].Y();
fPoly[1].Y()=b1;
}
if(ce2) {
fPoly[fNp-1].X() = fPoly[fNp-2].X() = fPoly[fNp-3].X();
fPoly[fNp-1].Y()=e2;
if(ce1)
fPoly[fNp-2].Y()=e1;
else
fPoly[fNp-2].Y()=
(fPoly[fNp-3].Y()-fPoly[fNp-4].Y())
/(fPoly[fNp-3].X()-fPoly[fNp-4].X());
} else if(ce1) {
fPoly[fNp-1].X() = fPoly[fNp-2].X();
fPoly[fNp-1].Y()=e1;
}
}
Int_t TSpline5::FindX(Double_t x) const
{
Int_t klow=0;
if(x<=fXmin) klow=0;
else if(x>=fXmax) klow=fNp-1;
else {
if(fKstep) {
klow = TMath::Min(Int_t((x-fXmin)/fDelta),fNp-1);
} else {
Int_t khig=fNp-1, khalf;
while(khig-klow>1)
if(x>fPoly[khalf=(klow+khig)/2].X())
klow=khalf;
else
khig=khalf;
}
if(!(fPoly[klow].X()<=x && x<=fPoly[klow+1].X()))
Error("Eval",
"Binary search failed x(%d) = %f < %f < x(%d) = %f\n",
klow,fPoly[klow].X(),x,fPoly[klow+1].X());
}
return klow;
}
Double_t TSpline5::Eval(Double_t x) const
{
Int_t klow=FindX(x);
return fPoly[klow].Eval(x);
}
Double_t TSpline5::Derivative(Double_t x) const
{
Int_t klow=FindX(x);
return fPoly[klow].Derivative(x);
}
void TSpline5::SaveAs(const char *filename, Option_t * ) const
{
ofstream *f = new ofstream(filename,ios::out);
if (f == 0 || gSystem->AccessPathName(filename,kWritePermission)) {
Error("SaveAs","Cannot open file:%s\n",filename);
return;
}
char buffer[512];
Int_t nch = strlen(filename);
sprintf(buffer,"double %s",filename);
char *dot = strstr(buffer,".");
if (dot) *dot = 0;
strcat(buffer,"(double x) {\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," const int fNp = %d, fKstep = %d;\n",fNp,fKstep);
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," const double fDelta = %g, fXmin = %g, fXmax = %g;\n",fDelta,fXmin,fXmax);
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," const double fX[%d] = {",fNp);
nch = strlen(buffer); f->write(buffer,nch);
buffer[0] = 0;
Int_t i;
char numb[20];
for (i=0;i<fNp;i++) {
sprintf(numb," %g,",fPoly[i].X());
nch = strlen(numb);
if (i == fNp-1) numb[nch-1]=0;
strcat(buffer,numb);
if (i%5 == 4 || i == fNp-1) {
nch = strlen(buffer); f->write(buffer,nch);
if (i != fNp-1) sprintf(buffer,"\n ");
}
}
sprintf(buffer," };\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," const double fY[%d] = {",fNp);
nch = strlen(buffer); f->write(buffer,nch);
buffer[0] = 0;
for (i=0;i<fNp;i++) {
sprintf(numb," %g,",fPoly[i].Y());
nch = strlen(numb);
if (i == fNp-1) numb[nch-1]=0;
strcat(buffer,numb);
if (i%5 == 4 || i == fNp-1) {
nch = strlen(buffer); f->write(buffer,nch);
if (i != fNp-1) sprintf(buffer,"\n ");
}
}
sprintf(buffer," };\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," const double fB[%d] = {",fNp);
nch = strlen(buffer); f->write(buffer,nch);
buffer[0] = 0;
for (i=0;i<fNp;i++) {
sprintf(numb," %g,",fPoly[i].B());
nch = strlen(numb);
if (i == fNp-1) numb[nch-1]=0;
strcat(buffer,numb);
if (i%5 == 4 || i == fNp-1) {
nch = strlen(buffer); f->write(buffer,nch);
if (i != fNp-1) sprintf(buffer,"\n ");
}
}
sprintf(buffer," };\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," const double fC[%d] = {",fNp);
nch = strlen(buffer); f->write(buffer,nch);
buffer[0] = 0;
for (i=0;i<fNp;i++) {
sprintf(numb," %g,",fPoly[i].C());
nch = strlen(numb);
if (i == fNp-1) numb[nch-1]=0;
strcat(buffer,numb);
if (i%5 == 4 || i == fNp-1) {
nch = strlen(buffer); f->write(buffer,nch);
if (i != fNp-1) sprintf(buffer,"\n ");
}
}
sprintf(buffer," };\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," const double fD[%d] = {",fNp);
nch = strlen(buffer); f->write(buffer,nch);
buffer[0] = 0;
for (i=0;i<fNp;i++) {
sprintf(numb," %g,",fPoly[i].D());
nch = strlen(numb);
if (i == fNp-1) numb[nch-1]=0;
strcat(buffer,numb);
if (i%5 == 4 || i == fNp-1) {
nch = strlen(buffer); f->write(buffer,nch);
if (i != fNp-1) sprintf(buffer,"\n ");
}
}
sprintf(buffer," };\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," const double fE[%d] = {",fNp);
nch = strlen(buffer); f->write(buffer,nch);
buffer[0] = 0;
for (i=0;i<fNp;i++) {
sprintf(numb," %g,",fPoly[i].E());
nch = strlen(numb);
if (i == fNp-1) numb[nch-1]=0;
strcat(buffer,numb);
if (i%5 == 4 || i == fNp-1) {
nch = strlen(buffer); f->write(buffer,nch);
if (i != fNp-1) sprintf(buffer,"\n ");
}
}
sprintf(buffer," };\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," const double fF[%d] = {",fNp);
nch = strlen(buffer); f->write(buffer,nch);
buffer[0] = 0;
for (i=0;i<fNp;i++) {
sprintf(numb," %g,",fPoly[i].F());
nch = strlen(numb);
if (i == fNp-1) numb[nch-1]=0;
strcat(buffer,numb);
if (i%5 == 4 || i == fNp-1) {
nch = strlen(buffer); f->write(buffer,nch);
if (i != fNp-1) sprintf(buffer,"\n ");
}
}
sprintf(buffer," };\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," int klow=0;\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," // If out of boundaries, extrapolate. It may be badly wrong\n");
sprintf(buffer," if(x<=fXmin) klow=0;\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," else if(x>=fXmax) klow=fNp-1;\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," else {\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," if(fKstep) {\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," // Equidistant knots, use histogramming\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," klow = int((x-fXmin)/fDelta);\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," if (klow < fNp-1) klow = fNp-1;\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," } else {\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," int khig=fNp-1, khalf;\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," // Non equidistant knots, binary search\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," while(khig-klow>1)\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," if(x>fX[khalf=(klow+khig)/2]) klow=khalf;\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," else khig=khalf;\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," }\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," }\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," // Evaluate now\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," double dx=x-fX[klow];\n");
nch = strlen(buffer); f->write(buffer,nch);
sprintf(buffer," return (fY[klow]+dx*(fB[klow]+dx*(fC[klow]+dx*(fD[klow]+dx*(fE[klow]+dx*fF[klow])))));\n");
nch = strlen(buffer); f->write(buffer,nch);
f->write("}\n",2);
if (f) { f->close(); delete f;}
}
void TSpline5::BuildCoeff()
{
Int_t i, m;
Double_t pqqr, p, q, r, s, t, u, v,
b1, p2, p3, q2, q3, r2, pq, pr, qr;
if (fNp <= 2) {
return;
}
m = fNp-2;
q = fPoly[1].X()-fPoly[0].X();
r = fPoly[2].X()-fPoly[1].X();
q2 = q*q;
r2 = r*r;
qr = q+r;
fPoly[0].D() = fPoly[0].E() = 0;
if (q) fPoly[1].D() = q*6.*q2/(qr*qr);
else fPoly[1].D() = 0;
if (m > 1) {
for (i = 1; i < m; ++i) {
p = q;
q = r;
r = fPoly[i+2].X()-fPoly[i+1].X();
p2 = q2;
q2 = r2;
r2 = r*r;
pq = qr;
qr = q+r;
if (q) {
q3 = q2*q;
pr = p*r;
pqqr = pq*qr;
fPoly[i+1].D() = q3*6./(qr*qr);
fPoly[i].D() += (q+q)*(pr*15.*pr+(p+r)*q
*(pr* 20.+q2*7.)+q2*
((p2+r2)*8.+pr*21.+q2+q2))/(pqqr*pqqr);
fPoly[i-1].D() += q3*6./(pq*pq);
fPoly[i].E() = q2*(p*qr+pq*3.*(qr+r+r))/(pqqr*qr);
fPoly[i-1].E() += q2*(r*pq+qr*3.*(pq+p+p))/(pqqr*pq);
fPoly[i-1].F() = q3/pqqr;
} else
fPoly[i+1].D() = fPoly[i].E() = fPoly[i-1].F() = 0;
}
}
if (r) fPoly[m-1].D() += r*6.*r2/(qr*qr);
for (i = 1; i < fNp; ++i) {
if (fPoly[i].X() != fPoly[i-1].X()) {
fPoly[i].B() =
(fPoly[i].Y()-fPoly[i-1].Y())/(fPoly[i].X()-fPoly[i-1].X());
} else {
fPoly[i].B() = fPoly[i].Y();
fPoly[i].Y() = fPoly[i-1].Y();
}
}
for (i = 2; i < fNp; ++i) {
if (fPoly[i].X() != fPoly[i-2].X()) {
fPoly[i].C() =
(fPoly[i].B()-fPoly[i-1].B())/(fPoly[i].X()-fPoly[i-2].X());
} else {
fPoly[i].C() = fPoly[i].B()*.5;
fPoly[i].B() = fPoly[i-1].B();
}
}
if (m > 1) {
p=fPoly[0].C()=fPoly[m-1].E()=fPoly[0].F()
=fPoly[m-2].F()=fPoly[m-1].F()=0;
fPoly[1].C() = fPoly[3].C()-fPoly[2].C();
fPoly[1].D() = 1./fPoly[1].D();
if (m > 2) {
for (i = 2; i < m; ++i) {
q = fPoly[i-1].D()*fPoly[i-1].E();
fPoly[i].D() = 1./(fPoly[i].D()-p*fPoly[i-2].F()-q*fPoly[i-1].E());
fPoly[i].E() -= q*fPoly[i-1].F();
fPoly[i].C() = fPoly[i+2].C()-fPoly[i+1].C()-p*fPoly[i-2].C()
-q*fPoly[i-1].C();
p = fPoly[i-1].D()*fPoly[i-1].F();
}
}
}
fPoly[fNp-2].C() = fPoly[fNp-1].C() = 0;
if (fNp > 3)
for (i=fNp-3; i > 0; --i)
fPoly[i].C() = (fPoly[i].C()-fPoly[i].E()*fPoly[i+1].C()
-fPoly[i].F()*fPoly[i+2].C())*fPoly[i].D();
m = fNp-1;
q = fPoly[1].X()-fPoly[0].X();
r = fPoly[2].X()-fPoly[1].X();
b1 = fPoly[1].B();
q3 = q*q*q;
qr = q+r;
if (qr) {
v = fPoly[1].C()/qr;
t = v;
} else
v = t = 0;
if (q) fPoly[0].F() = v/q;
else fPoly[0].F() = 0;
for (i = 1; i < m; ++i) {
p = q;
q = r;
if (i != m-1) r = fPoly[i+2].X()-fPoly[i+1].X();
else r = 0;
p3 = q3;
q3 = q*q*q;
pq = qr;
qr = q+r;
s = t;
if (qr) t = (fPoly[i+1].C()-fPoly[i].C())/qr;
else t = 0;
u = v;
v = t-s;
if (pq) {
fPoly[i].F() = fPoly[i-1].F();
if (q) fPoly[i].F() = v/q;
fPoly[i].E() = s*5.;
fPoly[i].D() = (fPoly[i].C()-q*s)*10;
fPoly[i].C() =
fPoly[i].D()*(p-q)+(fPoly[i+1].B()-fPoly[i].B()+(u-fPoly[i].E())*
p3-(v+fPoly[i].E())*q3)/pq;
fPoly[i].B() = (p*(fPoly[i+1].B()-v*q3)+q*(fPoly[i].B()-u*p3))/pq-p
*q*(fPoly[i].D()+fPoly[i].E()*(q-p));
} else {
fPoly[i].C() = fPoly[i-1].C();
fPoly[i].D() = fPoly[i].E() = fPoly[i].F() = 0;
}
}
p = fPoly[1].X()-fPoly[0].X();
s = fPoly[0].F()*p*p*p;
fPoly[0].E() = fPoly[0].D() = 0;
fPoly[0].C() = fPoly[1].C()-s*10;
fPoly[0].B() = b1-(fPoly[0].C()+s)*p;
q = fPoly[fNp-1].X()-fPoly[fNp-2].X();
t = fPoly[fNp-2].F()*q*q*q;
fPoly[fNp-1].E() = fPoly[fNp-1].D() = 0;
fPoly[fNp-1].C() = fPoly[fNp-2].C()+t*10;
fPoly[fNp-1].B() += (fPoly[fNp-1].C()-t)*q;
}
void TSpline5::Test()
{
Double_t hx;
Double_t diff[5];
Double_t a[1200], c[6];
Int_t i, j, k, m, n;
Double_t p, x[200], y[200], z;
Int_t jj, mm, nn;
Int_t mm1, nm1;
Double_t com[5];
printf("1 TEST OF TSpline5 WITH NONEQUIDISTANT KNOTS\n");
n = 5;
x[0] = -3;
x[1] = -1;
x[2] = 0;
x[3] = 3;
x[4] = 4;
y[0] = 7;
y[1] = 11;
y[2] = 26;
y[3] = 56;
y[4] = 29;
m = 3;
mm = m << 1;
mm1 = mm-1;
printf("\n-N = %3d M =%2d\n",n,m);
TSpline5 *spline = new TSpline5("Test",x,y,n);
for (i = 0; i < n; ++i)
spline->GetCoeff(i,hx, a[i],a[i+200],a[i+400],
a[i+600],a[i+800],a[i+1000]);
delete spline;
for (i = 0; i < mm1; ++i) diff[i] = com[i] = 0;
for (k = 0; k < n; ++k) {
for (i = 0; i < mm; ++i) c[i] = a[k+i*200];
printf(" ---------------------------------------%3d --------------------------------------------\n",k+1);
printf("%12.8f\n",x[k]);
if (k == n-1) {
printf("%16.8f\n",c[0]);
} else {
for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
printf("\n");
for (i = 0; i < mm1; ++i)
if ((z=TMath::Abs(a[k+i*200])) > com[i]) com[i] = z;
z = x[k+1]-x[k];
for (i = 1; i < mm; ++i)
for (jj = i; jj < mm; ++jj) {
j = mm+i-jj;
c[j-2] = c[j-1]*z+c[j-2];
}
for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
printf("\n");
for (i = 0; i < mm1; ++i)
if (!(k >= n-2 && i != 0))
if((z = TMath::Abs(c[i]-a[k+1+i*200]))
> diff[i]) diff[i] = z;
}
}
printf(" MAXIMUM ABSOLUTE VALUES OF DIFFERENCES \n");
for (i = 0; i < mm1; ++i) printf("%18.9E",diff[i]);
printf("\n");
printf(" MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS \n");
if (TMath::Abs(c[0]) > com[0])
com[0] = TMath::Abs(c[0]);
for (i = 0; i < mm1; ++i) printf("%16.8f",com[i]);
printf("\n");
m = 3;
for (n = 10; n <= 100; n += 10) {
mm = m << 1;
mm1 = mm-1;
nm1 = n-1;
for (i = 0; i < nm1; i += 2) {
x[i] = i+1;
x[i+1] = i+2;
y[i] = 1;
y[i+1] = 0;
}
if (n % 2 != 0) {
x[n-1] = n;
y[n-1] = 1;
}
printf("\n-N = %3d M =%2d\n",n,m);
spline = new TSpline5("Test",x,y,n);
for (i = 0; i < n; ++i)
spline->GetCoeff(i,hx,a[i],a[i+200],a[i+400],
a[i+600],a[i+800],a[i+1000]);
delete spline;
for (i = 0; i < mm1; ++i)
diff[i] = com[i] = 0;
for (k = 0; k < n; ++k) {
for (i = 0; i < mm; ++i)
c[i] = a[k+i*200];
if (n < 11) {
printf(" ---------------------------------------%3d --------------------------------------------\n",k+1);
printf("%12.8f\n",x[k]);
if (k == n-1) printf("%16.8f\n",c[0]);
}
if (k == n-1) break;
if (n <= 10) {
for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
printf("\n");
}
for (i = 0; i < mm1; ++i)
if ((z=TMath::Abs(a[k+i*200])) > com[i])
com[i] = z;
z = x[k+1]-x[k];
for (i = 1; i < mm; ++i)
for (jj = i; jj < mm; ++jj) {
j = mm+i-jj;
c[j-2] = c[j-1]*z+c[j-2];
}
if (n <= 10) {
for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
printf("\n");
}
for (i = 0; i < mm1; ++i)
if (!(k >= n-2 && i != 0))
if ((z = TMath::Abs(c[i]-a[k+1+i*200]))
> diff[i]) diff[i] = z;
}
printf(" MAXIMUM ABSOLUTE VALUES OF DIFFERENCES \n");
for (i = 0; i < mm1; ++i) printf("%18.9E",diff[i]);
printf("\n");
printf(" MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS \n");
if (TMath::Abs(c[0]) > com[0])
com[0] = TMath::Abs(c[0]);
for (i = 0; i < mm1; ++i) printf("%16.8E",com[i]);
printf("\n");
}
printf("1 TEST OF TSpline5 WITH NONEQUIDISTANT DOUBLE KNOTS\n");
n = 5;
x[0] = -3;
x[1] = -3;
x[2] = -1;
x[3] = -1;
x[4] = 0;
x[5] = 0;
x[6] = 3;
x[7] = 3;
x[8] = 4;
x[9] = 4;
y[0] = 7;
y[1] = 2;
y[2] = 11;
y[3] = 15;
y[4] = 26;
y[5] = 10;
y[6] = 56;
y[7] = -27;
y[8] = 29;
y[9] = -30;
m = 3;
nn = n << 1;
mm = m << 1;
mm1 = mm-1;
printf("-N = %3d M =%2d\n",n,m);
spline = new TSpline5("Test",x,y,nn);
for (i = 0; i < nn; ++i)
spline->GetCoeff(i,hx,a[i],a[i+200],a[i+400],
a[i+600],a[i+800],a[i+1000]);
delete spline;
for (i = 0; i < mm1; ++i)
diff[i] = com[i] = 0;
for (k = 0; k < nn; ++k) {
for (i = 0; i < mm; ++i)
c[i] = a[k+i*200];
printf(" ---------------------------------------%3d --------------------------------------------\n",k+1);
printf("%12.8f\n",x[k]);
if (k == nn-1) {
printf("%16.8f\n",c[0]);
break;
}
for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
printf("\n");
for (i = 0; i < mm1; ++i)
if ((z=TMath::Abs(a[k+i*200])) > com[i]) com[i] = z;
z = x[k+1]-x[k];
for (i = 1; i < mm; ++i)
for (jj = i; jj < mm; ++jj) {
j = mm+i-jj;
c[j-2] = c[j-1]*z+c[j-2];
}
for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
printf("\n");
for (i = 0; i < mm1; ++i)
if (!(k >= nn-2 && i != 0))
if ((z = TMath::Abs(c[i]-a[k+1+i*200]))
> diff[i]) diff[i] = z;
}
printf(" MAXIMUM ABSOLUTE VALUES OF DIFFERENCES \n");
for (i = 1; i <= mm1; ++i) {
printf("%18.9E",diff[i-1]);
}
printf("\n");
if (TMath::Abs(c[0]) > com[0])
com[0] = TMath::Abs(c[0]);
printf(" MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS \n");
for (i = 0; i < mm1; ++i) printf("%16.8f",com[i]);
printf("\n");
m = 3;
for (n = 10; n <= 100; n += 10) {
nn = n << 1;
mm = m << 1;
mm1 = mm-1;
p = 0;
for (i = 0; i < n; ++i) {
p += TMath::Abs(TMath::Sin(i+1));
x[(i << 1)] = p;
x[(i << 1)+1] = p;
y[(i << 1)] = TMath::Cos(i+1)-.5;
y[(i << 1)+1] = TMath::Cos((i << 1)+2)-.5;
}
printf("-N = %3d M =%2d\n",n,m);
spline = new TSpline5("Test",x,y,nn);
for (i = 0; i < nn; ++i)
spline->GetCoeff(i,hx,a[i],a[i+200],a[i+400],
a[i+600],a[i+800],a[i+1000]);
delete spline;
for (i = 0; i < mm1; ++i)
diff[i] = com[i] = 0;
for (k = 0; k < nn; ++k) {
for (i = 0; i < mm; ++i)
c[i] = a[k+i*200];
if (n < 11) {
printf(" ---------------------------------------%3d --------------------------------------------\n",k+1);
printf("%12.8f\n",x[k]);
if (k == nn-1) printf("%16.8f\n",c[0]);
}
if (k == nn-1) break;
if (n <= 10) {
for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
printf("\n");
}
for (i = 0; i < mm1; ++i)
if ((z=TMath::Abs(a[k+i*200])) > com[i]) com[i] = z;
z = x[k+1]-x[k];
for (i = 1; i < mm; ++i) {
for (jj = i; jj < mm; ++jj) {
j = mm+i-jj;
c[j-2] = c[j-1]*z+c[j-2];
}
}
if (n <= 10) {
for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
printf("\n");
}
for (i = 0; i < mm1; ++i)
if (!(k >= nn-2 && i != 0))
if ((z = TMath::Abs(c[i]-a[k+1+i*200]))
> diff[i]) diff[i] = z;
}
printf(" MAXIMUM ABSOLUTE VALUES OF DIFFERENCES \n");
for (i = 0; i < mm1; ++i) printf("%18.9E",diff[i]);
printf("\n");
printf(" MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS \n");
if (TMath::Abs(c[0]) > com[0])
com[0] = TMath::Abs(c[0]);
for (i = 0; i < mm1; ++i) printf("%18.9E",com[i]);
printf("\n");
}
printf("1 TEST OF TSpline5 WITH NONEQUIDISTANT KNOTS,\n");
printf(" ONE DOUBLE, ONE TRIPLE KNOT\n");
n = 8;
x[0] = -3;
x[1] = -1;
x[2] = -1;
x[3] = 0;
x[4] = 3;
x[5] = 3;
x[6] = 3;
x[7] = 4;
y[0] = 7;
y[1] = 11;
y[2] = 15;
y[3] = 26;
y[4] = 56;
y[5] = -30;
y[6] = -7;
y[7] = 29;
m = 3;
mm = m << 1;
mm1 = mm-1;
printf("-N = %3d M =%2d\n",n,m);
spline=new TSpline5("Test",x,y,n);
for (i = 0; i < n; ++i)
spline->GetCoeff(i,hx,a[i],a[i+200],a[i+400],
a[i+600],a[i+800],a[i+1000]);
delete spline;
for (i = 0; i < mm1; ++i)
diff[i] = com[i] = 0;
for (k = 0; k < n; ++k) {
for (i = 0; i < mm; ++i)
c[i] = a[k+i*200];
printf(" ---------------------------------------%3d --------------------------------------------\n",k+1);
printf("%12.8f\n",x[k]);
if (k == n-1) {
printf("%16.8f\n",c[0]);
break;
}
for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
printf("\n");
for (i = 0; i < mm1; ++i)
if ((z=TMath::Abs(a[k+i*200])) > com[i]) com[i] = z;
z = x[k+1]-x[k];
for (i = 1; i < mm; ++i)
for (jj = i; jj < mm; ++jj) {
j = mm+i-jj;
c[j-2] = c[j-1]*z+c[j-2];
}
for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
printf("\n");
for (i = 0; i < mm1; ++i)
if (!(k >= n-2 && i != 0))
if ((z = TMath::Abs(c[i]-a[k+1+i*200]))
> diff[i]) diff[i] = z;
}
printf(" MAXIMUM ABSOLUTE VALUES OF DIFFERENCES \n");
for (i = 0; i < mm1; ++i) printf("%18.9E",diff[i]);
printf("\n");
printf(" MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS \n");
if (TMath::Abs(c[0]) > com[0])
com[0] = TMath::Abs(c[0]);
for (i = 0; i < mm1; ++i) printf("%16.8f",com[i]);
printf("\n");
printf("1 TEST OF TSpline5 WITH NONEQUIDISTANT KNOTS,\n");
printf(" TWO DOUBLE, ONE TRIPLE KNOT\n");
n = 10;
x[0] = 0;
x[1] = 2;
x[2] = 2;
x[3] = 3;
x[4] = 3;
x[5] = 3;
x[6] = 5;
x[7] = 8;
x[8] = 9;
x[9] = 9;
y[0] = 163;
y[1] = 237;
y[2] = -127;
y[3] = 119;
y[4] = -65;
y[5] = 192;
y[6] = 293;
y[7] = 326;
y[8] = 0;
y[9] = -414;
m = 3;
mm = m << 1;
mm1 = mm-1;
printf("-N = %3d M =%2d\n",n,m);
spline = new TSpline5("Test",x,y,n);
for (i = 0; i < n; ++i)
spline->GetCoeff(i,hx,a[i],a[i+200],a[i+400],
a[i+600],a[i+800],a[i+1000]);
delete spline;
for (i = 0; i < mm1; ++i)
diff[i] = com[i] = 0;
for (k = 0; k < n; ++k) {
for (i = 0; i < mm; ++i)
c[i] = a[k+i*200];
printf(" ---------------------------------------%3d --------------------------------------------\n",k+1);
printf("%12.8f\n",x[k]);
if (k == n-1) {
printf("%16.8f\n",c[0]);
break;
}
for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
printf("\n");
for (i = 0; i < mm1; ++i)
if ((z=TMath::Abs(a[k+i*200])) > com[i]) com[i] = z;
z = x[k+1]-x[k];
for (i = 1; i < mm; ++i)
for (jj = i; jj < mm; ++jj) {
j = mm+i-jj;
c[j-2] = c[j-1]*z+c[j-2];
}
for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
printf("\n");
for (i = 0; i < mm1; ++i)
if (!(k >= n-2 && i != 0))
if((z = TMath::Abs(c[i]-a[k+1+i*200]))
> diff[i]) diff[i] = z;
}
printf(" MAXIMUM ABSOLUTE VALUES OF DIFFERENCES \n");
for (i = 0; i < mm1; ++i) printf("%18.9E",diff[i]);
printf("\n");
printf(" MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS \n");
if (TMath::Abs(c[0]) > com[0])
com[0] = TMath::Abs(c[0]);
for (i = 0; i < mm1; ++i) printf("%16.8f",com[i]);
printf("\n");
}
void TSpline5::Streamer(TBuffer &R__b)
{
if (R__b.IsReading()) {
UInt_t R__s, R__c;
Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
if (R__v > 1) {
TSpline5::Class()->ReadBuffer(R__b, this, R__v, R__s, R__c);
return;
}
TSpline::Streamer(R__b);
if (fNp > 0) {
fPoly = new TSplinePoly5[fNp];
for(Int_t i=0; i<fNp; ++i) {
fPoly[i].Streamer(R__b);
}
}
} else {
TSpline5::Class()->WriteBuffer(R__b,this);
}
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.