Logo ROOT  

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