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