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