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