Logo ROOT  

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