SplineStructs.h File Reference

Contains structures and helper functions for handling spline representations of systematic parameters in the MaCh3. More...

#include "Parameters/ParameterHandlerGeneric.h"
#include "Parameters/ParameterStructs.h"
#include <cmath>

Include dependency graph for SplineStructs.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	FastSplineInfo
	CW: Add a struct to hold info about the splinified xsec parameters and help with FindSplineSegment. More...
class	TResponseFunction_red
	KS: A reduced ResponseFunction Generic function used for evaluating weight. More...
class	TF1_red
	CW: A reduced TF1 class only. Only saves parameters for each TF1 and how many parameters each parameter set has. More...
class	TSpline3_red
	CW: Reduced TSpline3 class. More...

Functions
void	ApplyKnotWeightCap (TGraph *xsecgraph, const int splineParsIndex, ParameterHandlerGeneric *XsecCov)
	EM: Apply capping to knot weight for specified spline parameter. param graph needs to have been set in xsecgraph array first.
void	ApplyKnotWeightCapTSpline3 (TSpline3 *&Spline, const int splineParsIndex, ParameterHandlerGeneric *XsecCov)
	EM: Apply capping to knot weight for specified spline parameter. param graph needs to have been set in Spline array first.
bool	isFlat (TSpline3_red *&spl)
	CW: Helper function used in the constructor, tests to see if the spline is flat.
std::vector< std::vector< TSpline3_red * > >	ReduceTSpline3 (std::vector< std::vector< TSpline3 * > > &MasterSpline)
	CW: Reduced the TSpline3 to TSpline3_red.
std::vector< std::vector< TF1_red * > >	ReduceTF1 (std::vector< std::vector< TF1 * > > &MasterSpline)
	CW: Reduced the TF1 to TF1_red.
TResponseFunction_red *	CreateResponseFunction (TGraph *&graph, const RespFuncType SplineRespFuncType, const SplineInterpolation SplineInterpolationType, const std::string &Title)
	KS: Create Response Function using TGraph.

Detailed Description

Contains structures and helper functions for handling spline representations of systematic parameters in the MaCh3.

Author

Clarence Wret

Kamil Skwarczynski

Definition in file SplineStructs.h.

Function Documentation

ApplyKnotWeightCap()

void ApplyKnotWeightCap ( TGraph *  xsecgraph, const int  splineParsIndex, ParameterHandlerGeneric *  XsecCov  )

inline

EM: Apply capping to knot weight for specified spline parameter. param graph needs to have been set in xsecgraph array first.

Definition at line 51 of file SplineStructs.h.

                                                                                                               {
// ***************************************************************************
  if(xsecgraph == nullptr){
    MACH3LOG_ERROR("hmmm looks like you're trying to apply capping for spline parameter {} but it hasn't been set in xsecgraph yet", XsecCov->GetParFancyName(splineParsIndex));
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  // EM: cap the weights of the knots if specified in the config
  if(
    (XsecCov->GetParSplineKnotUpperBound(splineParsIndex) != M3::DefSplineKnotUpBound)
    || (XsecCov->GetParSplineKnotLowerBound(splineParsIndex) != M3::DefSplineKnotLowBound))
  {
    for(int knotId = 0; knotId < xsecgraph->GetN(); knotId++){
      double x,y;
 
      // EM: get the x and y of the point. Also double check that the requested point was valid just to be super safe
      if( xsecgraph->GetPoint(knotId, x, y) == -1) {
        MACH3LOG_ERROR("Invalid knot requested: {}", knotId);
        throw MaCh3Exception(__FILE__ , __LINE__ );
      }
 
      y = std::min(y, XsecCov->GetParSplineKnotUpperBound(splineParsIndex));
      y = std::max(y, XsecCov->GetParSplineKnotLowerBound(splineParsIndex));
 
      xsecgraph->SetPoint(knotId, x, y);
    }
 
    // EM: check if our cap made the spline flat, if so we set to 1 knot to avoid problems later on
    bool isFlat = true;
    for(int knotId = 0; knotId < xsecgraph->GetN(); knotId++){
      double x,y;
 
      // EM: get the x and y of the point. Also double check that the requested point was valid just to be super safe
      if( xsecgraph->GetPoint(knotId, x, y) == -1) {
        MACH3LOG_ERROR("Invalid knot requested: {}", knotId);
        throw MaCh3Exception(__FILE__ , __LINE__ );
      }
      if(std::abs(y - 1.0) > 1e-5) isFlat = false;
    }
 
    if(isFlat){
      xsecgraph->Set(1);
      xsecgraph->SetPoint(0, 0.0, 1.0);
    }
  }
}

ApplyKnotWeightCapTSpline3()

void ApplyKnotWeightCapTSpline3 ( TSpline3 *&  Spline, const int  splineParsIndex, ParameterHandlerGeneric *  XsecCov  )

inline

EM: Apply capping to knot weight for specified spline parameter. param graph needs to have been set in Spline array first.

Definition at line 100 of file SplineStructs.h.

                                                                                                                       {
// ***************************************************************************
  if(Spline == nullptr) {
    MACH3LOG_ERROR("hmmm looks like you're trying to apply capping for spline parameter {} but it hasn't been set in Spline yet", XsecCov->GetParFancyName(splineParsIndex));
    throw MaCh3Exception(__FILE__ , __LINE__ );
  }
 
  std::string oldName = Spline->GetName();
  // EM: cap the weights of the knots if specified in the config
  if((XsecCov->GetParSplineKnotUpperBound(splineParsIndex) != M3::DefSplineKnotUpBound) || (XsecCov->GetParSplineKnotLowerBound(splineParsIndex) != M3::DefSplineKnotLowBound))
  {
    const int NValues = Spline->GetNp();
    std::vector<double> XVals(NValues);
    std::vector<double> YVals(NValues);
    for(int knotId = 0; knotId < NValues; knotId++){
      double x,y;
      // Extract X and Y
      Spline->GetKnot(knotId, x, y);
 
      y = std::min(y, XsecCov->GetParSplineKnotUpperBound(splineParsIndex));
      y = std::max(y, XsecCov->GetParSplineKnotLowerBound(splineParsIndex));
 
      XVals[knotId] = x;
      YVals[knotId] = y;
    }
    delete Spline;
    // If we capped we have to make new TSpline3 to recalculate coefficients
    Spline = new TSpline3(oldName.c_str(), XVals.data(), YVals.data(), NValues);
  }
}

CreateResponseFunction()

TResponseFunction_red * CreateResponseFunction ( TGraph *&  graph, const RespFuncType  SplineRespFuncType, const SplineInterpolation  SplineInterpolationType, const std::string &  Title  )

inline

KS: Create Response Function using TGraph.

Parameters

graph	This holds weights for each "spline" knot
SplineRespFuncType	Type of response function, whether we use Spline or TF1
SplineInterpolationType	Interpolation type for example Monotonic Akima etc.
Title	title you want for ROOT object, isn't very useful as in MaCh3 we usually convert later to something light weight to not waste RAM.

Definition at line 799 of file SplineStructs.h.

                                                                             {
// *********************************
  TResponseFunction_red* RespFunc = nullptr;
 
  if (graph && graph->GetN() > 1)
  {
    if(SplineRespFuncType == kTSpline3_red)
    {
      // Here's the TSpline3
      TSpline3* spline = nullptr;
      TSpline3_red *spline_red = nullptr;
 
      // Create the TSpline3* from the TGraph* and build the coefficients
      spline = new TSpline3(Title.c_str(), graph);
      spline->SetNameTitle(Title.c_str(), Title.c_str());
 
      // Make the reduced TSpline3 format and delete the old spline
      spline_red = new TSpline3_red(spline, SplineInterpolationType);
 
      RespFunc = spline_red;
    }
    else if(SplineRespFuncType == kTF1_red)
    {
      // The TF1 object we build from fitting the TGraph
      TF1 *Fitter = nullptr;
      TF1_red *tf1_red = nullptr;
 
      if(graph->GetN() != 2) {
        MACH3LOG_ERROR("Trying to make TF1 from more than 2 knots.  Knots = {}", graph->GetN());
        MACH3LOG_ERROR("Currently support only linear with 2 knots :(");
        throw MaCh3Exception(__FILE__ , __LINE__ );
      }
 
      // Try simple linear function
      Fitter = new TF1(Title.c_str(), "([1]+[0]*x)", graph->GetX()[0], graph->GetX()[graph->GetN()-1]);
      //CW: For 2p2h shape C and O we can't fit a polynomial: try a linear combination of two linear functions around 0
      //Fitter = new TF1(Title.c_str(), "(x<=0)*(1+[0]*x)+(x>0)*([1]*x+1)", graph->GetX()[0], graph->GetX()[graph->GetN()-1]);
      // Fit 5hd order polynomial for all other parameters
      //Fitter = new TF1(Title.c_str(), "1+[0]*x+[1]*x*x+[2]*x*x*x+[3]*x*x*x*x+[4]*x*x*x*x*x", graph->GetX()[0], graph->GetX()[graph->GetN()-1]);
      //Pseudo Heaviside for Pauli Blocking
      //Fitter = new TF1(Title.c_str(), "(x <= 0)*(1+[0]*x) + (1 >= x)*(x > 0)*(1+[1]*x) + (x > 1)*([3]+[2]*x)", graph->GetX()[0], graph->GetX()[graph->GetN()-1]);
 
      // Fit the TF1 to the graph
      graph->Fit(Fitter, "Q0");
      // Make the reduced TF1 if we want
      tf1_red = new TF1_red(Fitter);
 
      RespFunc = tf1_red;
    }
    else
    {
      MACH3LOG_ERROR("Unsupported response function type");
      throw MaCh3Exception(__FILE__ , __LINE__ );
    }
  }
  else
  {
    RespFunc = nullptr;
  }
  return RespFunc;
}

isFlat()

bool isFlat ( TSpline3_red *&  spl )

inline

CW: Helper function used in the constructor, tests to see if the spline is flat.

Parameters

spl	pointer to TSpline3_red that will be checked

Definition at line 701 of file SplineStructs.h.

                                       {
// *****************************************
  int Np = spl->GetNp();
  M3::float_t x, y, b, c, d;
  // Go through spline segment parameters,
  // Get y values for each spline knot,
  // Every knot must evaluate to 1.0 to create a flat spline
  for(int i = 0; i < Np; i++) {
    spl->GetCoeff(i, x, y, b, c, d);
    if (y != 1) {
      return false;
    }
  }
  return true;
}

ReduceTF1()

std::vector< std::vector< TF1_red * > > ReduceTF1 ( std::vector< std::vector< TF1 * > > &  MasterSpline )

inline

CW: Reduced the TF1 to TF1_red.

Parameters

MasterSpline

Vector of TF1_red pointers which we strip back

Definition at line 758 of file SplineStructs.h.

                                                                                            {
// *********************************
  std::vector<std::vector<TF1*> >::iterator OuterIt;
  std::vector<TF1*>::iterator InnerIt;
 
  // The return vector
  std::vector<std::vector<TF1_red*> > ReducedVector;
  ReducedVector.reserve(MasterSpline.size());
 
  // Loop over each parameter
  int OuterCounter = 0;
  for (OuterIt = MasterSpline.begin(); OuterIt != MasterSpline.end(); ++OuterIt, ++OuterCounter) {
    // Make the temp vector
    std::vector<TF1_red*> TempVector;
    TempVector.reserve(OuterIt->size());
    int InnerCounter = 0;
    // Loop over each TSpline3 pointer
    for (InnerIt = OuterIt->begin(); InnerIt != OuterIt->end(); ++InnerIt, ++InnerCounter) {
      // Here's our delicious TSpline3 object
      TF1* spline = (*InnerIt);
      // Now make the reduced TSpline3 pointer (which deleted TSpline3)
      TF1_red* red = NULL;
      if (spline != NULL) {
        red = new TF1_red(spline);
        (*InnerIt) = spline;
      }
      // Push back onto new vector
      TempVector.push_back(red);
    } // End inner for loop
    ReducedVector.push_back(TempVector);
  } // End outer for loop
  // Now have the reduced vector
  return ReducedVector;
}

ReduceTSpline3()

std::vector< std::vector< TSpline3_red * > > ReduceTSpline3 ( std::vector< std::vector< TSpline3 * > > &  MasterSpline )

inline

CW: Reduced the TSpline3 to TSpline3_red.

Parameters

MasterSpline

Vector of TSpline3_red pointers which we strip back

Definition at line 720 of file SplineStructs.h.

                                                                                                           {
// *********************************
  std::vector<std::vector<TSpline3*> >::iterator OuterIt;
  std::vector<TSpline3*>::iterator InnerIt;
 
  // The return vector
  std::vector<std::vector<TSpline3_red*> > ReducedVector;
  ReducedVector.reserve(MasterSpline.size());
 
  // Loop over each parameter
  int OuterCounter = 0;
  for (OuterIt = MasterSpline.begin(); OuterIt != MasterSpline.end(); ++OuterIt, ++OuterCounter) {
    // Make the temp vector
    std::vector<TSpline3_red*> TempVector;
    TempVector.reserve(OuterIt->size());
    int InnerCounter = 0;
    // Loop over each TSpline3 pointer
    for (InnerIt = OuterIt->begin(); InnerIt != OuterIt->end(); ++InnerIt, ++InnerCounter) {
      // Here's our delicious TSpline3 object
      TSpline3 *spline = (*InnerIt);
      // Now make the reduced TSpline3 pointer
      TSpline3_red *red = NULL;
      if (spline != NULL) {
        red = new TSpline3_red(spline);
        (*InnerIt) = spline;
      }
      // Push back onto new vector
      TempVector.push_back(red);
    } // End inner for loop
    ReducedVector.push_back(TempVector);
  } // End outer for loop
  // Now have the reduced vector
  return ReducedVector;
}