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