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