MaCh3  2.2.3
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 [25]. 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 267 of file SplineStructs.h.

Constructor & Destructor Documentation

◆ TSpline3_red() [1/3]

TSpline3_red::TSpline3_red ( )
inline

Empty constructor.

Definition at line 271 of file SplineStructs.h.

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

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

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

◆ ~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 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.

Note
See root/hist/hist/src/TSpline3::FindX(double) or FindSplineSegment

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:25

◆ SetFunc()

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

Set the function [25].

Definition at line 313 of file SplineStructs.h.

313  {
314  nPoints = M3::int_t(spline->GetNp());
315  if (Par != nullptr) {
316  for (int i = 0; i < nPoints; ++i) {
317  delete[] Par[i];
318  Par[i] = nullptr;
319  }
320  delete[] Par;
321  Par = nullptr;
322  }
323  if (XPos != nullptr) delete[] XPos;
324  if (YResp != nullptr) delete[] YResp;
325  // Save the parameters for each knot
326  Par = new M3::float_t*[nPoints];
327  // Save the positions of the knots
328  XPos = new M3::float_t[nPoints];
329  // Save the y response at each knot
330  YResp = new M3::float_t[nPoints];
331 
332  //KS: Default TSpline3 ROOT implementation
333  if(InterPolation == kTSpline3)
334  {
335  for (int i = 0; i < nPoints; ++i) {
336  // 3 is the size of the TSpline3 coefficients
337  Par[i] = new M3::float_t[3];
338  double x = -999.99, y = -999.99, b = -999.99, c = -999.99, d = -999.99;
339  spline->GetCoeff(i, x, y, b, c, d);
340  XPos[i] = M3::float_t(x);
341  YResp[i] = M3::float_t(y);
342  Par[i][0] = M3::float_t(b);
343  Par[i][1] = M3::float_t(c);
344  Par[i][2] = M3::float_t(d);
345  }
346  }
347  //CW: Reduce to use linear spline interpolation for certain parameters
348  // 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
349  // Get which parameters should be linear from the fit manager
350  // Convert the spline number to global xsec parameter
351  // Loop over the splines points
352  // KS: kLinearFunc should be used with TF1, this is just as safety
353  else if(InterPolation == kLinear || InterPolation == kLinearFunc)
354  {
355  for (int k = 0; k < nPoints; ++k) {
356  Par[k] = new M3::float_t[3];
357 
358  Double_t x1, y1, b1, c1, d1, x2, y2, b2, c2, d2 = 0;
359  spline->GetCoeff(k, x1, y1, b1, c1, d1);
360 
361  double tempb = 0;
362  if (k == nPoints - 1) {
363  tempb = Par[k-1][0];
364  } else {
365  spline->GetCoeff(k + 1, x2, y2, b2, c2, d2);
366  tempb = (y2-y1)/(x2-x1);
367  }
368  XPos[k] = M3::float_t(x1);
369  YResp[k] = M3::float_t(y1);
370  Par[k][0] = M3::float_t(tempb); // linear slope
371  Par[k][1] = M3::float_t(0);
372  Par[k][2] = M3::float_t(0);
373  }
374  }
375  //EM: Akima spline is similar to regular cubic spline but is allowed to be discontinuous in 2nd derivative and coefficients in any segment
376  // only depend on th 2 nearest points on either side
377  else if(InterPolation == kAkima)
378  {
379  // get the knot values for the spline
380  for (int i = 0; i < nPoints; ++i) {
381  // 3 is the size of the TSpline3 coefficients
382  Par[i] = new M3::float_t[3];
383 
384  double x = -999.99, y = -999.99;
385  spline->GetKnot(i, x, y);
386 
387  XPos[i] = M3::float_t(x);
388  YResp[i] = M3::float_t(y);
389  }
390 
391  M3::float_t* mvals = new M3::float_t[nPoints + 3];
392  M3::float_t* svals = new M3::float_t[nPoints + 1];
393 
394  for (int i = -2; i <= nPoints; ++i) {
395  // if segment is first or last or 2nd to first or last, needs to be dealt with slightly differently;
396  // need to estimate the values for additinal points which would lie outside of the spline
397  if(i ==-2){
398  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]));
399  }
400  else if(i==-1){
401  mvals[i+2] = M3::float_t(2.0 * (YResp[1] - YResp[0]) / (XPos[1] - XPos[0]) - (YResp[2] - YResp[1]) / (XPos[2] - XPos[1]));
402  }
403  else if(i==nPoints){
404  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]));
405  }
406  else if(i == nPoints - 1){
407  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]));
408  }
409  //standard internal segment
410  else{
411  mvals[i+2] = (YResp[i+1] - YResp[i])/ (XPos[i+1] - XPos[i]);
412  }
413  }
414 
415  for(int i =2; i<=nPoints+2; i++){
416  if (std::abs(mvals[i+1] - mvals[i]) + std::abs(mvals[i-1] - mvals[i-2]) != 0.0){
417  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]));
418  }
419  else{svals[i-2] = mvals[i];}
420  }
421 
422  // calculate the coefficients for the spline
423  for(int i = 0; i <nPoints; i++){
424  M3::float_t b, c, d = M3::float_t(-999.999);
425 
426  b = svals[i];
427  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]);
428  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]));
429 
430  Par[i][0] = b;
431  Par[i][1] = c;
432  Par[i][2] = d;
433  }
434 
435  // check the input spline for linear segments, if there are any then overwrite the calculated coefficients
436  // 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
437  for(int i = 0; i <nPoints-1; i++){
438  double x = -999.99, y = -999.99, b = -999.99, c = -999.99, d = -999.99;
439  spline->GetCoeff(i, x, y, b, c, d);
440 
441  if((c == 0.0 && d == 0.0)){
442  Par[i][0] = M3::float_t(b);
443  Par[i][1] = M3::float_t(0.0);
444  Par[i][2] = M3::float_t(0.0);
445  }
446  }
447  delete[] mvals;
448  delete[] svals;
449  }
450  //EM: Monotone spline is similar to regular cubic spline but enforce the condition that the interpolated value at any point
451  // must be between its two nearest knots, DOES NOT make the entire spline monotonic, only the segments
452  else if(InterPolation == kMonotonic)
453  {
454  // values of the secants at each point (for calculating monotone spline)
455  M3::float_t * Secants = new M3::float_t[nPoints -1];
456  // values of the tangens at each point (for calculating monotone spline)
457  M3::float_t * Tangents = new M3::float_t[nPoints];
458 
459  // get the knot values for the spline
460  for (int i = 0; i < nPoints; ++i) {
461  // 3 is the size of the TSpline3 coefficients
462  Par[i] = new M3::float_t[3];
463 
464  double x = -999.99, y = -999.99;
465  spline->GetKnot(i, x, y);
466 
467  XPos[i] = M3::float_t(x);
468  YResp[i] = M3::float_t(y);
469 
470  Tangents[i] = 0.0;
471  }
472 
473  // deal with the case of two points (just do linear interpolation between them)
474  if (nPoints == 2){
475  Par[0][0] = (YResp[1] - YResp[0]) / (XPos[1] - XPos[0]);
476  Par[0][1] = 0.0;
477  Par[0][2] = 0.0;
478  // extra "virtual" segment at end to make Par array shape fit with knot arrays shapes
479  Par[1][1] = 0.0;
480  Par[1][2] = 0.0;
481 
482  return;
483  } // if nPoints !=2 do full monotonic spline treatment:
484  else
485  {
486  // first pass over knots to calculate the secants
487  for (int i = 0; i < nPoints-1; ++i) {
488  Secants[i] = (YResp[i+1] - YResp[i]) / (XPos[i+1] - XPos[i]);
489  MACH3LOG_TRACE("Secant {}: {}", i, Secants[i]);
490  }
491 
492  Tangents[0] = Secants[0];
493  Tangents[nPoints-1] = Secants[nPoints -2];
494 
495  M3::float_t alpha;
496  M3::float_t beta;
497 
498  // second pass over knots to calculate tangents
499  for (int i = 1; i < nPoints-1; ++i) {
500  if ((Secants[i-1] >= 0.0 && Secants[i] >= 0.0) | (Secants[i-1] < 0.0 && Secants[i] < 0.0)){ //check for same sign
501  Tangents[i] = M3::float_t((Secants[i-1] + Secants[i]) /2.0);
502  }
503  }
504 
505  // third pass over knots to rescale tangents
506  for (int i = 0; i < nPoints-1; ++i) {
507  if (Secants[i] == 0.0){
508  Tangents[i] = 0.0;
509  Tangents[i+1] = 0.0;
510  }
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:23
@ 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 for MaCh3 errors.
int int_t
Definition: Core.h:31

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: