Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TSpline.cxx
Go to the documentation of this file.
1// @(#)root/hist:$Id$
2// Author: Federico Carminati 28/02/2000
3
4/*************************************************************************
5 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12/** \class TSpline
13 \ingroup Hist
14 Base class for spline implementation containing the Draw/Paint methods.
15*/
16
17#include "TROOT.h"
18#include "TGraph.h"
19#include "TBuffer.h"
20#include "TSpline.h"
21#include "TVirtualPad.h"
22#include "TH1.h"
23#include "TF1.h"
24#include "TSystem.h"
25#include "TMath.h"
26#include "strlcpy.h"
27#include "snprintf.h"
28#include <iostream>
29#include <fstream>
30
37
38////////////////////////////////////////////////////////////////////////////////
39/// Copy constructor.
40
42 TNamed(sp),
43 TAttLine(sp),
44 TAttFill(sp),
46 fDelta(sp.fDelta),
47 fXmin(sp.fXmin),
48 fXmax(sp.fXmax),
49 fNp(sp.fNp),
50 fKstep(sp.fKstep),
51 fHistogram(nullptr),
52 fGraph(nullptr),
53 fNpx(sp.fNpx)
54{
55}
56
57////////////////////////////////////////////////////////////////////////////////
58/// Destructor.
59
61{
62 if(fHistogram) delete fHistogram;
63 if(fGraph) delete fGraph;
64}
65
66////////////////////////////////////////////////////////////////////////////////
67/// Assignment operator.
68
70{
71 if(this!=&sp) {
73 TAttLine::operator=(sp);
74 TAttFill::operator=(sp);
75 TAttMarker::operator=(sp);
76 fDelta=sp.fDelta;
77 fXmin=sp.fXmin;
78 fXmax=sp.fXmax;
79 fNp=sp.fNp;
80 fKstep=sp.fKstep;
81 fHistogram=nullptr;
82 fGraph=nullptr;
83 fNpx=sp.fNpx;
84 }
85 return *this;
86}
87
88////////////////////////////////////////////////////////////////////////////////
89/// Draw this function with its current attributes.
90///
91/// Possible option values are:
92///
93/// - "SAME" superimpose on top of existing picture
94/// - "L" connect all computed points with a straight line
95/// - "C" connect all computed points with a smooth curve.
96/// - "P" add a polymarker at each knot
97///
98/// Note that the default value is "L". Therefore to draw on top
99/// of an existing picture, specify option "LSAME"
100
102{
103 TString opt = option;
104 opt.ToLower();
105 if (gPad && !opt.Contains("same")) gPad->Clear();
106
108}
109
110////////////////////////////////////////////////////////////////////////////////
111/// Compute distance from point px,py to a spline.
112
114{
115 if (!fHistogram) return 999;
116 return fHistogram->DistancetoPrimitive(px, py);
117}
118
119////////////////////////////////////////////////////////////////////////////////
120/// Execute action corresponding to one event.
121
123{
124 if (!fHistogram) return;
125 fHistogram->ExecuteEvent(event, px, py);
126}
127
128////////////////////////////////////////////////////////////////////////////////
129/// Paint this function with its current attributes.
130
132{
133 Int_t i;
134 Double_t xv;
135
136 TString opt = option;
137 opt.ToLower();
139 pmin = gPad->PadtoX(gPad->GetUxmin());
140 pmax = gPad->PadtoX(gPad->GetUxmax());
141 xmin = fXmin;
142 xmax = fXmax;
143 if (opt.Contains("same")) {
144 if (xmax < pmin) return; // Otto: completely outside
145 if (xmin > pmax) return;
146 if (xmin < pmin) xmin = pmin;
147 if (xmax > pmax) xmax = pmax;
148 } else {
149 gPad->Clear();
150 }
151
152 // Create a temporary histogram and fill each channel with the function value
153 if (fHistogram)
154 if ((!gPad->GetLogx() && fHistogram->TestBit(TH1::kLogX)) ||
155 (gPad->GetLogx() && !fHistogram->TestBit(TH1::kLogX)))
156 { delete fHistogram; fHistogram = nullptr;}
157
158 if (fHistogram) {
159 //if (xmin != fXmin || xmax != fXmax)
161 } else {
162 // if logx, we must bin in logx and not in x !!!
163 // otherwise if several decades, one gets crazy results
164 if (xmin > 0 && gPad->GetLogx()) {
165 Double_t *xbins = new Double_t[fNpx+1];
169 for (i=0;i<=fNpx;i++) {
170 xbins[i] = gPad->PadtoX(xlogmin+ i*dlogx);
171 }
172 fHistogram = new TH1F("Spline",GetTitle(),fNpx,xbins);
174 delete [] xbins;
175 } else {
176 fHistogram = new TH1F("Spline",GetTitle(),fNpx,xmin,xmax);
177 }
178 if (!fHistogram) return;
179 fHistogram->SetDirectory(nullptr);
180 }
181 for (i=1;i<=fNpx;i++) {
183 fHistogram->SetBinContent(i,this->Eval(xv));
184 }
185
186 // Copy Function attributes to histogram attributes
196
197 // Draw the histogram
198 // but first strip off the 'p' option if any
199 char *o = (char *) opt.Data();
200 Int_t j=0;
201 i=0;
202 Bool_t graph=kFALSE;
203 do
204 if(o[i]=='p') graph=kTRUE ; else o[j++]=o[i];
205 while(o[i++]);
206 if (opt.Length() == 0 ) fHistogram->Paint("lf");
207 else if (opt == "same") fHistogram->Paint("lfsame");
208 else fHistogram->Paint(opt.Data());
209
210 // Think about the graph, if demanded
211 if(graph) {
212 if(!fGraph) {
213 Double_t *xx = new Double_t[fNp];
214 Double_t *yy = new Double_t[fNp];
215 for(i=0; i<fNp; ++i)
216 GetKnot(i,xx[i],yy[i]);
217 fGraph=new TGraph(fNp,xx,yy);
218 delete [] xx;
219 delete [] yy;
220 }
224 fGraph->Paint("p");
225 }
226}
227
228////////////////////////////////////////////////////////////////////////////////
229/// Stream an object of class TSpline.
230
232{
233 if (R__b.IsReading()) {
235 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
236 if (R__v > 1) {
237 R__b.ReadClassBuffer(TSpline::Class(), this, R__v, R__s, R__c);
238 return;
239 }
240 //====process old versions before automatic schema evolution
245
246 fNp = 0;
247 /*
248 R__b >> fDelta;
249 R__b >> fXmin;
250 R__b >> fXmax;
251 R__b >> fNp;
252 R__b >> fKstep;
253 R__b >> fHistogram;
254 R__b >> fGraph;
255 R__b >> fNpx;
256 */
257 R__b.CheckByteCount(R__s, R__c, TSpline::IsA());
258 //====end of old versions
259
260 } else {
261 R__b.WriteClassBuffer(TSpline::Class(),this);
262 }
263}
264
265/** \class TSplinePoly
266 \ingroup Hist
267 Base class for TSpline knot.
268*/
269
270////////////////////////////////////////////////////////////////////////////////
271/// Assignment operator.
272
274{
275 if(this != &other) {
278 }
279 return *this;
280}
281
282////////////////////////////////////////////////////////////////////////////////
283/// Utility called by the copy constructors and = operator.
284
286{
287 fX = other.fX;
288 fY = other.fY;
289}
290
291/** \class TSplinePoly3
292 \ingroup Hist
293 Class for TSpline3 knot.
294*/
295
296////////////////////////////////////////////////////////////////////////////////
297/// Assignment operator.
298
300{
301 if(this != &other) {
304 }
305 return *this;
306}
307
308////////////////////////////////////////////////////////////////////////////////
309/// Utility called by the copy constructors and = operator.
310
312{
313 fB = other.fB;
314 fC = other.fC;
315 fD = other.fD;
316}
317
318/** \class TSplinePoly5
319 \ingroup Hist
320 Class for TSpline5 knot.
321*/
322
323////////////////////////////////////////////////////////////////////////////////
324/// Assignment operator.
325
327{
328 if(this != &other) {
331 }
332 return *this;
333}
334
335////////////////////////////////////////////////////////////////////////////////
336/// Utility called by the copy constructors and = operator.
337
339{
340 fB = other.fB;
341 fC = other.fC;
342 fD = other.fD;
343 fE = other.fE;
344 fF = other.fF;
345}
346
347/** \class TSpline3
348 \ingroup Hist
349 Class to create third splines to interpolate knots
350 Arbitrary conditions can be introduced for first and second
351 derivatives at beginning and ending points
352 */
353
354////////////////////////////////////////////////////////////////////////////////
355/// Third spline creator given an array of arbitrary knots in increasing
356/// abscissa order and possibly end point conditions.
357
358TSpline3::TSpline3(const char *title,
359 Double_t x[], Double_t y[], Int_t n, const char *opt,
361 TSpline(title,-1,x[0],x[n-1],n,kFALSE),
362 fValBeg(valbeg), fValEnd(valend), fBegCond(0), fEndCond(0)
363{
364 fName="Spline3";
365
366 // Set endpoint conditions
367 if(opt) SetCond(opt);
368
369 // Create the polynomial terms and fill
370 // them with node information
371 fPoly = new TSplinePoly3[n];
372 for (Int_t i=0; i<n; ++i) {
373 fPoly[i].X() = x[i];
374 fPoly[i].Y() = y[i];
375 }
376
377 // Build the spline coefficients
378 BuildCoeff();
379}
380
381////////////////////////////////////////////////////////////////////////////////
382/// Third spline creator given an array of
383/// arbitrary function values on equidistant n abscissa
384/// values from xmin to xmax and possibly end point conditions.
385
386TSpline3::TSpline3(const char *title,
388 Double_t y[], Int_t n, const char *opt,
390 TSpline(title,(xmax-xmin)/(n-1), xmin, xmax, n, kTRUE),
391 fValBeg(valbeg), fValEnd(valend),
392 fBegCond(0), fEndCond(0)
393{
394 fName="Spline3";
395
396 // Set endpoint conditions
397 if(opt) SetCond(opt);
398
399 // Create the polynomial terms and fill
400 // them with node information
401 fPoly = new TSplinePoly3[n];
402 for (Int_t i=0; i<n; ++i) {
403 fPoly[i].X() = fXmin+i*fDelta;
404 fPoly[i].Y() = y[i];
405 }
406
407 // Build the spline coefficients
408 BuildCoeff();
409}
410
411////////////////////////////////////////////////////////////////////////////////
412/// Third spline creator given an array of
413/// arbitrary abscissas in increasing order and a function
414/// to interpolate and possibly end point conditions.
415
416TSpline3::TSpline3(const char *title,
417 Double_t x[], const TF1 *func, Int_t n, const char *opt,
419 TSpline(title,-1, x[0], x[n-1], n, kFALSE),
420 fValBeg(valbeg), fValEnd(valend),
421 fBegCond(0), fEndCond(0)
422{
423 fName="Spline3";
424
425 // Set endpoint conditions
426 if(opt) SetCond(opt);
427
428 // Create the polynomial terms and fill
429 // them with node information
430 fPoly = new TSplinePoly3[n];
431 for (Int_t i=0; i<n; ++i) {
432 fPoly[i].X() = x[i];
433 fPoly[i].Y() = ((TF1*)func)->Eval(x[i]);
434 }
435
436 // Build the spline coefficients
437 BuildCoeff();
438}
439
440////////////////////////////////////////////////////////////////////////////////
441/// Third spline creator given a function to be
442/// evaluated on n equidistant abscissa points between xmin
443/// and xmax and possibly end point conditions.
444
445TSpline3::TSpline3(const char *title,
447 const TF1 *func, Int_t n, const char *opt,
449 TSpline(title,(xmax-xmin)/(n-1), xmin, xmax, n, kTRUE),
450 fValBeg(valbeg), fValEnd(valend),
451 fBegCond(0), fEndCond(0)
452{
453 fName="Spline3";
454
455 // Set endpoint conditions
456 if(opt) SetCond(opt);
457
458 // Create the polynomial terms and fill
459 // them with node information
460 fPoly = new TSplinePoly3[n];
461 //when func is null we return. In this case it is assumed that the spline
462 //points will be given later via SetPoint and SetPointCoeff
463 if (!func) {fKstep = kFALSE; fDelta = -1; return;}
464 for (Int_t i=0; i<n; ++i) {
466 fPoly[i].X() = x;
467 fPoly[i].Y() = ((TF1*)func)->Eval(x);
468 }
469
470 // Build the spline coefficients
471 BuildCoeff();
472}
473
474////////////////////////////////////////////////////////////////////////////////
475/// Third spline creator given a TGraph with
476/// abscissa in increasing order and possibly end
477/// point conditions.
478
479TSpline3::TSpline3(const char *title,
480 const TGraph *g, const char *opt,
482 TSpline(title,-1,0,0,g->GetN(),kFALSE),
483 fValBeg(valbeg), fValEnd(valend),
484 fBegCond(0), fEndCond(0)
485{
486 fName="Spline3";
487
488 // Set endpoint conditions
489 if(opt) SetCond(opt);
490
491 // Create the polynomial terms and fill
492 // them with node information
493 fPoly = new TSplinePoly3[fNp];
494 for (Int_t i=0; i<fNp; ++i) {
495 Double_t xx, yy;
496 g->GetPoint(i,xx,yy);
497 fPoly[i].X()=xx;
498 fPoly[i].Y()=yy;
499 }
500 fXmin = fPoly[0].X();
501 fXmax = fPoly[fNp-1].X();
502
503 // Build the spline coefficients
504 BuildCoeff();
505}
506
507////////////////////////////////////////////////////////////////////////////////
508/// Third spline creator given a TH1.
509
510TSpline3::TSpline3(const TH1 *h, const char *opt,
512 TSpline(h->GetTitle(),-1,0,0,h->GetNbinsX(),kFALSE),
513 fValBeg(valbeg), fValEnd(valend),
514 fBegCond(0), fEndCond(0)
515{
516 fName=h->GetName();
517
518 // Set endpoint conditions
519 if(opt) SetCond(opt);
520
521 // Create the polynomial terms and fill
522 // them with node information
523 fPoly = new TSplinePoly3[fNp];
524 for (Int_t i=0; i<fNp; ++i) {
525 fPoly[i].X()=h->GetXaxis()->GetBinCenter(i+1);
526 fPoly[i].Y()=h->GetBinContent(i+1);
527 }
528 fXmin = fPoly[0].X();
529 fXmax = fPoly[fNp-1].X();
530
531 // Build the spline coefficients
532 BuildCoeff();
533}
534
535////////////////////////////////////////////////////////////////////////////////
536/// Copy constructor.
537
539 : TSpline(sp3), fValBeg(sp3.fValBeg), fValEnd(sp3.fValEnd), fBegCond(sp3.fBegCond), fEndCond(sp3.fEndCond)
540{
541 if (fNp > 0)
542 fPoly = new TSplinePoly3[fNp];
543 for (Int_t i = 0; i < fNp; ++i)
544 fPoly[i] = sp3.fPoly[i];
545}
546
547////////////////////////////////////////////////////////////////////////////////
548/// Assignment operator.
549
551{
552 if(this!=&sp3) {
554 if (fPoly) {
555 delete[] fPoly;
556 fPoly = nullptr;
557 }
558 if (fNp > 0)
559 fPoly = new TSplinePoly3[fNp];
560 for (Int_t i=0; i<fNp; ++i)
561 fPoly[i] = sp3.fPoly[i];
562
563 fValBeg=sp3.fValBeg;
564 fValEnd=sp3.fValEnd;
565 fBegCond=sp3.fBegCond;
566 fEndCond=sp3.fEndCond;
567 }
568 return *this;
569}
570
571////////////////////////////////////////////////////////////////////////////////
572/// Check the boundary conditions.
573
574void TSpline3::SetCond(const char *opt)
575{
576 const char *b1 = strstr(opt,"b1");
577 const char *e1 = strstr(opt,"e1");
578 const char *b2 = strstr(opt,"b2");
579 const char *e2 = strstr(opt,"e2");
580 if (b1 && b2)
581 Error("SetCond","Cannot specify first and second derivative at first point");
582 if (e1 && e2)
583 Error("SetCond","Cannot specify first and second derivative at last point");
584 if (b1) fBegCond=1;
585 else if (b2) fBegCond=2;
586 if (e1) fEndCond=1;
587 else if (e2) fEndCond=2;
588}
589
590////////////////////////////////////////////////////////////////////////////////
591/// Test method for TSpline5
592///
593/// ~~~ {.cpp}
594/// n number of data points.
595/// m 2*m-1 is order of spline.
596/// m = 2 always for third spline.
597/// nn,nm1,mm,
598/// mm1,i,k,
599/// j,jj temporary integer variables.
600/// z,p temporary double precision variables.
601/// x[n] the sequence of knots.
602/// y[n] the prescribed function values at the knots.
603/// a[200][4] two dimensional array whose columns are
604/// the computed spline coefficients
605/// diff[3] maximum values of differences of values and
606/// derivatives to right and left of knots.
607/// com[3] maximum values of coefficients.
608/// ~~~
609///
610/// test of TSpline3 with non equidistant knots and
611/// equidistant knots follows.
612
614{
615 Double_t hx;
616 Double_t diff[3];
617 Double_t a[800], c[4];
618 Int_t i, j, k, m, n;
619 Double_t x[200], y[200], z;
620 Int_t jj, mm;
621 Int_t mm1, nm1;
622 Double_t com[3];
623 printf("1 TEST OF TSpline3 WITH NONEQUIDISTANT KNOTS\n");
624 n = 5;
625 x[0] = -3;
626 x[1] = -1;
627 x[2] = 0;
628 x[3] = 3;
629 x[4] = 4;
630 y[0] = 7;
631 y[1] = 11;
632 y[2] = 26;
633 y[3] = 56;
634 y[4] = 29;
635 m = 2;
636 mm = m << 1;
637 mm1 = mm-1;
638 printf("\n-N = %3d M =%2d\n",n,m);
639 TSpline3 *spline = new TSpline3("Test",x,y,n);
640 for (i = 0; i < n; ++i)
641 spline->GetCoeff(i,hx, a[i],a[i+200],a[i+400],a[i+600]);
642 delete spline;
643 for (i = 0; i < mm1; ++i) diff[i] = com[i] = 0;
644 for (k = 0; k < n; ++k) {
645 for (i = 0; i < mm; ++i) c[i] = a[k+i*200];
646 printf(" ---------------------------------------%3d --------------------------------------------\n",k+1);
647 printf("%12.8f\n",x[k]);
648 if (k == n-1) {
649 printf("%16.8f\n",c[0]);
650 } else {
651 for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
652 printf("\n");
653 for (i = 0; i < mm1; ++i)
654 if ((z=TMath::Abs(a[k+i*200])) > com[i]) com[i] = z;
655 z = x[k+1]-x[k];
656 for (i = 1; i < mm; ++i)
657 for (jj = i; jj < mm; ++jj) {
658 j = mm+i-jj;
659 c[j-2] = c[j-1]*z+c[j-2];
660 }
661 for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
662 printf("\n");
663 for (i = 0; i < mm1; ++i)
664 if (!(k >= n-2 && i != 0))
665 if((z = TMath::Abs(c[i]-a[k+1+i*200]))
666 > diff[i]) diff[i] = z;
667 }
668 }
669 printf(" MAXIMUM ABSOLUTE VALUES OF DIFFERENCES \n");
670 for (i = 0; i < mm1; ++i) printf("%18.9E",diff[i]);
671 printf("\n");
672 printf(" MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS \n");
673 if (TMath::Abs(c[0]) > com[0])
674 com[0] = TMath::Abs(c[0]);
675 for (i = 0; i < mm1; ++i) printf("%16.8f",com[i]);
676 printf("\n");
677 m = 2;
678 for (n = 10; n <= 100; n += 10) {
679 mm = m << 1;
680 mm1 = mm-1;
681 nm1 = n-1;
682 for (i = 0; i < nm1; i += 2) {
683 x[i] = i+1;
684 x[i+1] = i+2;
685 y[i] = 1;
686 y[i+1] = 0;
687 }
688 if (n % 2 != 0) {
689 x[n-1] = n;
690 y[n-1] = 1;
691 }
692 printf("\n-N = %3d M =%2d\n",n,m);
693 spline = new TSpline3("Test",x,y,n);
694 for (i = 0; i < n; ++i)
695 spline->GetCoeff(i,hx,a[i],a[i+200],a[i+400],a[i+600]);
696 delete spline;
697 for (i = 0; i < mm1; ++i)
698 diff[i] = com[i] = 0;
699 for (k = 0; k < n; ++k) {
700 for (i = 0; i < mm; ++i)
701 c[i] = a[k+i*200];
702 if (n < 11) {
703 printf(" ---------------------------------------%3d --------------------------------------------\n",k+1);
704 printf("%12.8f\n",x[k]);
705 if (k == n-1) printf("%16.8f\n",c[0]);
706 }
707 if (k == n-1) break;
708 if (n <= 10) {
709 for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
710 printf("\n");
711 }
712 for (i = 0; i < mm1; ++i)
713 if ((z=TMath::Abs(a[k+i*200])) > com[i])
714 com[i] = z;
715 z = x[k+1]-x[k];
716 for (i = 1; i < mm; ++i)
717 for (jj = i; jj < mm; ++jj) {
718 j = mm+i-jj;
719 c[j-2] = c[j-1]*z+c[j-2];
720 }
721 if (n <= 10) {
722 for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
723 printf("\n");
724 }
725 for (i = 0; i < mm1; ++i)
726 if (!(k >= n-2 && i != 0))
727 if ((z = TMath::Abs(c[i]-a[k+1+i*200]))
728 > diff[i]) diff[i] = z;
729 }
730 printf(" MAXIMUM ABSOLUTE VALUES OF DIFFERENCES \n");
731 for (i = 0; i < mm1; ++i) printf("%18.9E",diff[i]);
732 printf("\n");
733 printf(" MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS \n");
734 if (TMath::Abs(c[0]) > com[0])
735 com[0] = TMath::Abs(c[0]);
736 for (i = 0; i < mm1; ++i) printf("%16.8E",com[i]);
737 printf("\n");
738 }
739}
740
741////////////////////////////////////////////////////////////////////////////////
742/// Find X.
743
745{
746 Int_t klow=0, khig=fNp-1;
747 //
748 // If out of boundaries, extrapolate
749 // It may be badly wrong
750 if(x<=fXmin) klow=0;
751 else if(x>=fXmax) klow=khig;
752 else {
753 if(fKstep) {
754 //
755 // Equidistant knots, use histogramming
757 // Correction for rounding errors
758 if (x < fPoly[klow].X())
759 klow = TMath::Max(klow-1,0);
760 else if (klow < khig) {
761 if (x > fPoly[klow+1].X()) ++klow;
762 }
763 } else {
764 Int_t khalf;
765 //
766 // Non equidistant knots, binary search
767 while(khig-klow>1)
768 if(x>fPoly[khalf=(klow+khig)/2].X())
769 klow=khalf;
770 else
771 khig=khalf;
772 //
773 // This could be removed, sanity check
774 if(!(fPoly[klow].X()<=x && x<=fPoly[klow+1].X()))
775 Error("Eval",
776 "Binary search failed x(%d) = %f < x= %f < x(%d) = %f\n",
777 klow,fPoly[klow].X(),x,klow+1,fPoly[klow+1].X());
778 }
779 }
780 return klow;
781}
782
783////////////////////////////////////////////////////////////////////////////////
784/// Eval this spline at x.
785
787{
788 Int_t klow=FindX(x);
789 if (klow >= fNp-1 && fNp > 1) klow = fNp-2; //see: https://savannah.cern.ch/bugs/?71651
790 return fPoly[klow].Eval(x);
791}
792
793////////////////////////////////////////////////////////////////////////////////
794/// Derivative.
795
797{
798 Int_t klow=FindX(x);
799 if (klow >= fNp-1) klow = fNp-2; //see: https://savannah.cern.ch/bugs/?71651
800 return fPoly[klow].Derivative(x);
801}
802
803////////////////////////////////////////////////////////////////////////////////
804/// Write this spline as a C++ function that can be executed without ROOT
805/// the name of the function is the name of the file up to the "." if any.
806
807void TSpline3::SaveAs(const char *filename, Option_t * /*option*/) const
808{
809 //open the file
810 std::ofstream *f = new std::ofstream(filename,std::ios::out);
811 if (f == nullptr || gSystem->AccessPathName(filename,kWritePermission)) {
812 Error("SaveAs","Cannot open file:%s\n",filename);
813 return;
814 }
815
816 //write the function name and the spline constants
817 char buffer[512];
819 snprintf(buffer,512,"double %s",filename);
820 char *dot = strstr(buffer,".");
821 if (dot) *dot = 0;
822 strlcat(buffer,"(double x) {\n",512);
823 nch = strlen(buffer); f->write(buffer,nch);
824 snprintf(buffer,512," const int fNp = %d, fKstep = %d;\n",fNp,fKstep);
825 nch = strlen(buffer); f->write(buffer,nch);
826 snprintf(buffer,512," const double fDelta = %.17g, fXmin = %.17g, fXmax = %.17g;\n",fDelta,fXmin,fXmax);
827 nch = strlen(buffer); f->write(buffer,nch);
828
829 //write the spline coefficients
830 //array fX
831 snprintf(buffer,512," const double fX[%d] = {",fNp);
832 nch = strlen(buffer); f->write(buffer,nch);
833 buffer[0] = 0;
834 Int_t i;
835 char numb[30];
836 for (i=0;i<fNp;i++) {
837 snprintf(numb,30," %.17g,",fPoly[i].X());
838 nch = strlen(numb);
839 if (i == fNp-1) numb[nch-1]=0;
840 strlcat(buffer,numb,512);
841 if (i%5 == 4 || i == fNp-1) {
842 nch = strlen(buffer); f->write(buffer,nch);
843 if (i != fNp-1) snprintf(buffer,512,"\n ");
844 }
845 }
846 snprintf(buffer,512," };\n");
847 nch = strlen(buffer); f->write(buffer,nch);
848 //array fY
849 snprintf(buffer,512," const double fY[%d] = {",fNp);
850 nch = strlen(buffer); f->write(buffer,nch);
851 buffer[0] = 0;
852 for (i=0;i<fNp;i++) {
853 snprintf(numb,30," %.17g,",fPoly[i].Y());
854 nch = strlen(numb);
855 if (i == fNp-1) numb[nch-1]=0;
856 strlcat(buffer,numb,512);
857 if (i%5 == 4 || i == fNp-1) {
858 nch = strlen(buffer); f->write(buffer,nch);
859 if (i != fNp-1) snprintf(buffer,512,"\n ");
860 }
861 }
862 snprintf(buffer,512," };\n");
863 nch = strlen(buffer); f->write(buffer,nch);
864 //array fB
865 snprintf(buffer,512," const double fB[%d] = {",fNp);
866 nch = strlen(buffer); f->write(buffer,nch);
867 buffer[0] = 0;
868 for (i=0;i<fNp;i++) {
869 snprintf(numb,30," %.17g,",fPoly[i].B());
870 nch = strlen(numb);
871 if (i == fNp-1) numb[nch-1]=0;
872 strlcat(buffer,numb,512);
873 if (i%5 == 4 || i == fNp-1) {
874 nch = strlen(buffer); f->write(buffer,nch);
875 if (i != fNp-1) snprintf(buffer,512,"\n ");
876 }
877 }
878 snprintf(buffer,512," };\n");
879 nch = strlen(buffer); f->write(buffer,nch);
880 //array fC
881 snprintf(buffer,512," const double fC[%d] = {",fNp);
882 nch = strlen(buffer); f->write(buffer,nch);
883 buffer[0] = 0;
884 for (i=0;i<fNp;i++) {
885 snprintf(numb,30," %.17g,",fPoly[i].C());
886 nch = strlen(numb);
887 if (i == fNp-1) numb[nch-1]=0;
888 strlcat(buffer,numb,512);
889 if (i%5 == 4 || i == fNp-1) {
890 nch = strlen(buffer); f->write(buffer,nch);
891 if (i != fNp-1) snprintf(buffer,512,"\n ");
892 }
893 }
894 snprintf(buffer,512," };\n");
895 nch = strlen(buffer); f->write(buffer,nch);
896 //array fD
897 snprintf(buffer,512," const double fD[%d] = {",fNp);
898 nch = strlen(buffer); f->write(buffer,nch);
899 buffer[0] = 0;
900 for (i=0;i<fNp;i++) {
901 snprintf(numb,30," %.17g,",fPoly[i].D());
902 nch = strlen(numb);
903 if (i == fNp-1) numb[nch-1]=0;
904 strlcat(buffer,numb,512);
905 if (i%5 == 4 || i == fNp-1) {
906 nch = strlen(buffer); f->write(buffer,nch);
907 if (i != fNp-1) snprintf(buffer,512,"\n ");
908 }
909 }
910 snprintf(buffer,512," };\n");
911 nch = strlen(buffer); f->write(buffer,nch);
912
913 //generate code for the spline evaluation
914 snprintf(buffer,512," int klow=0;\n");
915 nch = strlen(buffer); f->write(buffer,nch);
916
917 snprintf(buffer,512," // If out of boundaries, extrapolate. It may be badly wrong\n");
918 snprintf(buffer,512," if(x<=fXmin) klow=0;\n");
919 nch = strlen(buffer); f->write(buffer,nch);
920 snprintf(buffer,512," else if(x>=fXmax) klow=fNp-1;\n");
921 nch = strlen(buffer); f->write(buffer,nch);
922 snprintf(buffer,512," else {\n");
923 nch = strlen(buffer); f->write(buffer,nch);
924 snprintf(buffer,512," if(fKstep) {\n");
925 nch = strlen(buffer); f->write(buffer,nch);
926
927 snprintf(buffer,512," // Equidistant knots, use histogramming\n");
928 nch = strlen(buffer); f->write(buffer,nch);
929 snprintf(buffer,512," klow = int((x-fXmin)/fDelta);\n");
930 nch = strlen(buffer); f->write(buffer,nch);
931 snprintf(buffer,512," if (klow > fNp-1) klow = fNp-1;\n");
932 nch = strlen(buffer); f->write(buffer,nch);
933 snprintf(buffer,512," } else {\n");
934 nch = strlen(buffer); f->write(buffer,nch);
935 snprintf(buffer,512," int khig=fNp-1, khalf;\n");
936 nch = strlen(buffer); f->write(buffer,nch);
937
938 snprintf(buffer,512," // Non equidistant knots, binary search\n");
939 nch = strlen(buffer); f->write(buffer,nch);
940 snprintf(buffer,512," while(khig-klow>1)\n");
941 nch = strlen(buffer); f->write(buffer,nch);
942 snprintf(buffer,512," if(x>fX[khalf=(klow+khig)/2]) klow=khalf;\n");
943 nch = strlen(buffer); f->write(buffer,nch);
944 snprintf(buffer,512," else khig=khalf;\n");
945 nch = strlen(buffer); f->write(buffer,nch);
946 snprintf(buffer,512," }\n");
947 nch = strlen(buffer); f->write(buffer,nch);
948 snprintf(buffer,512," }\n");
949 nch = strlen(buffer); f->write(buffer,nch);
950 snprintf(buffer,512," // Evaluate now\n");
951 nch = strlen(buffer); f->write(buffer,nch);
952 snprintf(buffer,512," double dx=x-fX[klow];\n");
953 nch = strlen(buffer); f->write(buffer,nch);
954 snprintf(buffer,512," return (fY[klow]+dx*(fB[klow]+dx*(fC[klow]+dx*fD[klow])));\n");
955 nch = strlen(buffer); f->write(buffer,nch);
956
957 //close file
958 f->write("}\n",2);
959
960 if (f) { f->close(); delete f;}
961}
962
963////////////////////////////////////////////////////////////////////////////////
964/// Save primitive as a C++ statement(s) on output stream out.
965
966void TSpline3::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
967{
968 SavePrimitiveConstructor(out, Class(), "spline3",
969 TString::Format("\"%s\", %g, %g, (TF1 *)nullptr, %d, \"\", %g, %g",
970 TString(GetTitle()).ReplaceSpecialCppChars().Data(), fXmin, fXmax, fNp,
971 fValBeg, fValEnd));
972
973 out << " spline3->SetName(\"" << TString(GetName()).ReplaceSpecialCppChars() << "\");\n";
974
975 SaveFillAttributes(out, "spline3", 0, 1001);
976 SaveLineAttributes(out, "spline3", 1, 1, 1);
977 SaveMarkerAttributes(out, "spline3", 1, 1, 1);
978 if (fNpx != 100)
979 out << " spline3->SetNpx(" << fNpx << ");\n";
980
981 for (Int_t i = 0; i < fNp; i++) {
982 out << " spline3->SetPoint(" << i << "," << fPoly[i].X() << "," << fPoly[i].Y() << ");\n";
983 out << " spline3->SetPointCoeff(" << i << "," << fPoly[i].B() << "," << fPoly[i].C() << "," << fPoly[i].D()
984 << ");\n";
985 }
986 out << " spline3->Draw(\"" << TString(option).ReplaceSpecialCppChars() << "\");\n";
987}
988
989////////////////////////////////////////////////////////////////////////////////
990/// Set point number i.
991
993{
994 if (i < 0 || i >= fNp) return;
995 fPoly[i].X()= x;
996 fPoly[i].Y()= y;
997}
998
999////////////////////////////////////////////////////////////////////////////////
1000/// Set point coefficient number i.
1001
1003{
1004 if (i < 0 || i >= fNp) return;
1005 fPoly[i].B()= b;
1006 fPoly[i].C()= c;
1007 fPoly[i].D()= d;
1008}
1009
1010////////////////////////////////////////////////////////////////////////////////
1011/// Build coefficients.
1012///
1013/// ~~~ {.cpp}
1014/// subroutine cubspl ( tau, c, n, ibcbeg, ibcend )
1015/// from * a practical guide to splines * by c. de boor
1016/// ************************ input ***************************
1017/// n = number of data points. assumed to be .ge. 2.
1018/// (tau(i), c(1,i), i=1,...,n) = abscissae and ordinates of the
1019/// data points. tau is assumed to be strictly increasing.
1020/// ibcbeg, ibcend = boundary condition indicators, and
1021/// c(2,1), c(2,n) = boundary condition information. specifically,
1022/// ibcbeg = 0 means no boundary condition at tau(1) is given.
1023/// in this case, the not-a-knot condition is used, i.e. the
1024/// jump in the third derivative across tau(2) is forced to
1025/// zero, thus the first and the second cubic polynomial pieces
1026/// are made to coincide.)
1027/// ibcbeg = 1 means that the slope at tau(1) is made to equal
1028/// c(2,1), supplied by input.
1029/// ibcbeg = 2 means that the second derivative at tau(1) is
1030/// made to equal c(2,1), supplied by input.
1031/// ibcend = 0, 1, or 2 has analogous meaning concerning the
1032/// boundary condition at tau(n), with the additional infor-
1033/// mation taken from c(2,n).
1034/// *********************** output **************************
1035/// c(j,i), j=1,...,4; i=1,...,l (= n-1) = the polynomial coefficients
1036/// of the cubic interpolating spline with interior knots (or
1037/// joints) tau(2), ..., tau(n-1). precisely, in the interval
1038/// (tau(i), tau(i+1)), the spline f is given by
1039/// f(x) = c(1,i)+h*(c(2,i)+h*(c(3,i)+h*c(4,i)/3.)/2.)
1040/// where h = x - tau(i). the function program *ppvalu* may be
1041/// used to evaluate f or its derivatives from tau,c, l = n-1,
1042/// and k=4.
1043/// ~~~
1044
1046{
1047 Int_t i, j, l, m;
1049 //***** a tridiagonal linear system for the unknown slopes s(i) of
1050 // f at tau(i), i=1,...,n, is generated and then solved by gauss elim-
1051 // ination, with s(i) ending up in c(2,i), all i.
1052 // c(3,.) and c(4,.) are used initially for temporary storage.
1053 l = fNp-1;
1054 // compute first differences of x sequence and store in C also,
1055 // compute first divided difference of data and store in D.
1056 for (m=1; m<fNp ; ++m) {
1057 fPoly[m].C() = fPoly[m].X() - fPoly[m-1].X();
1058 fPoly[m].D() = (fPoly[m].Y() - fPoly[m-1].Y())/fPoly[m].C();
1059 }
1060 // construct first equation from the boundary condition, of the form
1061 // D[0]*s[0] + C[0]*s[1] = B[0]
1062 if(fBegCond==0) {
1063 if(fNp == 2) {
1064 // no condition at left end and n = 2.
1065 fPoly[0].D() = 1.;
1066 fPoly[0].C() = 1.;
1067 fPoly[0].B() = 2.*fPoly[1].D();
1068 } else {
1069 // not-a-knot condition at left end and n .gt. 2.
1070 fPoly[0].D() = fPoly[2].C();
1071 fPoly[0].C() = fPoly[1].C() + fPoly[2].C();
1072 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();
1073 }
1074 } else if (fBegCond==1) {
1075 // slope prescribed at left end.
1076 fPoly[0].B() = fValBeg;
1077 fPoly[0].D() = 1.;
1078 fPoly[0].C() = 0.;
1079 } else if (fBegCond==2) {
1080 // second derivative prescribed at left end.
1081 fPoly[0].D() = 2.;
1082 fPoly[0].C() = 1.;
1083 fPoly[0].B() = 3.*fPoly[1].D() - fPoly[1].C()/2.*fValBeg;
1084 }
1085 if(fNp > 2) {
1086 // if there are interior knots, generate the corresp. equations and car-
1087 // ry out the forward pass of gauss elimination, after which the m-th
1088 // equation reads D[m]*s[m] + C[m]*s[m+1] = B[m].
1089 for (m=1; m<l; ++m) {
1090 g = -fPoly[m+1].C()/fPoly[m-1].D();
1091 fPoly[m].B() = g*fPoly[m-1].B() + 3.*(fPoly[m].C()*fPoly[m+1].D()+fPoly[m+1].C()*fPoly[m].D());
1092 fPoly[m].D() = g*fPoly[m-1].C() + 2.*(fPoly[m].C() + fPoly[m+1].C());
1093 }
1094 // construct last equation from the second boundary condition, of the form
1095 // (-g*D[n-2])*s[n-2] + D[n-1]*s[n-1] = B[n-1]
1096 // if slope is prescribed at right end, one can go directly to back-
1097 // substitution, since c array happens to be set up just right for it
1098 // at this point.
1099 if(fEndCond == 0) {
1100 if (fNp > 3 || fBegCond != 0) {
1101 // not-a-knot and n .ge. 3, and either n.gt.3 or also not-a-knot at
1102 // left end point.
1103 g = fPoly[fNp-2].C() + fPoly[fNp-1].C();
1104 fPoly[fNp-1].B() = ((fPoly[fNp-1].C()+2.*g)*fPoly[fNp-1].D()*fPoly[fNp-2].C()
1105 + fPoly[fNp-1].C()*fPoly[fNp-1].C()*(fPoly[fNp-2].Y()-fPoly[fNp-3].Y())/fPoly[fNp-2].C())/g;
1106 g = -g/fPoly[fNp-2].D();
1107 fPoly[fNp-1].D() = fPoly[fNp-2].C();
1108 } else {
1109 // either (n=3 and not-a-knot also at left) or (n=2 and not not-a-
1110 // knot at left end point).
1111 fPoly[fNp-1].B() = 2.*fPoly[fNp-1].D();
1112 fPoly[fNp-1].D() = 1.;
1113 g = -1./fPoly[fNp-2].D();
1114 }
1115 } else if (fEndCond == 1) {
1116 fPoly[fNp-1].B() = fValEnd;
1117 goto L30;
1118 } else if (fEndCond == 2) {
1119 // second derivative prescribed at right endpoint.
1120 fPoly[fNp-1].B() = 3.*fPoly[fNp-1].D() + fPoly[fNp-1].C()/2.*fValEnd;
1121 fPoly[fNp-1].D() = 2.;
1122 g = -1./fPoly[fNp-2].D();
1123 }
1124 } else {
1125 if(fEndCond == 0) {
1126 if (fBegCond > 0) {
1127 // either (n=3 and not-a-knot also at left) or (n=2 and not not-a-
1128 // knot at left end point).
1129 fPoly[fNp-1].B() = 2.*fPoly[fNp-1].D();
1130 fPoly[fNp-1].D() = 1.;
1131 g = -1./fPoly[fNp-2].D();
1132 } else {
1133 // not-a-knot at right endpoint and at left endpoint and n = 2.
1134 fPoly[fNp-1].B() = fPoly[fNp-1].D();
1135 goto L30;
1136 }
1137 } else if(fEndCond == 1) {
1138 fPoly[fNp-1].B() = fValEnd;
1139 goto L30;
1140 } else if(fEndCond == 2) {
1141 // second derivative prescribed at right endpoint.
1142 fPoly[fNp-1].B() = 3.*fPoly[fNp-1].D() + fPoly[fNp-1].C()/2.*fValEnd;
1143 fPoly[fNp-1].D() = 2.;
1144 g = -1./fPoly[fNp-2].D();
1145 }
1146 }
1147 // complete forward pass of gauss elimination.
1148 fPoly[fNp-1].D() = g*fPoly[fNp-2].C() + fPoly[fNp-1].D();
1149 fPoly[fNp-1].B() = (g*fPoly[fNp-2].B() + fPoly[fNp-1].B())/fPoly[fNp-1].D();
1150 // carry out back substitution
1151L30: j = l-1;
1152 do {
1153 fPoly[j].B() = (fPoly[j].B() - fPoly[j].C()*fPoly[j+1].B())/fPoly[j].D();
1154 --j;
1155 } while (j>=0);
1156 //****** generate cubic coefficients in each interval, i.e., the deriv.s
1157 // at its left endpoint, from value and slope at its endpoints.
1158 for (i=1; i<fNp; ++i) {
1159 dtau = fPoly[i].C();
1160 divdf1 = (fPoly[i].Y() - fPoly[i-1].Y())/dtau;
1161 divdf3 = fPoly[i-1].B() + fPoly[i].B() - 2.*divdf1;
1162 fPoly[i-1].C() = (divdf1 - fPoly[i-1].B() - divdf3)/dtau;
1163 fPoly[i-1].D() = (divdf3/dtau)/dtau;
1164 }
1165}
1166
1167////////////////////////////////////////////////////////////////////////////////
1168/// Stream an object of class TSpline3.
1169
1171{
1172 if (R__b.IsReading()) {
1173 UInt_t R__s, R__c;
1174 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
1175 if (R__v > 1) {
1176 R__b.ReadClassBuffer(TSpline3::Class(), this, R__v, R__s, R__c);
1177 return;
1178 }
1179 //====process old versions before automatic schema evolution
1181 if (fNp > 0) {
1182 fPoly = new TSplinePoly3[fNp];
1183 for(Int_t i=0; i<fNp; ++i) {
1184 fPoly[i].Streamer(R__b);
1185 }
1186 }
1187 // R__b >> fPoly;
1188 R__b >> fValBeg;
1189 R__b >> fValEnd;
1190 R__b >> fBegCond;
1191 R__b >> fEndCond;
1192 } else {
1193 R__b.WriteClassBuffer(TSpline3::Class(),this);
1194 }
1195}
1196
1197/** \class TSpline5
1198 \ingroup Hist
1199 Class to create quintic natural splines to interpolate knots
1200 Arbitrary conditions can be introduced for first and second
1201 derivatives using double knots (see BuildCoeff) for more on this.
1202 Double knots are automatically introduced at ending points
1203 */
1204
1205////////////////////////////////////////////////////////////////////////////////
1206/// Quintic natural spline creator given an array of
1207/// arbitrary knots in increasing abscissa order and
1208/// possibly end point conditions.
1209
1210TSpline5::TSpline5(const char *title,
1211 Double_t x[], Double_t y[], Int_t n,
1212 const char *opt, Double_t b1, Double_t e1,
1213 Double_t b2, Double_t e2) :
1214 TSpline(title,-1, x[0], x[n-1], n, kFALSE)
1215{
1216 Int_t beg, end;
1217 const char *cb1, *ce1, *cb2, *ce2;
1218 fName="Spline5";
1219
1220 // Check endpoint conditions
1222
1223 // Create the polynomial terms and fill
1224 // them with node information
1225 fPoly = new TSplinePoly5[fNp];
1226 for (Int_t i=0; i<n; ++i) {
1227 fPoly[i+beg].X() = x[i];
1228 fPoly[i+beg].Y() = y[i];
1229 }
1230
1231 // Set the double knots at boundaries
1232 SetBoundaries(b1,e1,b2,e2,cb1,ce1,cb2,ce2);
1233
1234 // Build the spline coefficients
1235 BuildCoeff();
1236}
1237
1238////////////////////////////////////////////////////////////////////////////////
1239/// Quintic natural spline creator given an array of
1240/// arbitrary function values on equidistant n abscissa
1241/// values from xmin to xmax and possibly end point conditions.
1242
1243TSpline5::TSpline5(const char *title,
1245 Double_t y[], Int_t n,
1246 const char *opt, Double_t b1, Double_t e1,
1247 Double_t b2, Double_t e2) :
1248 TSpline(title,(xmax-xmin)/(n-1), xmin, xmax, n, kTRUE)
1249{
1250 Int_t beg, end;
1251 const char *cb1, *ce1, *cb2, *ce2;
1252 fName="Spline5";
1253
1254 // Check endpoint conditions
1256
1257 // Create the polynomial terms and fill
1258 // them with node information
1259 fPoly = new TSplinePoly5[fNp];
1260 for (Int_t i=0; i<n; ++i) {
1261 fPoly[i+beg].X() = fXmin+i*fDelta;
1262 fPoly[i+beg].Y() = y[i];
1263 }
1264
1265 // Set the double knots at boundaries
1266 SetBoundaries(b1,e1,b2,e2,cb1,ce1,cb2,ce2);
1267
1268 // Build the spline coefficients
1269 BuildCoeff();
1270}
1271
1272////////////////////////////////////////////////////////////////////////////////
1273/// Quintic natural spline creator given an array of
1274/// arbitrary abscissas in increasing order and a function
1275/// to interpolate and possibly end point conditions.
1276
1277TSpline5::TSpline5(const char *title,
1278 Double_t x[], const TF1 *func, Int_t n,
1279 const char *opt, Double_t b1, Double_t e1,
1280 Double_t b2, Double_t e2) :
1281 TSpline(title,-1, x[0], x[n-1], n, kFALSE)
1282{
1283 Int_t beg, end;
1284 const char *cb1, *ce1, *cb2, *ce2;
1285 fName="Spline5";
1286
1287 // Check endpoint conditions
1289
1290 // Create the polynomial terms and fill
1291 // them with node information
1292 fPoly = new TSplinePoly5[fNp];
1293 for (Int_t i=0; i<n; i++) {
1294 fPoly[i+beg].X() = x[i];
1295 fPoly[i+beg].Y() = ((TF1*)func)->Eval(x[i]);
1296 }
1297
1298 // Set the double knots at boundaries
1299 SetBoundaries(b1,e1,b2,e2,cb1,ce1,cb2,ce2);
1300
1301 // Build the spline coefficients
1302 BuildCoeff();
1303}
1304
1305////////////////////////////////////////////////////////////////////////////////
1306/// Quintic natural spline creator given a function to be
1307/// evaluated on n equidistant abscissa points between xmin
1308/// and xmax and possibly end point conditions.
1309
1310TSpline5::TSpline5(const char *title,
1312 const TF1 *func, Int_t n,
1313 const char *opt, Double_t b1, Double_t e1,
1314 Double_t b2, Double_t e2) :
1315 TSpline(title,(xmax-xmin)/(n-1), xmin, xmax, n, kTRUE)
1316{
1317 Int_t beg, end;
1318 const char *cb1, *ce1, *cb2, *ce2;
1319 fName="Spline5";
1320
1321 // Check endpoint conditions
1323
1324 // Create the polynomial terms and fill
1325 // them with node information
1326 fPoly = new TSplinePoly5[fNp];
1327 for (Int_t i=0; i<n; ++i) {
1329 fPoly[i+beg].X() = x;
1330 if (func) fPoly[i+beg].Y() = ((TF1*)func)->Eval(x);
1331 }
1332 if (!func) {fDelta = -1; fKstep = kFALSE;}
1333
1334 // Set the double knots at boundaries
1335 SetBoundaries(b1,e1,b2,e2,cb1,ce1,cb2,ce2);
1336
1337 // Build the spline coefficients
1338 if (func) BuildCoeff();
1339}
1340
1341////////////////////////////////////////////////////////////////////////////////
1342/// Quintic natural spline creator given a TGraph with
1343/// abscissa in increasing order and possibly end
1344/// point conditions.
1345
1346TSpline5::TSpline5(const char *title,
1347 const TGraph *g,
1348 const char *opt, Double_t b1, Double_t e1,
1349 Double_t b2, Double_t e2) :
1350 TSpline(title,-1,0,0,g->GetN(),kFALSE)
1351{
1352 Int_t beg, end;
1353 const char *cb1, *ce1, *cb2, *ce2;
1354 fName="Spline5";
1355
1356 // Check endpoint conditions
1358
1359 // Create the polynomial terms and fill
1360 // them with node information
1361 fPoly = new TSplinePoly5[fNp];
1362 for (Int_t i=0; i<fNp-beg; ++i) {
1363 Double_t xx, yy;
1364 g->GetPoint(i,xx,yy);
1365 fPoly[i+beg].X()=xx;
1366 fPoly[i+beg].Y()=yy;
1367 }
1368
1369 // Set the double knots at boundaries
1370 SetBoundaries(b1,e1,b2,e2,cb1,ce1,cb2,ce2);
1371 fXmin = fPoly[0].X();
1372 fXmax = fPoly[fNp-1].X();
1373
1374 // Build the spline coefficients
1375 BuildCoeff();
1376}
1377
1378////////////////////////////////////////////////////////////////////////////////
1379/// Quintic natural spline creator given a TH1.
1380
1382 const char *opt, Double_t b1, Double_t e1,
1383 Double_t b2, Double_t e2) :
1384 TSpline(h->GetTitle(),-1,0,0,h->GetNbinsX(),kFALSE)
1385{
1386 Int_t beg, end;
1387 const char *cb1, *ce1, *cb2, *ce2;
1388 fName=h->GetName();
1389
1390 // Check endpoint conditions
1392
1393 // Create the polynomial terms and fill
1394 // them with node information
1395 fPoly = new TSplinePoly5[fNp];
1396 for (Int_t i=0; i<fNp-beg; ++i) {
1397 fPoly[i+beg].X()=h->GetXaxis()->GetBinCenter(i+1);
1398 fPoly[i+beg].Y()=h->GetBinContent(i+1);
1399 }
1400
1401 // Set the double knots at boundaries
1402 SetBoundaries(b1,e1,b2,e2,cb1,ce1,cb2,ce2);
1403 fXmin = fPoly[0].X();
1404 fXmax = fPoly[fNp-1].X();
1405
1406 // Build the spline coefficients
1407 BuildCoeff();
1408}
1409
1410////////////////////////////////////////////////////////////////////////////////
1411/// Copy constructor.
1412
1414{
1415 if (fNp > 0)
1416 fPoly = new TSplinePoly5[fNp];
1417 for (Int_t i = 0; i < fNp; ++i)
1418 fPoly[i] = sp5.fPoly[i];
1419}
1420
1421////////////////////////////////////////////////////////////////////////////////
1422/// Assignment operator.
1423
1425{
1426 if (this != &sp5) {
1428 if (fPoly) {
1429 delete[] fPoly;
1430 fPoly = nullptr;
1431 }
1432 if (fNp > 0)
1433 fPoly = new TSplinePoly5[fNp];
1434 for (Int_t i = 0; i < fNp; ++i)
1435 fPoly[i] = sp5.fPoly[i];
1436 }
1437 return *this;
1438}
1439
1440////////////////////////////////////////////////////////////////////////////////
1441/// Check the boundary conditions and the
1442/// amount of extra double knots needed.
1443
1444void TSpline5::BoundaryConditions(const char *opt,Int_t &beg,Int_t &end,
1445 const char *&cb1,const char *&ce1,
1446 const char *&cb2,const char *&ce2)
1447{
1448 cb1=ce1=cb2=ce2=nullptr;
1449 beg=end=0;
1450 if(opt) {
1451 cb1 = strstr(opt,"b1");
1452 ce1 = strstr(opt,"e1");
1453 cb2 = strstr(opt,"b2");
1454 ce2 = strstr(opt,"e2");
1455 if(cb2) {
1456 fNp=fNp+2;
1457 beg=2;
1458 } else if(cb1) {
1459 fNp=fNp+1;
1460 beg=1;
1461 }
1462 if(ce2) {
1463 fNp=fNp+2;
1464 end=2;
1465 } else if(ce1) {
1466 fNp=fNp+1;
1467 end=1;
1468 }
1469 }
1470}
1471
1472////////////////////////////////////////////////////////////////////////////////
1473/// Set the boundary conditions at double/triple knots.
1474
1476 const char *cb1, const char *ce1, const char *cb2,
1477 const char *ce2)
1478{
1479 if(cb2) {
1480
1481 // Second derivative at the beginning
1482 fPoly[0].X() = fPoly[1].X() = fPoly[2].X();
1483 fPoly[0].Y() = fPoly[2].Y();
1484 fPoly[2].Y()=b2;
1485
1486 // If first derivative not given, we take the finite
1487 // difference from first and second point... not recommended
1488 if(cb1)
1489 fPoly[1].Y()=b1;
1490 else
1491 fPoly[1].Y()=(fPoly[3].Y()-fPoly[0].Y())/(fPoly[3].X()-fPoly[2].X());
1492 } else if(cb1) {
1493
1494 // First derivative at the end
1495 fPoly[0].X() = fPoly[1].X();
1496 fPoly[0].Y() = fPoly[1].Y();
1497 fPoly[1].Y()=b1;
1498 }
1499 if(ce2) {
1500
1501 // Second derivative at the end
1502 fPoly[fNp-1].X() = fPoly[fNp-2].X() = fPoly[fNp-3].X();
1503 fPoly[fNp-1].Y()=e2;
1504
1505 // If first derivative not given, we take the finite
1506 // difference from first and second point... not recommended
1507 if(ce1)
1508 fPoly[fNp-2].Y()=e1;
1509 else
1510 fPoly[fNp-2].Y()=
1511 (fPoly[fNp-3].Y()-fPoly[fNp-4].Y())
1512 /(fPoly[fNp-3].X()-fPoly[fNp-4].X());
1513 } else if(ce1) {
1514
1515 // First derivative at the end
1516 fPoly[fNp-1].X() = fPoly[fNp-2].X();
1517 fPoly[fNp-1].Y()=e1;
1518 }
1519}
1520
1521////////////////////////////////////////////////////////////////////////////////
1522/// Find X.
1523
1525{
1526 Int_t klow=0;
1527
1528 // If out of boundaries, extrapolate
1529 // It may be badly wrong
1530 if(x<=fXmin) klow=0;
1531 else if(x>=fXmax) klow=fNp-1;
1532 else {
1533 if(fKstep) {
1534
1535 // Equidistant knots, use histogramming
1537 } else {
1538 Int_t khig=fNp-1, khalf;
1539
1540 // Non equidistant knots, binary search
1541 while(khig-klow>1)
1542 if(x>fPoly[khalf=(klow+khig)/2].X())
1543 klow=khalf;
1544 else
1545 khig=khalf;
1546 }
1547
1548 // This could be removed, sanity check
1549 if(!(fPoly[klow].X()<=x && x<=fPoly[klow+1].X()))
1550 Error("Eval",
1551 "Binary search failed x(%d) = %f < x(%d) = %f\n",
1552 klow,fPoly[klow].X(),klow+1,fPoly[klow+1].X());
1553 }
1554 return klow;
1555}
1556
1557////////////////////////////////////////////////////////////////////////////////
1558/// Eval this spline at x.
1559
1561{
1562 Int_t klow=FindX(x);
1563 return fPoly[klow].Eval(x);
1564}
1565
1566////////////////////////////////////////////////////////////////////////////////
1567/// Derivative.
1568
1570{
1571 Int_t klow=FindX(x);
1572 return fPoly[klow].Derivative(x);
1573}
1574
1575////////////////////////////////////////////////////////////////////////////////
1576/// Write this spline as a C++ function that can be executed without ROOT
1577/// the name of the function is the name of the file up to the "." if any.
1578
1579void TSpline5::SaveAs(const char *filename, Option_t * /*option*/) const
1580{
1581 //open the file
1582 std::ofstream *f = new std::ofstream(filename,std::ios::out);
1583 if (f == nullptr || gSystem->AccessPathName(filename,kWritePermission)) {
1584 Error("SaveAs","Cannot open file:%s\n",filename);
1585 return;
1586 }
1587
1588 //write the function name and the spline constants
1589 char buffer[512];
1591 snprintf(buffer,512,"double %s",filename);
1592 char *dot = strstr(buffer,".");
1593 if (dot) *dot = 0;
1594 strlcat(buffer,"(double x) {\n",512);
1595 nch = strlen(buffer); f->write(buffer,nch);
1596 snprintf(buffer,512," const int fNp = %d, fKstep = %d;\n",fNp,fKstep);
1597 nch = strlen(buffer); f->write(buffer,nch);
1598 snprintf(buffer,512," const double fDelta = %.17g, fXmin = %.17g, fXmax = %.17g;\n",fDelta,fXmin,fXmax);
1599 nch = strlen(buffer); f->write(buffer,nch);
1600
1601 //write the spline coefficients
1602 //array fX
1603 snprintf(buffer,512," const double fX[%d] = {",fNp);
1604 nch = strlen(buffer); f->write(buffer,nch);
1605 buffer[0] = 0;
1606 Int_t i;
1607 char numb[30];
1608 for (i=0;i<fNp;i++) {
1609 snprintf(numb,30," %.17g,",fPoly[i].X());
1610 nch = strlen(numb);
1611 if (i == fNp-1) numb[nch-1]=0;
1612 strlcat(buffer,numb,512);
1613 if (i%5 == 4 || i == fNp-1) {
1614 nch = strlen(buffer); f->write(buffer,nch);
1615 if (i != fNp-1) snprintf(buffer,512,"\n ");
1616 }
1617 }
1618 snprintf(buffer,512," };\n");
1619 nch = strlen(buffer); f->write(buffer,nch);
1620 //array fY
1621 snprintf(buffer,512," const double fY[%d] = {",fNp);
1622 nch = strlen(buffer); f->write(buffer,nch);
1623 buffer[0] = 0;
1624 for (i=0;i<fNp;i++) {
1625 snprintf(numb,30," %.17g,",fPoly[i].Y());
1626 nch = strlen(numb);
1627 if (i == fNp-1) numb[nch-1]=0;
1628 strlcat(buffer,numb,512);
1629 if (i%5 == 4 || i == fNp-1) {
1630 nch = strlen(buffer); f->write(buffer,nch);
1631 if (i != fNp-1) snprintf(buffer,512,"\n ");
1632 }
1633 }
1634 snprintf(buffer,512," };\n");
1635 nch = strlen(buffer); f->write(buffer,nch);
1636 //array fB
1637 snprintf(buffer,512," const double fB[%d] = {",fNp);
1638 nch = strlen(buffer); f->write(buffer,nch);
1639 buffer[0] = 0;
1640 for (i=0;i<fNp;i++) {
1641 snprintf(numb,30," %.17g,",fPoly[i].B());
1642 nch = strlen(numb);
1643 if (i == fNp-1) numb[nch-1]=0;
1644 strlcat(buffer,numb,512);
1645 if (i%5 == 4 || i == fNp-1) {
1646 nch = strlen(buffer); f->write(buffer,nch);
1647 if (i != fNp-1) snprintf(buffer,512,"\n ");
1648 }
1649 }
1650 snprintf(buffer,512," };\n");
1651 nch = strlen(buffer); f->write(buffer,nch);
1652 //array fC
1653 snprintf(buffer,512," const double fC[%d] = {",fNp);
1654 nch = strlen(buffer); f->write(buffer,nch);
1655 buffer[0] = 0;
1656 for (i=0;i<fNp;i++) {
1657 snprintf(numb,30," %.17g,",fPoly[i].C());
1658 nch = strlen(numb);
1659 if (i == fNp-1) numb[nch-1]=0;
1660 strlcat(buffer,numb,512);
1661 if (i%5 == 4 || i == fNp-1) {
1662 nch = strlen(buffer); f->write(buffer,nch);
1663 if (i != fNp-1) snprintf(buffer,512,"\n ");
1664 }
1665 }
1666 snprintf(buffer,512," };\n");
1667 nch = strlen(buffer); f->write(buffer,nch);
1668 //array fD
1669 snprintf(buffer,512," const double fD[%d] = {",fNp);
1670 nch = strlen(buffer); f->write(buffer,nch);
1671 buffer[0] = 0;
1672 for (i=0;i<fNp;i++) {
1673 snprintf(numb,30," %.17g,",fPoly[i].D());
1674 nch = strlen(numb);
1675 if (i == fNp-1) numb[nch-1]=0;
1676 strlcat(buffer,numb,512);
1677 if (i%5 == 4 || i == fNp-1) {
1678 nch = strlen(buffer); f->write(buffer,nch);
1679 if (i != fNp-1) snprintf(buffer,512,"\n ");
1680 }
1681 }
1682 snprintf(buffer,512," };\n");
1683 nch = strlen(buffer); f->write(buffer,nch);
1684 //array fE
1685 snprintf(buffer,512," const double fE[%d] = {",fNp);
1686 nch = strlen(buffer); f->write(buffer,nch);
1687 buffer[0] = 0;
1688 for (i=0;i<fNp;i++) {
1689 snprintf(numb,30," %.17g,",fPoly[i].E());
1690 nch = strlen(numb);
1691 if (i == fNp-1) numb[nch-1]=0;
1692 strlcat(buffer,numb,512);
1693 if (i%5 == 4 || i == fNp-1) {
1694 nch = strlen(buffer); f->write(buffer,nch);
1695 if (i != fNp-1) snprintf(buffer,512,"\n ");
1696 }
1697 }
1698 snprintf(buffer,512," };\n");
1699 nch = strlen(buffer); f->write(buffer,nch);
1700 //array fF
1701 snprintf(buffer,512," const double fF[%d] = {",fNp);
1702 nch = strlen(buffer); f->write(buffer,nch);
1703 buffer[0] = 0;
1704 for (i=0;i<fNp;i++) {
1705 snprintf(numb,30," %.17g,",fPoly[i].F());
1706 nch = strlen(numb);
1707 if (i == fNp-1) numb[nch-1]=0;
1708 strlcat(buffer,numb,512);
1709 if (i%5 == 4 || i == fNp-1) {
1710 nch = strlen(buffer); f->write(buffer,nch);
1711 if (i != fNp-1) snprintf(buffer,512,"\n ");
1712 }
1713 }
1714 snprintf(buffer,512," };\n");
1715 nch = strlen(buffer); f->write(buffer,nch);
1716
1717 //generate code for the spline evaluation
1718 snprintf(buffer,512," int klow=0;\n");
1719 nch = strlen(buffer); f->write(buffer,nch);
1720
1721 snprintf(buffer,512," // If out of boundaries, extrapolate. It may be badly wrong\n");
1722 snprintf(buffer,512," if(x<=fXmin) klow=0;\n");
1723 nch = strlen(buffer); f->write(buffer,nch);
1724 snprintf(buffer,512," else if(x>=fXmax) klow=fNp-1;\n");
1725 nch = strlen(buffer); f->write(buffer,nch);
1726 snprintf(buffer,512," else {\n");
1727 nch = strlen(buffer); f->write(buffer,nch);
1728 snprintf(buffer,512," if(fKstep) {\n");
1729 nch = strlen(buffer); f->write(buffer,nch);
1730
1731 snprintf(buffer,512," // Equidistant knots, use histogramming\n");
1732 nch = strlen(buffer); f->write(buffer,nch);
1733 snprintf(buffer,512," klow = int((x-fXmin)/fDelta);\n");
1734 nch = strlen(buffer); f->write(buffer,nch);
1735 snprintf(buffer,512," if (klow > fNp-1) klow = fNp-1;\n");
1736 nch = strlen(buffer); f->write(buffer,nch);
1737 snprintf(buffer,512," } else {\n");
1738 nch = strlen(buffer); f->write(buffer,nch);
1739 snprintf(buffer,512," int khig=fNp-1, khalf;\n");
1740 nch = strlen(buffer); f->write(buffer,nch);
1741
1742 snprintf(buffer,512," // Non equidistant knots, binary search\n");
1743 nch = strlen(buffer); f->write(buffer,nch);
1744 snprintf(buffer,512," while(khig-klow>1)\n");
1745 nch = strlen(buffer); f->write(buffer,nch);
1746 snprintf(buffer,512," if(x>fX[khalf=(klow+khig)/2]) klow=khalf;\n");
1747 nch = strlen(buffer); f->write(buffer,nch);
1748 snprintf(buffer,512," else khig=khalf;\n");
1749 nch = strlen(buffer); f->write(buffer,nch);
1750 snprintf(buffer,512," }\n");
1751 nch = strlen(buffer); f->write(buffer,nch);
1752 snprintf(buffer,512," }\n");
1753 nch = strlen(buffer); f->write(buffer,nch);
1754 snprintf(buffer,512," // Evaluate now\n");
1755 nch = strlen(buffer); f->write(buffer,nch);
1756 snprintf(buffer,512," double dx=x-fX[klow];\n");
1757 nch = strlen(buffer); f->write(buffer,nch);
1758 snprintf(buffer,512," return (fY[klow]+dx*(fB[klow]+dx*(fC[klow]+dx*(fD[klow]+dx*(fE[klow]+dx*fF[klow])))));\n");
1759 nch = strlen(buffer); f->write(buffer,nch);
1760
1761 //close file
1762 f->write("}\n",2);
1763
1764 if (f) { f->close(); delete f;}
1765}
1766
1767////////////////////////////////////////////////////////////////////////////////
1768/// Save primitive as a C++ statement(s) on output stream out.
1769
1770void TSpline5::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
1771{
1772 Double_t b1 = fPoly[1].Y();
1773 Double_t e1 = fPoly[fNp - 1].Y();
1774 Double_t b2 = fPoly[2].Y();
1775 Double_t e2 = fPoly[fNp - 1].Y();
1776
1777 SavePrimitiveConstructor(out, Class(), "spline5",
1778 TString::Format("\"%s\", %g, %g, (TF1 *)nullptr, %d, \"\", %g, %g, %g, %g",
1779 TString(GetTitle()).ReplaceSpecialCppChars().Data(), fXmin, fXmax, fNp, b1,
1780 e1, b2, e2));
1781
1782 out << " spline5->SetName(\"" << TString(GetName()).ReplaceSpecialCppChars() << "\");\n";
1783
1784 SaveFillAttributes(out, "spline5", 0, 1001);
1785 SaveLineAttributes(out, "spline5", 1, 1, 1);
1786 SaveMarkerAttributes(out, "spline5", 1, 1, 1);
1787 if (fNpx != 100)
1788 out << " spline5->SetNpx(" << fNpx << ");\n";
1789
1790 for (Int_t i = 0; i < fNp; i++) {
1791 out << " spline5->SetPoint(" << i << "," << fPoly[i].X() << "," << fPoly[i].Y() << ");\n";
1792 out << " spline5->SetPointCoeff(" << i << "," << fPoly[i].B() << "," << fPoly[i].C() << "," << fPoly[i].D()
1793 << "," << fPoly[i].E() << "," << fPoly[i].F() << ");\n";
1794 }
1795
1796 out << " spline5->Draw(\"" << TString(option).ReplaceSpecialCppChars() << "\");\n";
1797}
1798
1799////////////////////////////////////////////////////////////////////////////////
1800/// Set point number i.
1801
1803{
1804
1805 if (i < 0 || i >= fNp) return;
1806 fPoly[i].X()= x;
1807 fPoly[i].Y()= y;
1808}
1809
1810////////////////////////////////////////////////////////////////////////////////
1811/// Set point coefficient number i.
1812
1815{
1816 if (i < 0 || i >= fNp) return;
1817 fPoly[i].B()= b;
1818 fPoly[i].C()= c;
1819 fPoly[i].D()= d;
1820 fPoly[i].E()= e;
1821 fPoly[i].F()= f;
1822}
1823
1824////////////////////////////////////////////////////////////////////////////////
1825/// Algorithm 600, collected algorithms from acm.
1826///
1827/// algorithm appeared in acm-trans. math. software, vol.9, no. 2,
1828/// jun., 1983, p. 258-259.
1829///
1830/// TSpline5 computes the coefficients of a quintic natural quintic spli
1831/// s(x) with knots x(i) interpolating there to given function values:
1832/// ~~~ {.cpp}
1833/// s(x(i)) = y(i) for i = 1,2, ..., n.
1834/// ~~~
1835/// in each interval (x(i),x(i+1)) the spline function s(xx) is a
1836/// polynomial of fifth degree:
1837/// ~~~ {.cpp}
1838/// s(xx) = ((((f(i)*p+e(i))*p+d(i))*p+c(i))*p+b(i))*p+y(i) (*)
1839/// = ((((-f(i)*q+e(i+1))*q-d(i+1))*q+c(i+1))*q-b(i+1))*q+y(i+1)
1840/// ~~~
1841/// where p = xx - x(i) and q = x(i+1) - xx.
1842/// (note the first subscript in the second expression.)
1843/// the different polynomials are pieced together so that s(x) and
1844/// its derivatives up to s"" are continuous.
1845///
1846/// ### input:
1847///
1848/// n number of data points, (at least three, i.e. n > 2)
1849/// x(1:n) the strictly increasing or decreasing sequence of
1850/// knots. the spacing must be such that the fifth power
1851/// of x(i+1) - x(i) can be formed without overflow or
1852/// underflow of exponents.
1853/// y(1:n) the prescribed function values at the knots.
1854///
1855/// ### output:
1856///
1857/// b,c,d,e,f the computed spline coefficients as in (*).
1858/// (1:n) specifically
1859/// b(i) = s'(x(i)), c(i) = s"(x(i))/2, d(i) = s"'(x(i))/6,
1860/// e(i) = s""(x(i))/24, f(i) = s""'(x(i))/120.
1861/// f(n) is neither used nor altered. the five arrays
1862/// b,c,d,e,f must always be distinct.
1863///
1864/// ### option:
1865///
1866/// it is possible to specify values for the first and second
1867/// derivatives of the spline function at arbitrarily many knots.
1868/// this is done by relaxing the requirement that the sequence of
1869/// knots be strictly increasing or decreasing. specifically:
1870///
1871/// ~~~ {.cpp}
1872/// if x(j) = x(j+1) then s(x(j)) = y(j) and s'(x(j)) = y(j+1),
1873/// if x(j) = x(j+1) = x(j+2) then in addition s"(x(j)) = y(j+2).
1874/// ~~~
1875///
1876/// note that s""(x) is discontinuous at a double knot and, in
1877/// addition, s"'(x) is discontinuous at a triple knot. the
1878/// subroutine assigns y(i) to y(i+1) in these cases and also to
1879/// y(i+2) at a triple knot. the representation (*) remains
1880/// valid in each open interval (x(i),x(i+1)). at a double knot,
1881/// x(j) = x(j+1), the output coefficients have the following values:
1882/// ~~~ {.cpp}
1883/// y(j) = s(x(j)) = y(j+1)
1884/// b(j) = s'(x(j)) = b(j+1)
1885/// c(j) = s"(x(j))/2 = c(j+1)
1886/// d(j) = s"'(x(j))/6 = d(j+1)
1887/// e(j) = s""(x(j)-0)/24 e(j+1) = s""(x(j)+0)/24
1888/// f(j) = s""'(x(j)-0)/120 f(j+1) = s""'(x(j)+0)/120
1889/// ~~~
1890/// at a triple knot, x(j) = x(j+1) = x(j+2), the output
1891/// coefficients have the following values:
1892/// ~~~ {.cpp}
1893/// y(j) = s(x(j)) = y(j+1) = y(j+2)
1894/// b(j) = s'(x(j)) = b(j+1) = b(j+2)
1895/// c(j) = s"(x(j))/2 = c(j+1) = c(j+2)
1896/// d(j) = s"'((x(j)-0)/6 d(j+1) = 0 d(j+2) = s"'(x(j)+0)/6
1897/// e(j) = s""(x(j)-0)/24 e(j+1) = 0 e(j+2) = s""(x(j)+0)/24
1898/// f(j) = s""'(x(j)-0)/120 f(j+1) = 0 f(j+2) = s""'(x(j)+0)/120
1899/// ~~~
1900
1902{
1903 Int_t i, m;
1904 Double_t pqqr, p, q, r, s, t, u, v,
1905 b1, p2, p3, q2, q3, r2, pq, pr, qr;
1906
1907 if (fNp <= 2) {
1908 return;
1909 }
1910
1911 // coefficients of a positive definite, pentadiagonal matrix,
1912 // stored in D, E, F from 1 to n-3.
1913 m = fNp-2;
1914 q = fPoly[1].X()-fPoly[0].X();
1915 r = fPoly[2].X()-fPoly[1].X();
1916 q2 = q*q;
1917 r2 = r*r;
1918 qr = q+r;
1919 fPoly[0].D() = fPoly[0].E() = 0;
1920 if (q) fPoly[1].D() = q*6.*q2/(qr*qr);
1921 else fPoly[1].D() = 0;
1922
1923 if (m > 1) {
1924 for (i = 1; i < m; ++i) {
1925 p = q;
1926 q = r;
1927 r = fPoly[i+2].X()-fPoly[i+1].X();
1928 p2 = q2;
1929 q2 = r2;
1930 r2 = r*r;
1931 pq = qr;
1932 qr = q+r;
1933 if (q) {
1934 q3 = q2*q;
1935 pr = p*r;
1936 pqqr = pq*qr;
1937 fPoly[i+1].D() = q3*6./(qr*qr);
1938 fPoly[i].D() += (q+q)*(pr*15.*pr+(p+r)*q
1939 *(pr* 20.+q2*7.)+q2*
1940 ((p2+r2)*8.+pr*21.+q2+q2))/(pqqr*pqqr);
1941 fPoly[i-1].D() += q3*6./(pq*pq);
1942 fPoly[i].E() = q2*(p*qr+pq*3.*(qr+r+r))/(pqqr*qr);
1943 fPoly[i-1].E() += q2*(r*pq+qr*3.*(pq+p+p))/(pqqr*pq);
1944 fPoly[i-1].F() = q3/pqqr;
1945 } else
1946 fPoly[i+1].D() = fPoly[i].E() = fPoly[i-1].F() = 0;
1947 }
1948 }
1949 if (r) fPoly[m-1].D() += r*6.*r2/(qr*qr);
1950
1951 // First and second order divided differences of the given function
1952 // values, stored in b from 2 to n and in c from 3 to n
1953 // respectively. care is taken of double and triple knots.
1954 for (i = 1; i < fNp; ++i) {
1955 if (fPoly[i].X() != fPoly[i-1].X()) {
1956 fPoly[i].B() =
1957 (fPoly[i].Y()-fPoly[i-1].Y())/(fPoly[i].X()-fPoly[i-1].X());
1958 } else {
1959 fPoly[i].B() = fPoly[i].Y();
1960 fPoly[i].Y() = fPoly[i-1].Y();
1961 }
1962 }
1963 for (i = 2; i < fNp; ++i) {
1964 if (fPoly[i].X() != fPoly[i-2].X()) {
1965 fPoly[i].C() =
1966 (fPoly[i].B()-fPoly[i-1].B())/(fPoly[i].X()-fPoly[i-2].X());
1967 } else {
1968 fPoly[i].C() = fPoly[i].B()*.5;
1969 fPoly[i].B() = fPoly[i-1].B();
1970 }
1971 }
1972
1973 // Solve the linear system with c(i+2) - c(i+1) as right-hand side. */
1974 if (m > 1) {
1975 p=fPoly[0].C()=fPoly[m-1].E()=fPoly[0].F()
1976 =fPoly[m-2].F()=fPoly[m-1].F()=0;
1977 fPoly[1].C() = fPoly[3].C()-fPoly[2].C();
1978 fPoly[1].D() = 1./fPoly[1].D();
1979
1980 if (m > 2) {
1981 for (i = 2; i < m; ++i) {
1982 q = fPoly[i-1].D()*fPoly[i-1].E();
1983 fPoly[i].D() = 1./(fPoly[i].D()-p*fPoly[i-2].F()-q*fPoly[i-1].E());
1984 fPoly[i].E() -= q*fPoly[i-1].F();
1985 fPoly[i].C() = fPoly[i+2].C()-fPoly[i+1].C()-p*fPoly[i-2].C()
1986 -q*fPoly[i-1].C();
1987 p = fPoly[i-1].D()*fPoly[i-1].F();
1988 }
1989 }
1990 }
1991
1992 fPoly[fNp-2].C() = fPoly[fNp-1].C() = 0;
1993 if (fNp > 3)
1994 for (i=fNp-3; i > 0; --i)
1995 fPoly[i].C() = (fPoly[i].C()-fPoly[i].E()*fPoly[i+1].C()
1996 -fPoly[i].F()*fPoly[i+2].C())*fPoly[i].D();
1997
1998 // Integrate the third derivative of s(x)
1999 m = fNp-1;
2000 q = fPoly[1].X()-fPoly[0].X();
2001 r = fPoly[2].X()-fPoly[1].X();
2002 b1 = fPoly[1].B();
2003 q3 = q*q*q;
2004 qr = q+r;
2005 if (qr) {
2006 v = fPoly[1].C()/qr;
2007 t = v;
2008 } else
2009 v = t = 0;
2010 if (q) fPoly[0].F() = v/q;
2011 else fPoly[0].F() = 0;
2012 for (i = 1; i < m; ++i) {
2013 p = q;
2014 q = r;
2015 if (i != m-1) r = fPoly[i+2].X()-fPoly[i+1].X();
2016 else r = 0;
2017 p3 = q3;
2018 q3 = q*q*q;
2019 pq = qr;
2020 qr = q+r;
2021 s = t;
2022 if (qr) t = (fPoly[i+1].C()-fPoly[i].C())/qr;
2023 else t = 0;
2024 u = v;
2025 v = t-s;
2026 if (pq) {
2027 fPoly[i].F() = fPoly[i-1].F();
2028 if (q) fPoly[i].F() = v/q;
2029 fPoly[i].E() = s*5.;
2030 fPoly[i].D() = (fPoly[i].C()-q*s)*10;
2031 fPoly[i].C() =
2032 fPoly[i].D()*(p-q)+(fPoly[i+1].B()-fPoly[i].B()+(u-fPoly[i].E())*
2033 p3-(v+fPoly[i].E())*q3)/pq;
2034 fPoly[i].B() = (p*(fPoly[i+1].B()-v*q3)+q*(fPoly[i].B()-u*p3))/pq-p
2035 *q*(fPoly[i].D()+fPoly[i].E()*(q-p));
2036 } else {
2037 fPoly[i].C() = fPoly[i-1].C();
2038 fPoly[i].D() = fPoly[i].E() = fPoly[i].F() = 0;
2039 }
2040 }
2041
2042 // End points x(1) and x(n)
2043 p = fPoly[1].X()-fPoly[0].X();
2044 s = fPoly[0].F()*p*p*p;
2045 fPoly[0].E() = fPoly[0].D() = 0;
2046 fPoly[0].C() = fPoly[1].C()-s*10;
2047 fPoly[0].B() = b1-(fPoly[0].C()+s)*p;
2048
2049 q = fPoly[fNp-1].X()-fPoly[fNp-2].X();
2050 t = fPoly[fNp-2].F()*q*q*q;
2051 fPoly[fNp-1].E() = fPoly[fNp-1].D() = 0;
2052 fPoly[fNp-1].C() = fPoly[fNp-2].C()+t*10;
2053 fPoly[fNp-1].B() += (fPoly[fNp-1].C()-t)*q;
2054}
2055
2056////////////////////////////////////////////////////////////////////////////////
2057/// Test method for TSpline5
2058///
2059/// ~~~ {.cpp}
2060/// n number of data points.
2061/// m 2*m-1 is order of spline.
2062/// m = 3 always for quintic spline.
2063/// nn,nm1,mm,
2064/// mm1,i,k,
2065/// j,jj temporary integer variables.
2066/// z,p temporary double precision variables.
2067/// x[n] the sequence of knots.
2068/// y[n] the prescribed function values at the knots.
2069/// a[200][6] two dimensional array whose columns are
2070/// the computed spline coefficients
2071/// diff[5] maximum values of differences of values and
2072/// derivatives to right and left of knots.
2073/// com[5] maximum values of coefficients.
2074/// ~~~
2075///
2076/// test of TSpline5 with non equidistant knots and
2077/// equidistant knots follows.
2078
2080{
2081 Double_t hx;
2082 Double_t diff[5];
2083 Double_t a[1200], c[6];
2084 Int_t i, j, k, m, n;
2085 Double_t p, x[200], y[200], z;
2086 Int_t jj, mm, nn;
2087 Int_t mm1, nm1;
2088 Double_t com[5];
2089
2090 printf("1 TEST OF TSpline5 WITH NONEQUIDISTANT KNOTS\n");
2091 n = 5;
2092 x[0] = -3;
2093 x[1] = -1;
2094 x[2] = 0;
2095 x[3] = 3;
2096 x[4] = 4;
2097 y[0] = 7;
2098 y[1] = 11;
2099 y[2] = 26;
2100 y[3] = 56;
2101 y[4] = 29;
2102 m = 3;
2103 mm = m << 1;
2104 mm1 = mm-1;
2105 printf("\n-N = %3d M =%2d\n",n,m);
2106 TSpline5 *spline = new TSpline5("Test",x,y,n);
2107 for (i = 0; i < n; ++i)
2108 spline->GetCoeff(i,hx, a[i],a[i+200],a[i+400],
2109 a[i+600],a[i+800],a[i+1000]);
2110 delete spline;
2111 for (i = 0; i < mm1; ++i) diff[i] = com[i] = 0;
2112 for (k = 0; k < n; ++k) {
2113 for (i = 0; i < mm; ++i) c[i] = a[k+i*200];
2114 printf(" ---------------------------------------%3d --------------------------------------------\n",k+1);
2115 printf("%12.8f\n",x[k]);
2116 if (k == n-1) {
2117 printf("%16.8f\n",c[0]);
2118 } else {
2119 for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
2120 printf("\n");
2121 for (i = 0; i < mm1; ++i)
2122 if ((z=TMath::Abs(a[k+i*200])) > com[i]) com[i] = z;
2123 z = x[k+1]-x[k];
2124 for (i = 1; i < mm; ++i)
2125 for (jj = i; jj < mm; ++jj) {
2126 j = mm+i-jj;
2127 c[j-2] = c[j-1]*z+c[j-2];
2128 }
2129 for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
2130 printf("\n");
2131 for (i = 0; i < mm1; ++i)
2132 if (!(k >= n-2 && i != 0))
2133 if((z = TMath::Abs(c[i]-a[k+1+i*200]))
2134 > diff[i]) diff[i] = z;
2135 }
2136 }
2137 printf(" MAXIMUM ABSOLUTE VALUES OF DIFFERENCES \n");
2138 for (i = 0; i < mm1; ++i) printf("%18.9E",diff[i]);
2139 printf("\n");
2140 printf(" MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS \n");
2141 if (TMath::Abs(c[0]) > com[0])
2142 com[0] = TMath::Abs(c[0]);
2143 for (i = 0; i < mm1; ++i) printf("%16.8f",com[i]);
2144 printf("\n");
2145 m = 3;
2146 for (n = 10; n <= 100; n += 10) {
2147 mm = m << 1;
2148 mm1 = mm-1;
2149 nm1 = n-1;
2150 for (i = 0; i < nm1; i += 2) {
2151 x[i] = i+1;
2152 x[i+1] = i+2;
2153 y[i] = 1;
2154 y[i+1] = 0;
2155 }
2156 if (n % 2 != 0) {
2157 x[n-1] = n;
2158 y[n-1] = 1;
2159 }
2160 printf("\n-N = %3d M =%2d\n",n,m);
2161 spline = new TSpline5("Test",x,y,n);
2162 for (i = 0; i < n; ++i)
2163 spline->GetCoeff(i,hx,a[i],a[i+200],a[i+400],
2164 a[i+600],a[i+800],a[i+1000]);
2165 delete spline;
2166 for (i = 0; i < mm1; ++i)
2167 diff[i] = com[i] = 0;
2168 for (k = 0; k < n; ++k) {
2169 for (i = 0; i < mm; ++i)
2170 c[i] = a[k+i*200];
2171 if (n < 11) {
2172 printf(" ---------------------------------------%3d --------------------------------------------\n",k+1);
2173 printf("%12.8f\n",x[k]);
2174 if (k == n-1) printf("%16.8f\n",c[0]);
2175 }
2176 if (k == n-1) break;
2177 if (n <= 10) {
2178 for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
2179 printf("\n");
2180 }
2181 for (i = 0; i < mm1; ++i)
2182 if ((z=TMath::Abs(a[k+i*200])) > com[i])
2183 com[i] = z;
2184 z = x[k+1]-x[k];
2185 for (i = 1; i < mm; ++i)
2186 for (jj = i; jj < mm; ++jj) {
2187 j = mm+i-jj;
2188 c[j-2] = c[j-1]*z+c[j-2];
2189 }
2190 if (n <= 10) {
2191 for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
2192 printf("\n");
2193 }
2194 for (i = 0; i < mm1; ++i)
2195 if (!(k >= n-2 && i != 0))
2196 if ((z = TMath::Abs(c[i]-a[k+1+i*200]))
2197 > diff[i]) diff[i] = z;
2198 }
2199 printf(" MAXIMUM ABSOLUTE VALUES OF DIFFERENCES \n");
2200 for (i = 0; i < mm1; ++i) printf("%18.9E",diff[i]);
2201 printf("\n");
2202 printf(" MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS \n");
2203 if (TMath::Abs(c[0]) > com[0])
2204 com[0] = TMath::Abs(c[0]);
2205 for (i = 0; i < mm1; ++i) printf("%16.8E",com[i]);
2206 printf("\n");
2207 }
2208
2209 // Test of TSpline5 with non equidistant double knots follows
2210 printf("1 TEST OF TSpline5 WITH NONEQUIDISTANT DOUBLE KNOTS\n");
2211 n = 5;
2212 x[0] = -3;
2213 x[1] = -3;
2214 x[2] = -1;
2215 x[3] = -1;
2216 x[4] = 0;
2217 x[5] = 0;
2218 x[6] = 3;
2219 x[7] = 3;
2220 x[8] = 4;
2221 x[9] = 4;
2222 y[0] = 7;
2223 y[1] = 2;
2224 y[2] = 11;
2225 y[3] = 15;
2226 y[4] = 26;
2227 y[5] = 10;
2228 y[6] = 56;
2229 y[7] = -27;
2230 y[8] = 29;
2231 y[9] = -30;
2232 m = 3;
2233 nn = n << 1;
2234 mm = m << 1;
2235 mm1 = mm-1;
2236 printf("-N = %3d M =%2d\n",n,m);
2237 spline = new TSpline5("Test",x,y,nn);
2238 for (i = 0; i < nn; ++i)
2239 spline->GetCoeff(i,hx,a[i],a[i+200],a[i+400],
2240 a[i+600],a[i+800],a[i+1000]);
2241 delete spline;
2242 for (i = 0; i < mm1; ++i)
2243 diff[i] = com[i] = 0;
2244 for (k = 0; k < nn; ++k) {
2245 for (i = 0; i < mm; ++i)
2246 c[i] = a[k+i*200];
2247 printf(" ---------------------------------------%3d --------------------------------------------\n",k+1);
2248 printf("%12.8f\n",x[k]);
2249 if (k == nn-1) {
2250 printf("%16.8f\n",c[0]);
2251 break;
2252 }
2253 for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
2254 printf("\n");
2255 for (i = 0; i < mm1; ++i)
2256 if ((z=TMath::Abs(a[k+i*200])) > com[i]) com[i] = z;
2257 z = x[k+1]-x[k];
2258 for (i = 1; i < mm; ++i)
2259 for (jj = i; jj < mm; ++jj) {
2260 j = mm+i-jj;
2261 c[j-2] = c[j-1]*z+c[j-2];
2262 }
2263 for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
2264 printf("\n");
2265 for (i = 0; i < mm1; ++i)
2266 if (!(k >= nn-2 && i != 0))
2267 if ((z = TMath::Abs(c[i]-a[k+1+i*200]))
2268 > diff[i]) diff[i] = z;
2269 }
2270 printf(" MAXIMUM ABSOLUTE VALUES OF DIFFERENCES \n");
2271 for (i = 1; i <= mm1; ++i) {
2272 printf("%18.9E",diff[i-1]);
2273 }
2274 printf("\n");
2275 if (TMath::Abs(c[0]) > com[0])
2276 com[0] = TMath::Abs(c[0]);
2277 printf(" MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS \n");
2278 for (i = 0; i < mm1; ++i) printf("%16.8f",com[i]);
2279 printf("\n");
2280 m = 3;
2281 for (n = 10; n <= 100; n += 10) {
2282 nn = n << 1;
2283 mm = m << 1;
2284 mm1 = mm-1;
2285 p = 0;
2286 for (i = 0; i < n; ++i) {
2287 p += TMath::Abs(TMath::Sin(i+1));
2288 x[(i << 1)] = p;
2289 x[(i << 1)+1] = p;
2290 y[(i << 1)] = TMath::Cos(i+1)-.5;
2291 y[(i << 1)+1] = TMath::Cos((i << 1)+2)-.5;
2292 }
2293 printf("-N = %3d M =%2d\n",n,m);
2294 spline = new TSpline5("Test",x,y,nn);
2295 for (i = 0; i < nn; ++i)
2296 spline->GetCoeff(i,hx,a[i],a[i+200],a[i+400],
2297 a[i+600],a[i+800],a[i+1000]);
2298 delete spline;
2299 for (i = 0; i < mm1; ++i)
2300 diff[i] = com[i] = 0;
2301 for (k = 0; k < nn; ++k) {
2302 for (i = 0; i < mm; ++i)
2303 c[i] = a[k+i*200];
2304 if (n < 11) {
2305 printf(" ---------------------------------------%3d --------------------------------------------\n",k+1);
2306 printf("%12.8f\n",x[k]);
2307 if (k == nn-1) printf("%16.8f\n",c[0]);
2308 }
2309 if (k == nn-1) break;
2310 if (n <= 10) {
2311 for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
2312 printf("\n");
2313 }
2314 for (i = 0; i < mm1; ++i)
2315 if ((z=TMath::Abs(a[k+i*200])) > com[i]) com[i] = z;
2316 z = x[k+1]-x[k];
2317 for (i = 1; i < mm; ++i) {
2318 for (jj = i; jj < mm; ++jj) {
2319 j = mm+i-jj;
2320 c[j-2] = c[j-1]*z+c[j-2];
2321 }
2322 }
2323 if (n <= 10) {
2324 for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
2325 printf("\n");
2326 }
2327 for (i = 0; i < mm1; ++i)
2328 if (!(k >= nn-2 && i != 0))
2329 if ((z = TMath::Abs(c[i]-a[k+1+i*200]))
2330 > diff[i]) diff[i] = z;
2331 }
2332 printf(" MAXIMUM ABSOLUTE VALUES OF DIFFERENCES \n");
2333 for (i = 0; i < mm1; ++i) printf("%18.9E",diff[i]);
2334 printf("\n");
2335 printf(" MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS \n");
2336 if (TMath::Abs(c[0]) > com[0])
2337 com[0] = TMath::Abs(c[0]);
2338 for (i = 0; i < mm1; ++i) printf("%18.9E",com[i]);
2339 printf("\n");
2340 }
2341
2342 // test of TSpline5 with non equidistant knots, one double knot,
2343 // one triple knot, follows.
2344 printf("1 TEST OF TSpline5 WITH NONEQUIDISTANT KNOTS,\n");
2345 printf(" ONE DOUBLE, ONE TRIPLE KNOT\n");
2346 n = 8;
2347 x[0] = -3;
2348 x[1] = -1;
2349 x[2] = -1;
2350 x[3] = 0;
2351 x[4] = 3;
2352 x[5] = 3;
2353 x[6] = 3;
2354 x[7] = 4;
2355 y[0] = 7;
2356 y[1] = 11;
2357 y[2] = 15;
2358 y[3] = 26;
2359 y[4] = 56;
2360 y[5] = -30;
2361 y[6] = -7;
2362 y[7] = 29;
2363 m = 3;
2364 mm = m << 1;
2365 mm1 = mm-1;
2366 printf("-N = %3d M =%2d\n",n,m);
2367 spline=new TSpline5("Test",x,y,n);
2368 for (i = 0; i < n; ++i)
2369 spline->GetCoeff(i,hx,a[i],a[i+200],a[i+400],
2370 a[i+600],a[i+800],a[i+1000]);
2371 delete spline;
2372 for (i = 0; i < mm1; ++i)
2373 diff[i] = com[i] = 0;
2374 for (k = 0; k < n; ++k) {
2375 for (i = 0; i < mm; ++i)
2376 c[i] = a[k+i*200];
2377 printf(" ---------------------------------------%3d --------------------------------------------\n",k+1);
2378 printf("%12.8f\n",x[k]);
2379 if (k == n-1) {
2380 printf("%16.8f\n",c[0]);
2381 break;
2382 }
2383 for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
2384 printf("\n");
2385 for (i = 0; i < mm1; ++i)
2386 if ((z=TMath::Abs(a[k+i*200])) > com[i]) com[i] = z;
2387 z = x[k+1]-x[k];
2388 for (i = 1; i < mm; ++i)
2389 for (jj = i; jj < mm; ++jj) {
2390 j = mm+i-jj;
2391 c[j-2] = c[j-1]*z+c[j-2];
2392 }
2393 for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
2394 printf("\n");
2395 for (i = 0; i < mm1; ++i)
2396 if (!(k >= n-2 && i != 0))
2397 if ((z = TMath::Abs(c[i]-a[k+1+i*200]))
2398 > diff[i]) diff[i] = z;
2399 }
2400 printf(" MAXIMUM ABSOLUTE VALUES OF DIFFERENCES \n");
2401 for (i = 0; i < mm1; ++i) printf("%18.9E",diff[i]);
2402 printf("\n");
2403 printf(" MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS \n");
2404 if (TMath::Abs(c[0]) > com[0])
2405 com[0] = TMath::Abs(c[0]);
2406 for (i = 0; i < mm1; ++i) printf("%16.8f",com[i]);
2407 printf("\n");
2408
2409 // Test of TSpline5 with non equidistant knots, two double knots,
2410 // one triple knot,follows.
2411 printf("1 TEST OF TSpline5 WITH NONEQUIDISTANT KNOTS,\n");
2412 printf(" TWO DOUBLE, ONE TRIPLE KNOT\n");
2413 n = 10;
2414 x[0] = 0;
2415 x[1] = 2;
2416 x[2] = 2;
2417 x[3] = 3;
2418 x[4] = 3;
2419 x[5] = 3;
2420 x[6] = 5;
2421 x[7] = 8;
2422 x[8] = 9;
2423 x[9] = 9;
2424 y[0] = 163;
2425 y[1] = 237;
2426 y[2] = -127;
2427 y[3] = 119;
2428 y[4] = -65;
2429 y[5] = 192;
2430 y[6] = 293;
2431 y[7] = 326;
2432 y[8] = 0;
2433 y[9] = -414;
2434 m = 3;
2435 mm = m << 1;
2436 mm1 = mm-1;
2437 printf("-N = %3d M =%2d\n",n,m);
2438 spline = new TSpline5("Test",x,y,n);
2439 for (i = 0; i < n; ++i)
2440 spline->GetCoeff(i,hx,a[i],a[i+200],a[i+400],
2441 a[i+600],a[i+800],a[i+1000]);
2442 delete spline;
2443 for (i = 0; i < mm1; ++i)
2444 diff[i] = com[i] = 0;
2445 for (k = 0; k < n; ++k) {
2446 for (i = 0; i < mm; ++i)
2447 c[i] = a[k+i*200];
2448 printf(" ---------------------------------------%3d --------------------------------------------\n",k+1);
2449 printf("%12.8f\n",x[k]);
2450 if (k == n-1) {
2451 printf("%16.8f\n",c[0]);
2452 break;
2453 }
2454 for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
2455 printf("\n");
2456 for (i = 0; i < mm1; ++i)
2457 if ((z=TMath::Abs(a[k+i*200])) > com[i]) com[i] = z;
2458 z = x[k+1]-x[k];
2459 for (i = 1; i < mm; ++i)
2460 for (jj = i; jj < mm; ++jj) {
2461 j = mm+i-jj;
2462 c[j-2] = c[j-1]*z+c[j-2];
2463 }
2464 for (i = 0; i < mm; ++i) printf("%16.8f",c[i]);
2465 printf("\n");
2466 for (i = 0; i < mm1; ++i)
2467 if (!(k >= n-2 && i != 0))
2468 if((z = TMath::Abs(c[i]-a[k+1+i*200]))
2469 > diff[i]) diff[i] = z;
2470 }
2471 printf(" MAXIMUM ABSOLUTE VALUES OF DIFFERENCES \n");
2472 for (i = 0; i < mm1; ++i) printf("%18.9E",diff[i]);
2473 printf("\n");
2474 printf(" MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS \n");
2475 if (TMath::Abs(c[0]) > com[0])
2476 com[0] = TMath::Abs(c[0]);
2477 for (i = 0; i < mm1; ++i) printf("%16.8f",com[i]);
2478 printf("\n");
2479}
2480
2481////////////////////////////////////////////////////////////////////////////////
2482/// Stream an object of class TSpline5.
2483
2485{
2486 if (R__b.IsReading()) {
2487 UInt_t R__s, R__c;
2488 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
2489 if (R__v > 1) {
2490 R__b.ReadClassBuffer(TSpline5::Class(), this, R__v, R__s, R__c);
2491 return;
2492 }
2493 //====process old versions before automatic schema evolution
2495 if (fNp > 0) {
2496 fPoly = new TSplinePoly5[fNp];
2497 for(Int_t i=0; i<fNp; ++i) {
2498 fPoly[i].Streamer(R__b);
2499 }
2500 }
2501 // R__b >> fPoly;
2502 } else {
2503 R__b.WriteClassBuffer(TSpline5::Class(),this);
2504 }
2505}
#define d(i)
Definition RSha256.hxx:102
#define b(i)
Definition RSha256.hxx:100
#define f(i)
Definition RSha256.hxx:104
#define c(i)
Definition RSha256.hxx:101
#define g(i)
Definition RSha256.hxx:105
#define a(i)
Definition RSha256.hxx:99
#define h(i)
Definition RSha256.hxx:106
#define e(i)
Definition RSha256.hxx:103
int Int_t
Definition RtypesCore.h:45
short Version_t
Definition RtypesCore.h:65
constexpr Bool_t kFALSE
Definition RtypesCore.h:94
constexpr Bool_t kTRUE
Definition RtypesCore.h:93
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:374
#define X(type, name)
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t option
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char filename
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
float xmin
float * q
float xmax
@ kWritePermission
Definition TSystem.h:54
R__EXTERN TSystem * gSystem
Definition TSystem.h:572
#define gPad
#define snprintf
Definition civetweb.c:1540
Fill Area Attributes class.
Definition TAttFill.h:20
virtual void Streamer(TBuffer &)
virtual Color_t GetFillColor() const
Return the fill area color.
Definition TAttFill.h:31
virtual Style_t GetFillStyle() const
Return the fill area style.
Definition TAttFill.h:32
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition TAttFill.h:38
virtual void SetFillStyle(Style_t fstyle)
Set the fill area style.
Definition TAttFill.h:40
virtual void SaveFillAttributes(std::ostream &out, const char *name, Int_t coldef=1, Int_t stydef=1001)
Save fill attributes as C++ statement(s) on output stream out.
Definition TAttFill.cxx:239
Line Attributes class.
Definition TAttLine.h:20
virtual void Streamer(TBuffer &)
virtual Color_t GetLineColor() const
Return the line color.
Definition TAttLine.h:35
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition TAttLine.h:44
virtual Width_t GetLineWidth() const
Return the line width.
Definition TAttLine.h:37
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition TAttLine.h:45
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition TAttLine.h:42
virtual Style_t GetLineStyle() const
Return the line style.
Definition TAttLine.h:36
virtual void SaveLineAttributes(std::ostream &out, const char *name, Int_t coldef=1, Int_t stydef=1, Int_t widdef=1)
Save line attributes as C++ statement(s) on output stream out.
Definition TAttLine.cxx:275
Marker Attributes class.
Definition TAttMarker.h:20
virtual void SaveMarkerAttributes(std::ostream &out, const char *name, Int_t coldef=1, Int_t stydef=1, Int_t sizdef=1)
Save line attributes as C++ statement(s) on output stream out.
virtual Style_t GetMarkerStyle() const
Return the marker style.
Definition TAttMarker.h:33
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition TAttMarker.h:39
virtual Color_t GetMarkerColor() const
Return the marker color.
Definition TAttMarker.h:32
virtual Size_t GetMarkerSize() const
Return the marker size.
Definition TAttMarker.h:34
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition TAttMarker.h:41
virtual void Streamer(TBuffer &)
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
Definition TAttMarker.h:46
virtual void SetLimits(Double_t xmin, Double_t xmax)
Definition TAxis.h:166
Buffer base class used for serializing objects.
Definition TBuffer.h:43
1-Dim function class
Definition TF1.h:233
A TGraph is an object made of two arrays X and Y with npoints each.
Definition TGraph.h:41
void Paint(Option_t *chopt="") override
Draw this graph with its current attributes.
Definition TGraph.cxx:1978
1-D histogram with a float per channel (see TH1 documentation)
Definition TH1.h:647
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:59
virtual void SetDirectory(TDirectory *dir)
By default, when a histogram is created, it is added to the list of histogram objects in the current ...
Definition TH1.cxx:8923
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition TH1.cxx:9127
Int_t DistancetoPrimitive(Int_t px, Int_t py) override
Compute distance from point px,py to a line.
Definition TH1.cxx:2794
@ kLogX
X-axis in log scale.
Definition TH1.h:179
@ kNoStats
Don't draw stats box.
Definition TH1.h:176
TAxis * GetXaxis()
Definition TH1.h:341
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition TH1.cxx:9208
void Paint(Option_t *option="") override
Control routine to paint any kind of histograms.
Definition TH1.cxx:6204
void ExecuteEvent(Int_t event, Int_t px, Int_t py) override
Execute action corresponding to one event.
Definition TH1.cxx:3211
The TNamed class is the base class for all named ROOT classes.
Definition TNamed.h:29
const char * GetName() const override
Returns name of object.
Definition TNamed.h:49
void Streamer(TBuffer &) override
Stream an object of class TObject.
const char * GetTitle() const override
Returns title of object.
Definition TNamed.h:50
TString fName
Definition TNamed.h:32
TNamed & operator=(const TNamed &rhs)
TNamed assignment operator.
Definition TNamed.cxx:51
TObject & operator=(const TObject &rhs)
TObject assignment operator.
Definition TObject.h:300
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition TObject.h:203
virtual void AppendPad(Option_t *option="")
Append graphics object to current pad.
Definition TObject.cxx:203
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:847
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:1054
static void SavePrimitiveConstructor(std::ostream &out, TClass *cl, const char *variable_name, const char *constructor_agrs="", Bool_t empty_line=kTRUE)
Save object constructor in the output stream "out".
Definition TObject.cxx:771
Class to create third splines to interpolate knots Arbitrary conditions can be introduced for first a...
Definition TSpline.h:182
void SaveAs(const char *filename="", Option_t *option="") const override
Write this spline as a C++ function that can be executed without ROOT the name of the function is the...
Definition TSpline.cxx:807
Int_t fEndCond
0=no end cond, 1=first derivative, 2=second derivative
Definition TSpline.h:188
Int_t fBegCond
0=no beg cond, 1=first derivative, 2=second derivative
Definition TSpline.h:187
Int_t FindX(Double_t x) const
Find X.
Definition TSpline.cxx:744
void SavePrimitive(std::ostream &out, Option_t *option="") override
Save primitive as a C++ statement(s) on output stream out.
Definition TSpline.cxx:966
static void Test()
Test method for TSpline5.
Definition TSpline.cxx:613
Double_t fValBeg
Initial value of first or second derivative.
Definition TSpline.h:185
void BuildCoeff() override
Build coefficients.
Definition TSpline.cxx:1045
void Streamer(TBuffer &) override
Stream an object of class TSpline3.
Definition TSpline.cxx:1170
Double_t Eval(Double_t x) const override
Eval this spline at x.
Definition TSpline.cxx:786
Double_t fValEnd
End value of first or second derivative.
Definition TSpline.h:186
void SetCond(const char *opt)
Check the boundary conditions.
Definition TSpline.cxx:574
Double_t Derivative(Double_t x) const
Derivative.
Definition TSpline.cxx:796
TSplinePoly3 * fPoly
[fNp] Array of polynomial terms
Definition TSpline.h:184
TSpline3 & operator=(const TSpline3 &)
Assignment operator.
Definition TSpline.cxx:550
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set point number i.
Definition TSpline.cxx:992
TSpline3()
Definition TSpline.h:194
virtual void SetPointCoeff(Int_t i, Double_t b, Double_t c, Double_t d)
Set point coefficient number i.
Definition TSpline.cxx:1002
static TClass * Class()
Class to create quintic natural splines to interpolate knots Arbitrary conditions can be introduced f...
Definition TSpline.h:238
static void Test()
Test method for TSpline5.
Definition TSpline.cxx:2079
Double_t Eval(Double_t x) const override
Eval this spline at x.
Definition TSpline.cxx:1560
static TClass * Class()
void Streamer(TBuffer &) override
Stream an object of class TSpline5.
Definition TSpline.cxx:2484
void BuildCoeff() override
Algorithm 600, collected algorithms from acm.
Definition TSpline.cxx:1901
TSplinePoly5 * fPoly
[fNp] Array of polynomial terms
Definition TSpline.h:240
virtual void SetPointCoeff(Int_t i, Double_t b, Double_t c, Double_t d, Double_t e, Double_t f)
Set point coefficient number i.
Definition TSpline.cxx:1813
void SaveAs(const char *filename="", Option_t *option="") const override
Write this spline as a C++ function that can be executed without ROOT the name of the function is the...
Definition TSpline.cxx:1579
TSpline5()
Definition TSpline.h:250
Double_t Derivative(Double_t x) const
Derivative.
Definition TSpline.cxx:1569
void BoundaryConditions(const char *opt, Int_t &beg, Int_t &end, const char *&cb1, const char *&ce1, const char *&cb2, const char *&ce2)
Check the boundary conditions and the amount of extra double knots needed.
Definition TSpline.cxx:1444
void SavePrimitive(std::ostream &out, Option_t *option="") override
Save primitive as a C++ statement(s) on output stream out.
Definition TSpline.cxx:1770
void 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)
Set the boundary conditions at double/triple knots.
Definition TSpline.cxx:1475
Int_t FindX(Double_t x) const
Find X.
Definition TSpline.cxx:1524
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set point number i.
Definition TSpline.cxx:1802
TSpline5 & operator=(const TSpline5 &)
Assignment operator.
Definition TSpline.cxx:1424
Class for TSpline3 knot.
Definition TSpline.h:105
Double_t fC
Second order expansion coefficient : fC*2! is the second derivative at x.
Definition TSpline.h:108
Double_t Eval(Double_t x) const override
Definition TSpline.h:121
Double_t & D()
Definition TSpline.h:120
Double_t fD
Third order expansion coefficient : fD*3! is the third derivative at x.
Definition TSpline.h:109
Double_t Derivative(Double_t x) const
Definition TSpline.h:126
Double_t & C()
Definition TSpline.h:119
Double_t fB
First order expansion coefficient : fB*1! is the first derivative at x.
Definition TSpline.h:107
Double_t & B()
Definition TSpline.h:118
void CopyPoly(TSplinePoly3 const &other)
Utility called by the copy constructors and = operator.
Definition TSpline.cxx:311
TSplinePoly3 & operator=(TSplinePoly3 const &other)
Assignment operator.
Definition TSpline.cxx:299
void Streamer(TBuffer &) override
Stream an object of class TObject.
Class for TSpline5 knot.
Definition TSpline.h:140
Double_t & C()
Definition TSpline.h:159
void CopyPoly(TSplinePoly5 const &other)
Utility called by the copy constructors and = operator.
Definition TSpline.cxx:338
TSplinePoly5 & operator=(TSplinePoly5 const &other)
Assignment operator.
Definition TSpline.cxx:326
Double_t fF
Fifth order expansion coefficient : fF*5! is the fifth derivative at x.
Definition TSpline.h:146
Double_t Derivative(Double_t x) const
Definition TSpline.h:168
void Streamer(TBuffer &) override
Stream an object of class TObject.
Double_t fB
First order expansion coefficient : fB*1! is the first derivative at x.
Definition TSpline.h:142
Double_t fC
Second order expansion coefficient : fC*2! is the second derivative at x.
Definition TSpline.h:143
Double_t & B()
Definition TSpline.h:158
Double_t & F()
Definition TSpline.h:162
Double_t fD
Third order expansion coefficient : fD*3! is the third derivative at x.
Definition TSpline.h:144
Double_t & E()
Definition TSpline.h:161
Double_t Eval(Double_t x) const override
Definition TSpline.h:163
Double_t fE
Fourth order expansion coefficient : fE*4! is the fourth derivative at x.
Definition TSpline.h:145
Double_t & D()
Definition TSpline.h:160
Base class for TSpline knot.
Definition TSpline.h:75
Double_t & Y()
Definition TSpline.h:88
Double_t & X()
Definition TSpline.h:87
void CopyPoly(TSplinePoly const &other)
Utility called by the copy constructors and = operator.
Definition TSpline.cxx:285
TSplinePoly & operator=(TSplinePoly const &other)
Assignment operator.
Definition TSpline.cxx:273
Double_t fY
Constant term.
Definition TSpline.h:78
Double_t fX
Abscissa.
Definition TSpline.h:77
Base class for spline implementation containing the Draw/Paint methods.
Definition TSpline.h:31
void Streamer(TBuffer &) override
Stream an object of class TSpline.
Definition TSpline.cxx:231
virtual Double_t Eval(Double_t x) const =0
void ExecuteEvent(Int_t event, Int_t px, Int_t py) override
Execute action corresponding to one event.
Definition TSpline.cxx:122
TGraph * fGraph
Graph for drawing the knots.
Definition TSpline.h:39
Int_t DistancetoPrimitive(Int_t px, Int_t py) override
Compute distance from point px,py to a spline.
Definition TSpline.cxx:113
Double_t fXmin
Minimum value of abscissa.
Definition TSpline.h:34
TClass * IsA() const override
Definition TSpline.h:69
TH1F * fHistogram
Temporary histogram.
Definition TSpline.h:38
TSpline()
Definition TSpline.h:47
void Draw(Option_t *option="") override
Draw this function with its current attributes.
Definition TSpline.cxx:101
Double_t fDelta
Distance between equidistant knots.
Definition TSpline.h:33
virtual void GetKnot(Int_t i, Double_t &x, Double_t &y) const =0
TSpline & operator=(const TSpline &)
Assignment operator.
Definition TSpline.cxx:69
void Paint(Option_t *option="") override
Paint this function with its current attributes.
Definition TSpline.cxx:131
Bool_t fKstep
True of equidistant knots.
Definition TSpline.h:37
Int_t fNp
Number of knots.
Definition TSpline.h:36
~TSpline() override
Destructor.
Definition TSpline.cxx:60
static TClass * Class()
Double_t fXmax
Maximum value of abscissa.
Definition TSpline.h:35
Int_t fNpx
Number of points used for graphical representation.
Definition TSpline.h:40
Basic string class.
Definition TString.h:139
Ssiz_t Length() const
Definition TString.h:417
void ToLower()
Change string to lower-case.
Definition TString.cxx:1182
TString & ReplaceSpecialCppChars()
Find special characters which are typically used in printf() calls and replace them by appropriate es...
Definition TString.cxx:1114
const char * Data() const
Definition TString.h:376
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
Definition TString.cxx:2378
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:632
virtual Bool_t AccessPathName(const char *path, EAccessMode mode=kFileExists)
Returns FALSE if one can access a file using the specified access mode.
Definition TSystem.cxx:1296
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
#define F(x, y, z)
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:250
Int_t FloorNint(Double_t x)
Returns the nearest integer of TMath::Floor(x).
Definition TMath.h:690
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:198
Double_t Cos(Double_t)
Returns the cosine of an angle of x radians.
Definition TMath.h:598
Double_t Sin(Double_t)
Returns the sine of an angle of x radians.
Definition TMath.h:592
Double_t Log10(Double_t x)
Returns the common (base-10) logarithm of x.
Definition TMath.h:766
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123
TMarker m
Definition textangle.C:8
TLine l
Definition textangle.C:4