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