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