Logo ROOT  
masterv6.34v6.32v6.30v6.28v6.26v6.24v6.22v6.20v6.18v6.16v6.14v6.12v6.10v6.08v6.06v6.04

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