MaCh3  2.5.0
Reference Guide
Public Member Functions | Protected Attributes | List of all members
TSpline3_red Class Reference

CW: Reduced TSpline3 class. More...

#include <Splines/SplineStructs.h>

Inheritance diagram for TSpline3_red:
[legend]
Collaboration diagram for TSpline3_red:
[legend]

Public Member Functions

 TSpline3_red ()
 Empty constructor. More...
 
 TSpline3_red (TSpline3 *&spline, SplineInterpolation InterPolation=kTSpline3)
 The constructor that takes a TSpline3 pointer and copies in to memory. More...
 
 TSpline3_red (M3::float_t *X, M3::float_t *Y, M3::int_t N, M3::float_t **P)
 constructor taking parameters More...
 
void SetFunc (TSpline3 *&spline, SplineInterpolation InterPolation=kTSpline3)
 Set the function [28]. More...
 
virtual ~TSpline3_red ()
 Empty destructor. More...
 
int FindX (double x)
 Find the segment relevant to this variation in x. More...
 
double Eval (double var) override
 CW: Evaluate the weight from a variation. More...
 
M3::int_t GetNp () override
 CW: Get the number of points. More...
 
void GetKnot (int i, M3::float_t &xtmp, M3::float_t &ytmp)
 
void GetCoeff (int segment, M3::float_t &x, M3::float_t &y, M3::float_t &b, M3::float_t &c, M3::float_t &d)
 CW: Get the coefficient of a given segment. More...
 
TSpline3 * ConstructTSpline3 ()
 CW: Make a TSpline3 from the reduced splines. More...
 
void Print () override
 Print detailed info. More...
 
- Public Member Functions inherited from TResponseFunction_red
 TResponseFunction_red ()
 Empty constructor. More...
 
virtual ~TResponseFunction_red ()
 Empty destructor. More...
 

Protected Attributes

M3::int_t nPoints
 Number of points/knot in TSpline3. More...
 
M3::float_t ** Par
 Always uses a third order polynomial, so hard-code the number of coefficients in implementation. More...
 
M3::float_tXPos
 Positions of each x for each knot. More...
 
M3::float_tYResp
 y-value for each knot More...
 

Detailed Description

CW: Reduced TSpline3 class.

Definition at line 266 of file SplineStructs.h.

Constructor & Destructor Documentation

◆ TSpline3_red() [1/3]

TSpline3_red::TSpline3_red ( )
inline

Empty constructor.

Definition at line 270 of file SplineStructs.h.

271  nPoints = 0;
272  Par = nullptr;
273  XPos = nullptr;
274  YResp = nullptr;
275  }
TResponseFunction_red()
Empty constructor.
M3::float_t ** Par
Always uses a third order polynomial, so hard-code the number of coefficients in implementation.
M3::int_t nPoints
Number of points/knot in TSpline3.
M3::float_t * XPos
Positions of each x for each knot.
M3::float_t * YResp
y-value for each knot

◆ TSpline3_red() [2/3]

TSpline3_red::TSpline3_red ( TSpline3 *&  spline,
SplineInterpolation  InterPolation = kTSpline3 
)
inline

The constructor that takes a TSpline3 pointer and copies in to memory.

Definition at line 278 of file SplineStructs.h.

279  Par = nullptr;
280  XPos = nullptr;
281  YResp = nullptr;
282  SetFunc(spline, InterPolation);
283  }
void SetFunc(TSpline3 *&spline, SplineInterpolation InterPolation=kTSpline3)
Set the function .

◆ TSpline3_red() [3/3]

TSpline3_red::TSpline3_red ( M3::float_t X,
M3::float_t Y,
M3::int_t  N,
M3::float_t **  P 
)
inline

constructor taking parameters

Definition at line 286 of file SplineStructs.h.

287  nPoints = N;
288  // Save the parameters for each knot
289  Par = new M3::float_t*[nPoints];
290  // Save the positions of the knots
291  XPos = new M3::float_t[nPoints];
292  // Save the y response at each knot
293  YResp = new M3::float_t[nPoints];
294  for(int j = 0; j < N; ++j){
295  Par[j] = new M3::float_t[3];
296  Par[j][0] = P[j][0];
297  Par[j][1] = P[j][1];
298  Par[j][2] = P[j][2];
299  XPos[j] = X[j];
300  YResp[j] = Y[j];
301 
302  if((Par[j][0] == -999) | (Par[j][1] ==-999) | (Par[j][2] ==-999) | (XPos[j] ==-999) | (YResp[j] ==-999)){
303  MACH3LOG_ERROR("******************* Bad parameter values when constructing TSpline3_red *********************");
304  MACH3LOG_ERROR("Passed values (i, x, y, b, c, d): {}, {}, {}, {}, {}, {}", j, X[j], Y[j], P[j][0], P[j][1], P[j][2]);
305  MACH3LOG_ERROR("Set values (i, x, y, b, c, d): {}, {}, {}, {}, {}, {}", j, XPos[j], YResp[j], Par[j][0], Par[j][1], Par[j][2]);
306  MACH3LOG_ERROR("*********************************************************************************************");
307  }
308  }
309  }
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:37
double float_t
Definition: Core.h:37

◆ ~TSpline3_red()

virtual TSpline3_red::~TSpline3_red ( )
inlinevirtual

Empty destructor.

Definition at line 656 of file SplineStructs.h.

656  {
657  if(Par != nullptr) {
658  for (int i = 0; i < nPoints; ++i) {
659  if (Par[i] != nullptr) {
660  delete[] Par[i];
661  }
662  }
663  delete[] Par;
664  }
665  if(XPos != nullptr) delete[] XPos;
666  if(YResp != nullptr) delete[] YResp;
667  Par = nullptr;
668  XPos = YResp = nullptr;
669  }

Member Function Documentation

◆ ConstructTSpline3()

TSpline3* TSpline3_red::ConstructTSpline3 ( )
inline

CW: Make a TSpline3 from the reduced splines.

Definition at line 738 of file SplineStructs.h.

738  {
739  // KS: Sadly ROOT only accepts double...
740  #ifdef _LOW_MEMORY_STRUCTS_
741  std::vector<Double_t> xPosDoubles(nPoints);
742  std::vector<Double_t> yPosDoubles(nPoints);
743  for (Int_t i = 0; i < nPoints; ++i) {
744  xPosDoubles[i] = static_cast<Double_t>(XPos[i]); // Convert float to double
745  yPosDoubles[i] = static_cast<Double_t>(YResp[i]); // Convert float to double
746  }
747  TSpline3 *spline = new TSpline3("Spline", xPosDoubles.data(), yPosDoubles.data(), static_cast<int>(nPoints));
748  #else
749  TSpline3 *spline = new TSpline3("Spline", XPos, YResp, nPoints);
750  #endif
751  for (Int_t i = 0; i < nPoints; ++i) {
752  spline->SetPointCoeff(i, Par[i][0], Par[i][1], Par[i][2]);
753  }
754 
755  return spline;
756  }

◆ Eval()

double TSpline3_red::Eval ( double  var)
inlineoverridevirtual

CW: Evaluate the weight from a variation.

Implements TResponseFunction_red.

Definition at line 708 of file SplineStructs.h.

708  {
709  // Get the segment for this variation
710  int segment = FindX(var);
711  // The get the coefficients for this variation
712  M3::float_t x = M3::float_t(-999.99), y = M3::float_t(-999.99), b = M3::float_t(-999.99), c = M3::float_t(-999.99), d = M3::float_t(-999.99);
713  GetCoeff(segment, x, y, b, c, d);
714  double dx = var - x;
715  // Evaluate the third order polynomial
716  double weight = y+dx*(b+dx*(c+d*dx));
717  return weight;
718  }
int FindX(double x)
Find the segment relevant to this variation in x.
void GetCoeff(int segment, M3::float_t &x, M3::float_t &y, M3::float_t &b, M3::float_t &c, M3::float_t &d)
CW: Get the coefficient of a given segment.

◆ FindX()

int TSpline3_red::FindX ( double  x)
inline

Find the segment relevant to this variation in x.

See also
TSpline3::FindX(double) for ROOT implementation
SplineBase::FindSplineSegment FindSplineSegment for the base class version

Definition at line 674 of file SplineStructs.h.

674  {
675  // The segment we're interested in (klow in ROOT code)
676  int segment = 0;
677  int kHigh = nPoints-1;
678  // If the variation is below the lowest saved spline point
679  if (x <= XPos[0]){
680  segment = 0;
681  // If the variation is above the highest saved spline point
682  } else if (x >= XPos[nPoints-1]) {
683  // Yes, the -2 is indeed correct, see TSpline.cxx:814 and //see: https://savannah.cern.ch/bugs/?71651
684  segment = kHigh;
685  // If the variation is between the maximum and minimum, perform a binary search
686  } else {
687  // The top point we've got
688  int kHalf = 0;
689  // While there is still a difference in the points (we haven't yet found the segment)
690  // This is a binary search, incrementing segment and decrementing kHalf until we've found the segment
691  while (kHigh - segment > 1) {
692  // Increment the half-step
693  kHalf = (segment + kHigh)/2;
694  // If our variation is above the kHalf, set the segment to kHalf
695  if (x > XPos[kHalf]) {
696  segment = kHalf;
697  // Else move kHigh down
698  } else {
699  kHigh = kHalf;
700  }
701  } // End the while: we've now done our binary search
702  } // End the else: we've now found our point
703  if (segment >= nPoints-1 && nPoints > 1) segment = nPoints-2;
704  return segment;
705  }

◆ GetCoeff()

void TSpline3_red::GetCoeff ( int  segment,
M3::float_t x,
M3::float_t y,
M3::float_t b,
M3::float_t c,
M3::float_t d 
)
inline

CW: Get the coefficient of a given segment.

Definition at line 729 of file SplineStructs.h.

729  {
730  b = Par[segment][0];
731  c = Par[segment][1];
732  d = Par[segment][2];
733  x = XPos[segment];
734  y = YResp[segment];
735  }

◆ GetKnot()

void TSpline3_red::GetKnot ( int  i,
M3::float_t xtmp,
M3::float_t ytmp 
)
inline

Definition at line 723 of file SplineStructs.h.

723  {
724  xtmp = XPos[i];
725  ytmp = YResp[i];
726  }

◆ GetNp()

M3::int_t TSpline3_red::GetNp ( )
inlineoverridevirtual

CW: Get the number of points.

Implements TResponseFunction_red.

Definition at line 721 of file SplineStructs.h.

721 { return nPoints; }

◆ Print()

void TSpline3_red::Print ( )
inlineoverridevirtual

Print detailed info.

Implements TResponseFunction_red.

Definition at line 759 of file SplineStructs.h.

759  {
760  MACH3LOG_INFO("Printing TSpline_red:");
761  MACH3LOG_INFO(" Nknots = {}", nPoints);
762  for (int i = 0; i < nPoints; ++i) {
763  MACH3LOG_INFO(" i = {} x = {} y = {} b = {} c = {} d = {}",
764  i, XPos[i], YResp[i], Par[i][0], Par[i][1], Par[i][2]);
765  }
766  }
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:35

◆ SetFunc()

void TSpline3_red::SetFunc ( TSpline3 *&  spline,
SplineInterpolation  InterPolation = kTSpline3 
)
inline

Set the function [28].

Definition at line 312 of file SplineStructs.h.

312  {
313  nPoints = M3::int_t(spline->GetNp());
314  if (Par != nullptr) {
315  for (int i = 0; i < nPoints; ++i) {
316  delete[] Par[i];
317  Par[i] = nullptr;
318  }
319  delete[] Par;
320  Par = nullptr;
321  }
322  if (XPos != nullptr) delete[] XPos;
323  if (YResp != nullptr) delete[] YResp;
324  // Save the parameters for each knot
325  Par = new M3::float_t*[nPoints];
326  // Save the positions of the knots
327  XPos = new M3::float_t[nPoints];
328  // Save the y response at each knot
329  YResp = new M3::float_t[nPoints];
330 
331  //KS: Default TSpline3 ROOT implementation
332  if(InterPolation == kTSpline3)
333  {
334  for (int i = 0; i < nPoints; ++i) {
335  // 3 is the size of the TSpline3 coefficients
336  Par[i] = new M3::float_t[3];
337  double x = -999.99, y = -999.99, b = -999.99, c = -999.99, d = -999.99;
338  spline->GetCoeff(i, x, y, b, c, d);
339  XPos[i] = M3::float_t(x);
340  YResp[i] = M3::float_t(y);
341  Par[i][0] = M3::float_t(b);
342  Par[i][1] = M3::float_t(c);
343  Par[i][2] = M3::float_t(d);
344  }
345  }
346  //CW: Reduce to use linear spline interpolation for certain parameters
347  // Not the most elegant way: use TSpline3 object but set coefficients to zero and recalculate spline points; the smart way (but more human intensive) would be to save memory here and simply not store the zeros at all
348  // Get which parameters should be linear from the fit manager
349  // Convert the spline number to global xsec parameter
350  // Loop over the splines points
351  // KS: kLinearFunc should be used with TF1, this is just as safety
352  else if(InterPolation == kLinear || InterPolation == kLinearFunc)
353  {
354  for (int k = 0; k < nPoints; ++k) {
355  Par[k] = new M3::float_t[3];
356 
357  Double_t x1, y1, b1, c1, d1, x2, y2, b2, c2, d2 = 0;
358  spline->GetCoeff(k, x1, y1, b1, c1, d1);
359 
360  double tempb = 0;
361  if (k == nPoints - 1) {
362  tempb = Par[k-1][0];
363  } else {
364  spline->GetCoeff(k + 1, x2, y2, b2, c2, d2);
365  tempb = (y2-y1)/(x2-x1);
366  }
367  XPos[k] = M3::float_t(x1);
368  YResp[k] = M3::float_t(y1);
369  Par[k][0] = M3::float_t(tempb); // linear slope
370  Par[k][1] = M3::float_t(0);
371  Par[k][2] = M3::float_t(0);
372  }
373  }
374  //EM: Akima spline is similar to regular cubic spline but is allowed to be discontinuous in 2nd derivative and coefficients in any segment
375  // only depend on th 2 nearest points on either side
376  else if(InterPolation == kAkima)
377  {
378  // get the knot values for the spline
379  for (int i = 0; i < nPoints; ++i) {
380  // 3 is the size of the TSpline3 coefficients
381  Par[i] = new M3::float_t[3];
382 
383  double x = -999.99, y = -999.99;
384  spline->GetKnot(i, x, y);
385 
386  XPos[i] = M3::float_t(x);
387  YResp[i] = M3::float_t(y);
388  }
389 
390  M3::float_t* mvals = new M3::float_t[nPoints + 3];
391  M3::float_t* svals = new M3::float_t[nPoints + 1];
392 
393  for (int i = -2; i <= nPoints; ++i) {
394  // if segment is first or last or 2nd to first or last, needs to be dealt with slightly differently;
395  // need to estimate the values for additinal points which would lie outside of the spline
396  if(i ==-2){
397  mvals[i+2] = M3::float_t(3.0 * (YResp[1] - YResp[0]) / (XPos[1] - XPos[0]) - 2.0*(YResp[2] - YResp[1]) / (XPos[2] - XPos[1]));
398  }
399  else if(i==-1){
400  mvals[i+2] = M3::float_t(2.0 * (YResp[1] - YResp[0]) / (XPos[1] - XPos[0]) - (YResp[2] - YResp[1]) / (XPos[2] - XPos[1]));
401  }
402  else if(i==nPoints){
403  mvals[i+2] = M3::float_t(3.0 * (YResp[nPoints-1] - YResp[nPoints-2]) / (XPos[nPoints-1] - XPos[nPoints-2]) - 2.0*(YResp[nPoints-2] - YResp[nPoints-3]) / (XPos[nPoints-2] - XPos[nPoints-3]));
404  }
405  else if(i == nPoints - 1){
406  mvals[i+2] = M3::float_t(2.0 * (YResp[nPoints-1] - YResp[nPoints-2]) / (XPos[nPoints-1] - XPos[nPoints-2]) - (YResp[nPoints-2] - YResp[nPoints-3]) / (XPos[nPoints-2] - XPos[nPoints-3]));
407  }
408  //standard internal segment
409  else{
410  mvals[i+2] = (YResp[i+1] - YResp[i])/ (XPos[i+1] - XPos[i]);
411  }
412  }
413 
414  for(int i =2; i<=nPoints+2; i++){
415  if (std::abs(mvals[i+1] - mvals[i]) + std::abs(mvals[i-1] - mvals[i-2]) != 0.0){
416  svals[i-2] = (std::abs(mvals[i+1] - mvals[i]) * mvals[i-1] + std::abs(mvals[i-1] - mvals[i-2]) *mvals[i]) / (std::abs(mvals[i+1] - mvals[i]) + std::abs(mvals[i-1] - mvals[i-2]));
417  }
418  else{svals[i-2] = mvals[i];}
419  }
420 
421  // calculate the coefficients for the spline
422  for(int i = 0; i <nPoints; i++){
423  M3::float_t b, c, d = M3::float_t(-999.999);
424 
425  b = svals[i];
426  c = M3::float_t(3.0* (YResp[i+1] - YResp[i]) / (XPos[i+1] - XPos[i]) -2.0 *svals[i] - svals[i +1]) /(XPos[i+1] - XPos[i]);
427  d = M3::float_t((svals[i + 1] +svals[i]) - 2.0*(YResp[i+1] - YResp[i]) / (XPos[i+1] - XPos[i])) / ((XPos[i+1] - XPos[i]) * (XPos[i+1] - XPos[i]));
428 
429  Par[i][0] = b;
430  Par[i][1] = c;
431  Par[i][2] = d;
432  }
433 
434  // check the input spline for linear segments, if there are any then overwrite the calculated coefficients
435  // this will pretty much only ever be the case if they are set to be linear in samplePDFND i.e. the user wants it to be linear
436  for(int i = 0; i <nPoints-1; i++){
437  double x = -999.99, y = -999.99, b = -999.99, c = -999.99, d = -999.99;
438  spline->GetCoeff(i, x, y, b, c, d);
439 
440  if((c == 0.0 && d == 0.0)){
441  Par[i][0] = M3::float_t(b);
442  Par[i][1] = M3::float_t(0.0);
443  Par[i][2] = M3::float_t(0.0);
444  }
445  }
446  delete[] mvals;
447  delete[] svals;
448  }
449  //EM: Monotone spline is similar to regular cubic spline but enforce the condition that the interpolated value at any point
450  // must be between its two nearest knots, DOES NOT make the entire spline monotonic, only the segments
451  else if(InterPolation == kMonotonic)
452  {
453  // values of the secants at each point (for calculating monotone spline)
454  M3::float_t * Secants = new M3::float_t[nPoints -1];
455  // values of the tangens at each point (for calculating monotone spline)
456  M3::float_t * Tangents = new M3::float_t[nPoints];
457 
458  // get the knot values for the spline
459  for (int i = 0; i < nPoints; ++i) {
460  // 3 is the size of the TSpline3 coefficients
461  Par[i] = new M3::float_t[3];
462 
463  double x = -999.99, y = -999.99;
464  spline->GetKnot(i, x, y);
465 
466  XPos[i] = M3::float_t(x);
467  YResp[i] = M3::float_t(y);
468 
469  Tangents[i] = 0.0;
470  }
471 
472  // deal with the case of two points (just do linear interpolation between them)
473  if (nPoints == 2){
474  Par[0][0] = (YResp[1] - YResp[0]) / (XPos[1] - XPos[0]);
475  Par[0][1] = 0.0;
476  Par[0][2] = 0.0;
477  // extra "virtual" segment at end to make Par array shape fit with knot arrays shapes
478  Par[1][1] = 0.0;
479  Par[1][2] = 0.0;
480 
481  return;
482  } // if nPoints !=2 do full monotonic spline treatment:
483  else
484  {
485  // first pass over knots to calculate the secants
486  for (int i = 0; i < nPoints-1; ++i) {
487  Secants[i] = (YResp[i+1] - YResp[i]) / (XPos[i+1] - XPos[i]);
488  MACH3LOG_TRACE("Secant {}: {}", i, Secants[i]);
489  }
490 
491  Tangents[0] = Secants[0];
492  Tangents[nPoints-1] = Secants[nPoints -2];
493 
494  M3::float_t alpha;
495  M3::float_t beta;
496 
497  // second pass over knots to calculate tangents
498  for (int i = 1; i < nPoints-1; ++i) {
499  if ((Secants[i-1] >= 0.0 && Secants[i] >= 0.0) | (Secants[i-1] < 0.0 && Secants[i] < 0.0)){ //check for same sign
500  Tangents[i] = M3::float_t((Secants[i-1] + Secants[i]) /2.0);
501  }
502  }
503 
504  // third pass over knots to rescale tangents
505  for (int i = 0; i < nPoints-1; ++i) {
506  if (Secants[i] == 0.0){
507  Tangents[i] = 0.0;
508  Tangents[i+1] = 0.0;
509  }
510 
511  else{
512  alpha = Tangents[i] / Secants[i];
513  beta = Tangents[i+1] / Secants[i];
514 
515  if (alpha <0.0){
516  Tangents[i] = 0.0;
517  }
518  if (beta < 0.0){
519  Tangents[i+1] = 0.0;
520  }
521 
522  if (alpha * alpha + beta * beta >9.0){
523  M3::float_t tau = M3::float_t(3.0 / std::sqrt(alpha * alpha + beta * beta));
524  Tangents[i] = tau * alpha * Secants[i];
525  Tangents[i+1] = tau * beta * Secants[i];
526  }
527  }
528  } // finished rescaling tangents
529  // fourth pass over knots to calculate the coefficients for the spline
530  M3::float_t dx;
531  for(int i = 0; i <nPoints-1; i++){
532  M3::float_t b, c, d = M3::float_t(-999.999);
533  dx = XPos[i+1] - XPos[i];
534 
535  b = Tangents[i] * dx;
536  c = M3::float_t(3.0* (YResp[i+1] - YResp[i]) -2.0 *dx * Tangents[i] - dx * Tangents[i +1]);
537  d = M3::float_t(2.0* (YResp[i] - YResp[i+1]) + dx * (Tangents[i] + Tangents[i+1]));
538 
539  Par[i][0] = b / dx;
540  Par[i][1] = c / (dx * dx);
541  Par[i][2] = d / (dx * dx * dx);
542 
543  if((Par[i][0] == -999) | (Par[i][1] == -999) | (Par[i][2] ==-999) | (Par[i][0] == -999.999) | (Par[i][1] == -999.999) | (Par[i][2] ==-999.999)){
544  MACH3LOG_INFO("Bad spline parameters for segment {}: (b, c, d) = {}, {}, {}. This will cause problems with GPU.",
545  i, Par[i][0], Par[i][1], Par[i][2]);
546  }
547  MACH3LOG_TRACE("b: {}", b);
548  MACH3LOG_TRACE("dx: {}, x_0: {}, x_1: {}", dx, XPos[i], XPos[i+1]);
549  MACH3LOG_TRACE(" y_0: {}, y_1: {}", YResp[i], YResp[i+1]);
550  }
551 
552  // include params for final "segment" outside of the spline so that par array fits with x and y arrays,
553  // should never actually get used but if not set then the GPU code gets very angry
554  Par[nPoints-1][0] = 0.0;
555  Par[nPoints-1][1] = 0.0;
556  Par[nPoints-1][2] = 0.0;
557 
558  // check the input spline for linear segments, if there are any then overwrite the calculated coefficients
559  // this will pretty much only ever be the case if they are set to be linear in samplePDFND i.e. the user wants it to be linear
560  for(int i = 0; i <nPoints-1; i++){
561  double x = -999.99, y = -999.99, b = -999.99, c = -999.99, d = -999.99;
562  spline->GetCoeff(i, x, y, b, c, d);
563 
564  if((c == 0.0 && d == 0.0)){
565  Par[i][0] = M3::float_t(b);
566  Par[i][1] = 0.0;
567  Par[i][2] = 0.0;
568  }
569  }
570  delete[] Secants;
571  delete[] Tangents;
572  } // end of if(nPoints !=2)
573  }
574  else if (InterPolation == kKochanekBartels)
575  {
576  // Allocate memory for tangents and coefficients
577  M3::float_t* Tangents = new M3::float_t[nPoints];
578  for (int i = 0; i < nPoints; ++i)
579  {
580  Par[i] = new M3::float_t[3];
581  double x = -999.99, y = -999.99;
582  spline->GetKnot(i, x, y);
583  XPos[i] = M3::float_t(x);
584  YResp[i] = M3::float_t(y);
585  Tangents[i] = 0.0;
586  }
587 
588  // KS: Setting all parameters to 0 gives a Catmull-Rom spline.
589  // Tension (T):
590  // - T = 0.0: Default smooth interpolation (Catmull-Rom).
591  // - T = 1.0: Maximum tension; curve becomes linear between knots (no curvature).
592  // - T = -1.0: Overshooting; creates loops or "bouncy" effects around knots.
593  constexpr M3::float_t T = 0.0;
594 
595  // Continuity (C):
596  // - C = 0.0: Default smooth transition (continuous first derivative).
597  // - C = 1.0: Sharp corner/crease at knot (discontinuous first derivative).
598  // - C = -1.0: Inverted curve; creates loops or kinks at knot.
599  constexpr M3::float_t C = 0.0;
600 
601  // Bias (B):
602  // - B = 0.0: Default symmetric curve around knot.
603  // - B = 1.0: Curve is "pulled" toward the next knot (leading effect).
604  // - B = -1.0: Curve is "pulled" toward the previous knot (lagging effect).
605  constexpr M3::float_t B = 0.0;
606 
607  // Calculate tangents for internal knots
608  for (int i = 1; i < nPoints - 1; ++i)
609  {
610  M3::float_t d0 = (YResp[i] - YResp[i-1]) / (XPos[i] - XPos[i-1]);
611  M3::float_t d1 = (YResp[i+1] - YResp[i]) / (XPos[i+1] - XPos[i]);
612 
613  M3::float_t term1 = (1.0 - T) * (1.0 + B) * (1.0 + C) * 0.5 * d0;
614  M3::float_t term2 = (1.0 - T) * (1.0 - B) * (1.0 - C) * 0.5 * d1;
615 
616  Tangents[i] = term1 + term2;
617  }
618 
619  // Boundary conditions (simple choice: secant slopes)
620  Tangents[0] = (YResp[1] - YResp[0]) / (XPos[1] - XPos[0]);
621  Tangents[nPoints-1]= (YResp[nPoints-1] - YResp[nPoints-2]) / (XPos[nPoints-1] - XPos[nPoints-2]);
622 
623  // Compute cubic coefficients for each segment
624  for (int i = 0; i < nPoints - 1; ++i)
625  {
626  M3::float_t dx = XPos[i+1] - XPos[i];
627  M3::float_t dy = YResp[i+1] - YResp[i];
628  M3::float_t m0 = Tangents[i];
629  M3::float_t m1 = Tangents[i+1];
630 
631  Par[i][0] = m0; // b
632  Par[i][1] = (3*dy/(dx*dx)) - (2*m0 + m1)/dx; // c
633  Par[i][2] = (m0 + m1 - 2*dy/dx) / (dx*dx); // d
634 
635  MACH3LOG_TRACE("KB segment {}: dx={}, dy={}, b={}, c={}, d={}",
636  i, dx, dy, Par[i][0], Par[i][1], Par[i][2]);
637  }
638  // Last "virtual" segment
639  Par[nPoints-1][0] = 0.0;
640  Par[nPoints-1][1] = 0.0;
641  Par[nPoints-1][2] = 0.0;
642 
643  delete[] Tangents;
644  }
645  else
646  {
647  MACH3LOG_ERROR("Unsupported interpolation type {}", static_cast<int>(InterPolation));
648  throw MaCh3Exception(__FILE__ , __LINE__ );
649  }
650 
651  delete spline;
652  spline = nullptr;
653  }
#define MACH3LOG_TRACE
Definition: MaCh3Logger.h:33
@ kTSpline3
Default TSpline3 interpolation.
@ kMonotonic
EM: DOES NOT make the entire spline monotonic, only the segments.
@ kKochanekBartels
KS: Kochanek-Bartels spline: allows local control of tension, continuity, and bias.
@ kLinear
Linear interpolation between knots.
@ kLinearFunc
Liner interpolation using TF1 not spline.
@ kAkima
EM: Akima spline iis allowed to be discontinuous in 2nd derivative and coefficients in any segment.
Custom exception class used throughout MaCh3.
int int_t
Definition: Core.h:38

Member Data Documentation

◆ nPoints

M3::int_t TSpline3_red::nPoints
protected

Number of points/knot in TSpline3.

Definition at line 770 of file SplineStructs.h.

◆ Par

M3::float_t** TSpline3_red::Par
protected

Always uses a third order polynomial, so hard-code the number of coefficients in implementation.

Definition at line 772 of file SplineStructs.h.

◆ XPos

M3::float_t* TSpline3_red::XPos
protected

Positions of each x for each knot.

Definition at line 774 of file SplineStructs.h.

◆ YResp

M3::float_t* TSpline3_red::YResp
protected

y-value for each knot

Definition at line 776 of file SplineStructs.h.


The documentation for this class was generated from the following file: