MaCh3  2.4.2
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 [26]. 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 265 of file SplineStructs.h.

Constructor & Destructor Documentation

◆ TSpline3_red() [1/3]

TSpline3_red::TSpline3_red ( )
inline

Empty constructor.

Definition at line 269 of file SplineStructs.h.

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

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

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

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

Member Function Documentation

◆ ConstructTSpline3()

TSpline3* TSpline3_red::ConstructTSpline3 ( )
inline

CW: Make a TSpline3 from the reduced splines.

Definition at line 737 of file SplineStructs.h.

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

◆ Eval()

double TSpline3_red::Eval ( double  var)
inlineoverridevirtual

CW: Evaluate the weight from a variation.

Implements TResponseFunction_red.

Definition at line 707 of file SplineStructs.h.

707  {
708  // Get the segment for this variation
709  int segment = FindX(var);
710  // The get the coefficients for this variation
711  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);
712  GetCoeff(segment, x, y, b, c, d);
713  double dx = var - x;
714  // Evaluate the third order polynomial
715  double weight = y+dx*(b+dx*(c+d*dx));
716  return weight;
717  }
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 673 of file SplineStructs.h.

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

◆ 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 728 of file SplineStructs.h.

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

◆ GetKnot()

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

Definition at line 722 of file SplineStructs.h.

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

◆ GetNp()

M3::int_t TSpline3_red::GetNp ( )
inlineoverridevirtual

CW: Get the number of points.

Implements TResponseFunction_red.

Definition at line 720 of file SplineStructs.h.

720 { return nPoints; }

◆ Print()

void TSpline3_red::Print ( )
inlineoverridevirtual

Print detailed info.

Implements TResponseFunction_red.

Definition at line 758 of file SplineStructs.h.

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

◆ SetFunc()

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

Set the function [26].

Definition at line 311 of file SplineStructs.h.

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

◆ XPos

M3::float_t* TSpline3_red::XPos
protected

Positions of each x for each knot.

Definition at line 773 of file SplineStructs.h.

◆ YResp

M3::float_t* TSpline3_red::YResp
protected

y-value for each knot

Definition at line 775 of file SplineStructs.h.


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