MaCh3  2.6.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 [29]. More...
 
virtual ~TSpline3_red ()
 Empty destructor. More...
 
int FindX (double x) const
 Find the segment relevant to this variation in x. More...
 
double Eval (const double var) const override
 CW: Evaluate the weight from a variation. More...
 
M3::int_t GetNp () const override
 CW: Get the number of points. More...
 
void GetKnot (int i, M3::float_t &xtmp, M3::float_t &ytmp) const
 
void GetCoeff (int segment, M3::float_t &x, M3::float_t &y, M3::float_t &b, M3::float_t &c, M3::float_t &d) const
 CW: Get the coefficient of a given segment. More...
 
TSpline3 * ConstructTSpline3 ()
 CW: Make a TSpline3 from the reduced splines. More...
 
void Print () const 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 268 of file SplineStructs.h.

Constructor & Destructor Documentation

◆ TSpline3_red() [1/3]

TSpline3_red::TSpline3_red ( )
inline

Empty constructor.

Definition at line 272 of file SplineStructs.h.

273  nPoints = 0;
274  Par = nullptr;
275  XPos = nullptr;
276  YResp = nullptr;
277  }
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 280 of file SplineStructs.h.

281  Par = nullptr;
282  XPos = nullptr;
283  YResp = nullptr;
284  SetFunc(spline, InterPolation);
285  }
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 288 of file SplineStructs.h.

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

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

Member Function Documentation

◆ ConstructTSpline3()

TSpline3* TSpline3_red::ConstructTSpline3 ( )
inline

CW: Make a TSpline3 from the reduced splines.

Definition at line 740 of file SplineStructs.h.

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

◆ Eval()

double TSpline3_red::Eval ( const double  var) const
inlineoverridevirtual

CW: Evaluate the weight from a variation.

Implements TResponseFunction_red.

Definition at line 709 of file SplineStructs.h.

709  {
710  // Get the segment for this variation
711  int segment = FindX(var);
712  // The get the coefficients for this variation
713  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);
714  GetCoeff(segment, x, y, b, c, d);
715  double dx = var - x;
716  // Evaluate the third order polynomial
717  double weight = y+dx*(b+dx*(c+d*dx));
718  return weight;
719  }
int FindX(double x) const
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) const
CW: Get the coefficient of a given segment.

◆ FindX()

int TSpline3_red::FindX ( double  x) const
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 675 of file SplineStructs.h.

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

◆ 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 
) const
inline

CW: Get the coefficient of a given segment.

Definition at line 730 of file SplineStructs.h.

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

◆ GetKnot()

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

Definition at line 724 of file SplineStructs.h.

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

◆ GetNp()

M3::int_t TSpline3_red::GetNp ( ) const
inlineoverridevirtual

CW: Get the number of points.

Implements TResponseFunction_red.

Definition at line 722 of file SplineStructs.h.

722 { return nPoints; }

◆ Print()

void TSpline3_red::Print ( ) const
inlineoverridevirtual

Print detailed info.

Implements TResponseFunction_red.

Definition at line 761 of file SplineStructs.h.

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

◆ SetFunc()

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

Set the function [29].

Definition at line 314 of file SplineStructs.h.

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

◆ XPos

M3::float_t* TSpline3_red::XPos
protected

Positions of each x for each knot.

Definition at line 776 of file SplineStructs.h.

◆ YResp

M3::float_t* TSpline3_red::YResp
protected

y-value for each knot

Definition at line 778 of file SplineStructs.h.


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