MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
Classes | Functions
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_redCreateResponseFunction (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.

51 {
52// ***************************************************************************
53 if(xsecgraph == nullptr){
54 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));
55 throw MaCh3Exception(__FILE__ , __LINE__ );
56 }
57
58 // EM: cap the weights of the knots if specified in the config
59 if(
60 (XsecCov->GetParSplineKnotUpperBound(splineParsIndex) != M3::DefSplineKnotUpBound)
61 || (XsecCov->GetParSplineKnotLowerBound(splineParsIndex) != M3::DefSplineKnotLowBound))
62 {
63 for(int knotId = 0; knotId < xsecgraph->GetN(); knotId++){
64 double x,y;
65
66 // EM: get the x and y of the point. Also double check that the requested point was valid just to be super safe
67 if( xsecgraph->GetPoint(knotId, x, y) == -1) {
68 MACH3LOG_ERROR("Invalid knot requested: {}", knotId);
69 throw MaCh3Exception(__FILE__ , __LINE__ );
70 }
71
72 y = std::min(y, XsecCov->GetParSplineKnotUpperBound(splineParsIndex));
73 y = std::max(y, XsecCov->GetParSplineKnotLowerBound(splineParsIndex));
74
75 xsecgraph->SetPoint(knotId, x, y);
76 }
77
78 // EM: check if our cap made the spline flat, if so we set to 1 knot to avoid problems later on
79 bool isFlat = true;
80 for(int knotId = 0; knotId < xsecgraph->GetN(); knotId++){
81 double x,y;
82
83 // EM: get the x and y of the point. Also double check that the requested point was valid just to be super safe
84 if( xsecgraph->GetPoint(knotId, x, y) == -1) {
85 MACH3LOG_ERROR("Invalid knot requested: {}", knotId);
86 throw MaCh3Exception(__FILE__ , __LINE__ );
87 }
88 if(std::abs(y - 1.0) > 1e-5) isFlat = false;
89 }
90
91 if(isFlat){
92 xsecgraph->Set(1);
93 xsecgraph->SetPoint(0, 0.0, 1.0);
94 }
95 }
96}
std::vector< bool > isFlat
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:25
Custom exception class for MaCh3 errors.
double GetParSplineKnotUpperBound(const int i) const
EM: value at which we cap spline knot weight.
std::string GetParFancyName(const int i) const
Get fancy name of the Parameter.
double GetParSplineKnotLowerBound(const int i) const
EM: value at which we cap spline knot weight.
static constexpr const double DefSplineKnotUpBound
Default value for spline knot capping, default mean not capping is being applied.
Definition: Core.h:76
static constexpr const double DefSplineKnotLowBound
Default value for spline knot capping, default mean not capping is being applied.
Definition: Core.h:78

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

100 {
101// ***************************************************************************
102 if(Spline == nullptr) {
103 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));
104 throw MaCh3Exception(__FILE__ , __LINE__ );
105 }
106
107 std::string oldName = Spline->GetName();
108 // EM: cap the weights of the knots if specified in the config
109 if((XsecCov->GetParSplineKnotUpperBound(splineParsIndex) != M3::DefSplineKnotUpBound) || (XsecCov->GetParSplineKnotLowerBound(splineParsIndex) != M3::DefSplineKnotLowBound))
110 {
111 const int NValues = Spline->GetNp();
112 std::vector<double> XVals(NValues);
113 std::vector<double> YVals(NValues);
114 for(int knotId = 0; knotId < NValues; knotId++){
115 double x,y;
116 // Extract X and Y
117 Spline->GetKnot(knotId, x, y);
118
119 y = std::min(y, XsecCov->GetParSplineKnotUpperBound(splineParsIndex));
120 y = std::max(y, XsecCov->GetParSplineKnotLowerBound(splineParsIndex));
121
122 XVals[knotId] = x;
123 YVals[knotId] = y;
124 }
125 delete Spline;
126 // If we capped we have to make new TSpline3 to recalculate coefficients
127 Spline = new TSpline3(oldName.c_str(), XVals.data(), YVals.data(), NValues);
128 }
129}

◆ CreateResponseFunction()

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

KS: Create Response Function using TGraph.

Parameters
graphThis holds weights for each "spline" knot
SplineRespFuncTypeType of response function, whether we use Spline or TF1
SplineInterpolationTypeInterpolation type for example Monotonic Akima etc.
Titletitle 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.

802 {
803// *********************************
804 TResponseFunction_red* RespFunc = nullptr;
805
806 if (graph && graph->GetN() > 1)
807 {
808 if(SplineRespFuncType == kTSpline3_red)
809 {
810 // Here's the TSpline3
811 TSpline3* spline = nullptr;
812 TSpline3_red *spline_red = nullptr;
813
814 // Create the TSpline3* from the TGraph* and build the coefficients
815 spline = new TSpline3(Title.c_str(), graph);
816 spline->SetNameTitle(Title.c_str(), Title.c_str());
817
818 // Make the reduced TSpline3 format and delete the old spline
819 spline_red = new TSpline3_red(spline, SplineInterpolationType);
820
821 RespFunc = spline_red;
822 }
823 else if(SplineRespFuncType == kTF1_red)
824 {
825 // The TF1 object we build from fitting the TGraph
826 TF1 *Fitter = nullptr;
827 TF1_red *tf1_red = nullptr;
828
829 if(graph->GetN() != 2) {
830 MACH3LOG_ERROR("Trying to make TF1 from more than 2 knots. Knots = {}", graph->GetN());
831 MACH3LOG_ERROR("Currently support only linear with 2 knots :(");
832 throw MaCh3Exception(__FILE__ , __LINE__ );
833 }
834
835 // Try simple linear function
836 Fitter = new TF1(Title.c_str(), "([1]+[0]*x)", graph->GetX()[0], graph->GetX()[graph->GetN()-1]);
837 //CW: For 2p2h shape C and O we can't fit a polynomial: try a linear combination of two linear functions around 0
838 //Fitter = new TF1(Title.c_str(), "(x<=0)*(1+[0]*x)+(x>0)*([1]*x+1)", graph->GetX()[0], graph->GetX()[graph->GetN()-1]);
839 // Fit 5hd order polynomial for all other parameters
840 //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]);
841 //Pseudo Heaviside for Pauli Blocking
842 //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]);
843
844 // Fit the TF1 to the graph
845 graph->Fit(Fitter, "Q0");
846 // Make the reduced TF1 if we want
847 tf1_red = new TF1_red(Fitter);
848
849 RespFunc = tf1_red;
850 }
851 else
852 {
853 MACH3LOG_ERROR("Unsupported response function type");
854 throw MaCh3Exception(__FILE__ , __LINE__ );
855 }
856 }
857 else
858 {
859 RespFunc = nullptr;
860 }
861 return RespFunc;
862}
@ kTF1_red
Uses TF1_red for interpolation.
@ kTSpline3_red
Uses TSpline3_red for interpolation.
CW: A reduced TF1 class only. Only saves parameters for each TF1 and how many parameters each paramet...
KS: A reduced ResponseFunction Generic function used for evaluating weight.
CW: Reduced TSpline3 class.

◆ isFlat()

bool isFlat ( TSpline3_red *&  spl)
inline

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

Parameters
splpointer to TSpline3_red that will be checked

Definition at line 701 of file SplineStructs.h.

701 {
702// *****************************************
703 int Np = spl->GetNp();
704 M3::float_t x, y, b, c, d;
705 // Go through spline segment parameters,
706 // Get y values for each spline knot,
707 // Every knot must evaluate to 1.0 to create a flat spline
708 for(int i = 0; i < Np; i++) {
709 spl->GetCoeff(i, x, y, b, c, d);
710 if (y != 1) {
711 return false;
712 }
713 }
714 return true;
715}
M3::int_t GetNp() override
CW: Get the number of points.
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.
double float_t
Definition: Core.h:28

◆ ReduceTF1()

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

CW: Reduced the TF1 to TF1_red.

Parameters
MasterSplineVector of TF1_red pointers which we strip back

Definition at line 758 of file SplineStructs.h.

758 {
759// *********************************
760 std::vector<std::vector<TF1*> >::iterator OuterIt;
761 std::vector<TF1*>::iterator InnerIt;
762
763 // The return vector
764 std::vector<std::vector<TF1_red*> > ReducedVector;
765 ReducedVector.reserve(MasterSpline.size());
766
767 // Loop over each parameter
768 int OuterCounter = 0;
769 for (OuterIt = MasterSpline.begin(); OuterIt != MasterSpline.end(); ++OuterIt, ++OuterCounter) {
770 // Make the temp vector
771 std::vector<TF1_red*> TempVector;
772 TempVector.reserve(OuterIt->size());
773 int InnerCounter = 0;
774 // Loop over each TSpline3 pointer
775 for (InnerIt = OuterIt->begin(); InnerIt != OuterIt->end(); ++InnerIt, ++InnerCounter) {
776 // Here's our delicious TSpline3 object
777 TF1* spline = (*InnerIt);
778 // Now make the reduced TSpline3 pointer (which deleted TSpline3)
779 TF1_red* red = NULL;
780 if (spline != NULL) {
781 red = new TF1_red(spline);
782 (*InnerIt) = spline;
783 }
784 // Push back onto new vector
785 TempVector.push_back(red);
786 } // End inner for loop
787 ReducedVector.push_back(TempVector);
788 } // End outer for loop
789 // Now have the reduced vector
790 return ReducedVector;
791}

◆ ReduceTSpline3()

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

CW: Reduced the TSpline3 to TSpline3_red.

Parameters
MasterSplineVector of TSpline3_red pointers which we strip back

Definition at line 720 of file SplineStructs.h.

720 {
721// *********************************
722 std::vector<std::vector<TSpline3*> >::iterator OuterIt;
723 std::vector<TSpline3*>::iterator InnerIt;
724
725 // The return vector
726 std::vector<std::vector<TSpline3_red*> > ReducedVector;
727 ReducedVector.reserve(MasterSpline.size());
728
729 // Loop over each parameter
730 int OuterCounter = 0;
731 for (OuterIt = MasterSpline.begin(); OuterIt != MasterSpline.end(); ++OuterIt, ++OuterCounter) {
732 // Make the temp vector
733 std::vector<TSpline3_red*> TempVector;
734 TempVector.reserve(OuterIt->size());
735 int InnerCounter = 0;
736 // Loop over each TSpline3 pointer
737 for (InnerIt = OuterIt->begin(); InnerIt != OuterIt->end(); ++InnerIt, ++InnerCounter) {
738 // Here's our delicious TSpline3 object
739 TSpline3 *spline = (*InnerIt);
740 // Now make the reduced TSpline3 pointer
741 TSpline3_red *red = NULL;
742 if (spline != NULL) {
743 red = new TSpline3_red(spline);
744 (*InnerIt) = spline;
745 }
746 // Push back onto new vector
747 TempVector.push_back(red);
748 } // End inner for loop
749 ReducedVector.push_back(TempVector);
750 } // End outer for loop
751 // Now have the reduced vector
752 return ReducedVector;
753}