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