9#pragma GCC diagnostic push
10#pragma GCC diagnostic ignored "-Wuseless-cast"
11#pragma GCC diagnostic ignored "-Wfloat-conversion"
12#pragma GCC diagnostic ignored "-Wnull-dereference"
39 std::vector<M3::float_t>
xPts;
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));
63 for(
int knotId = 0; knotId < xsecgraph->GetN(); knotId++){
67 if( xsecgraph->GetPoint(knotId, x, y) == -1) {
75 xsecgraph->SetPoint(knotId, x, y);
80 for(
int knotId = 0; knotId < xsecgraph->GetN(); knotId++){
84 if( xsecgraph->GetPoint(knotId, x, y) == -1) {
88 if(std::abs(y - 1.0) > 1e-5)
isFlat =
false;
93 xsecgraph->SetPoint(0, 0.0, 1.0);
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));
107 std::string oldName = Spline->GetName();
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++){
117 Spline->GetKnot(knotId, x, y);
127 Spline =
new TSpline3(oldName.c_str(), XVals.data(), YVals.data(), NValues);
141 virtual double Eval(
const double var)=0;
170 for (
int i = 0; i <
length; ++i) {
184 if (
Par != NULL)
delete[]
Par;
186 for (
int i = 0; i <
length; ++i) {
194 inline double Eval(
const double var)
override {
218 Par[Parameter] = Value;
224 MACH3LOG_ERROR(
"You requested parameter number {} but length is {} parameters", Parameter,
length);
228 return Par[Parameter];
242 for (
int i = 0; i <
length; i++) {
248 inline TF1*
ConstructTF1(
const std::string& function,
const int xmin,
const int xmax) {
249 TF1 *func =
new TF1(
"TF1", function.c_str(), xmin, xmax);
250 for(
int i = 0; i <
length; ++i) {
251 func->SetParameter(i,
Par[i]);
283 SetFunc(spline, InterPolation);
295 for(
int j = 0; j < N; ++j){
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(
"*********************************************************************************************");
315 for (
int i = 0; i <
nPoints; ++i) {
334 for (
int i = 0; i <
nPoints; ++i) {
337 double x = -999.99, y = -999.99, b = -999.99, c = -999.99, d = -999.99;
338 spline->GetCoeff(i, x, y, b, c, d);
354 for (
int k = 0; k <
nPoints; ++k) {
357 Double_t x1, y1, b1, c1, d1, x2, y2, b2, c2, d2 = 0;
358 spline->GetCoeff(k, x1, y1, b1, c1, d1);
359 spline->GetCoeff(k+1, x2, y2, b2, c2, d2);
360 double tempb = (y2-y1)/(x2-x1);
371 else if(InterPolation ==
kAkima)
374 for (
int i = 0; i <
nPoints; ++i) {
378 double x = -999.99, y = -999.99;
379 spline->GetKnot(i, x, y);
388 for (
int i = -2; i <=
nPoints; ++i) {
409 for(
int i =2; i<=
nPoints+2; i++){
410 if (std::abs(mvals[i+1] - mvals[i]) + std::abs(mvals[i-1] - mvals[i-2]) != 0.0){
411 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]));
413 else{svals[i-2] = mvals[i];}
417 for(
int i = 0; i <
nPoints; i++){
431 for(
int i = 0; i <
nPoints-1; i++){
432 double x = -999.99, y = -999.99, b = -999.99, c = -999.99, d = -999.99;
433 spline->GetCoeff(i, x, y, b, c, d);
435 if((c == 0.0 && d == 0.0)){
454 for (
int i = 0; i <
nPoints; ++i) {
458 double x = -999.99, y = -999.99;
459 spline->GetKnot(i, x, y);
481 for (
int i = 0; i <
nPoints-1; ++i) {
486 Tangents[0] = Secants[0];
493 for (
int i = 1; i <
nPoints-1; ++i) {
494 if ((Secants[i-1] >= 0.0 && Secants[i] >= 0.0) | (Secants[i-1] < 0.0 && Secants[i] < 0.0)){
495 Tangents[i] =
M3::float_t((Secants[i-1] + Secants[i]) /2.0);
500 for (
int i = 0; i <
nPoints-1; ++i) {
501 if (Secants[i] == 0.0){
507 alpha = Tangents[i] / Secants[i];
508 beta = Tangents[i+1] / Secants[i];
517 if (alpha * alpha + beta * beta >9.0){
519 Tangents[i] = tau * alpha * Secants[i];
520 Tangents[i+1] = tau * beta * Secants[i];
526 for(
int i = 0; i <
nPoints-1; i++){
530 b = Tangents[i] * dx;
535 Par[i][1] = c / (dx * dx);
536 Par[i][2] = d / (dx * dx * dx);
538 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)){
539 MACH3LOG_INFO(
"Bad spline parameters for segment {}: (b, c, d) = {}, {}, {}. This will cause problems with GPU.",
555 for(
int i = 0; i <
nPoints-1; i++){
556 double x = -999.99, y = -999.99, b = -999.99, c = -999.99, d = -999.99;
557 spline->GetCoeff(i, x, y, b, c, d);
559 if((c == 0.0 && d == 0.0)){
571 MACH3LOG_ERROR(
"Unsupported interpolation type {}",
static_cast<int>(InterPolation));
582 for (
int i = 0; i <
nPoints; ++i) {
583 if (
Par[i] != NULL) {
614 while (kHigh - segment > 1) {
616 kHalf = (segment + kHigh)/2;
618 if (x >
XPos[kHalf]) {
631 inline double Eval(
double var)
override {
633 int segment =
FindX(var);
639 double weight = y+dx*(b+dx*(c+d*dx));
663 #ifdef _LOW_MEMORY_STRUCTS_
664 std::vector<Double_t> xPosDoubles(
nPoints);
665 std::vector<Double_t> yPosDoubles(
nPoints);
666 for (Int_t i = 0; i <
nPoints; ++i) {
667 xPosDoubles[i] =
static_cast<Double_t
>(
XPos[i]);
668 yPosDoubles[i] =
static_cast<Double_t
>(
YResp[i]);
670 TSpline3 *spline =
new TSpline3(
"Spline", xPosDoubles.data(), yPosDoubles.data(),
static_cast<int>(
nPoints));
681 for (
int i = 0; i <
nPoints; ++i) {
703 int Np = spl->
GetNp();
708 for(
int i = 0; i < Np; i++) {
720inline std::vector<std::vector<TSpline3_red*> >
ReduceTSpline3(std::vector<std::vector<TSpline3*> > &MasterSpline) {
722 std::vector<std::vector<TSpline3*> >::iterator OuterIt;
723 std::vector<TSpline3*>::iterator InnerIt;
726 std::vector<std::vector<TSpline3_red*> > ReducedVector;
727 ReducedVector.reserve(MasterSpline.size());
730 int OuterCounter = 0;
731 for (OuterIt = MasterSpline.begin(); OuterIt != MasterSpline.end(); ++OuterIt, ++OuterCounter) {
733 std::vector<TSpline3_red*> TempVector;
734 TempVector.reserve(OuterIt->size());
735 int InnerCounter = 0;
737 for (InnerIt = OuterIt->begin(); InnerIt != OuterIt->end(); ++InnerIt, ++InnerCounter) {
739 TSpline3 *spline = (*InnerIt);
742 if (spline != NULL) {
747 TempVector.push_back(red);
749 ReducedVector.push_back(TempVector);
752 return ReducedVector;
758inline std::vector<std::vector<TF1_red*> >
ReduceTF1(std::vector<std::vector<TF1*> > &MasterSpline) {
760 std::vector<std::vector<TF1*> >::iterator OuterIt;
761 std::vector<TF1*>::iterator InnerIt;
764 std::vector<std::vector<TF1_red*> > ReducedVector;
765 ReducedVector.reserve(MasterSpline.size());
768 int OuterCounter = 0;
769 for (OuterIt = MasterSpline.begin(); OuterIt != MasterSpline.end(); ++OuterIt, ++OuterCounter) {
771 std::vector<TF1_red*> TempVector;
772 TempVector.reserve(OuterIt->size());
773 int InnerCounter = 0;
775 for (InnerIt = OuterIt->begin(); InnerIt != OuterIt->end(); ++InnerIt, ++InnerCounter) {
777 TF1* spline = (*InnerIt);
780 if (spline != NULL) {
785 TempVector.push_back(red);
787 ReducedVector.push_back(TempVector);
790 return ReducedVector;
802 const std::string& Title) {
806 if (graph && graph->GetN() > 1)
811 TSpline3* spline =
nullptr;
815 spline =
new TSpline3(Title.c_str(), graph);
816 spline->SetNameTitle(Title.c_str(), Title.c_str());
819 spline_red =
new TSpline3_red(spline, SplineInterpolationType);
821 RespFunc = spline_red;
823 else if(SplineRespFuncType ==
kTF1_red)
826 TF1 *Fitter =
nullptr;
829 if(graph->GetN() != 2) {
830 MACH3LOG_ERROR(
"Trying to make TF1 from more than 2 knots. Knots = {}", graph->GetN());
836 Fitter =
new TF1(Title.c_str(),
"([1]+[0]*x)", graph->GetX()[0], graph->GetX()[graph->GetN()-1]);
845 graph->Fit(Fitter,
"Q0");
864#pragma GCC diagnostic pop
std::vector< bool > isFlat
SplineInterpolation
Make an enum of the spline interpolation type.
@ kTSpline3
Default TSpline3 interpolation.
@ kMonotonic
EM: DOES NOT make the entire spline monotonic, only the segments.
@ 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.
RespFuncType
Make an enum of the spline interpolation type.
@ kTF1_red
Uses TF1_red for interpolation.
@ kTSpline3_red
Uses TSpline3_red for interpolation.
std::vector< std::vector< TSpline3_red * > > ReduceTSpline3(std::vector< std::vector< TSpline3 * > > &MasterSpline)
CW: Reduced the TSpline3 to TSpline3_red.
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 i...
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 i...
TResponseFunction_red * CreateResponseFunction(TGraph *&graph, const RespFuncType SplineRespFuncType, const SplineInterpolation SplineInterpolationType, const std::string &Title)
KS: Create Response Function using TGraph.
std::vector< std::vector< TF1_red * > > ReduceTF1(std::vector< std::vector< TF1 * > > &MasterSpline)
CW: Reduced the TF1 to TF1_red.
Custom exception class for MaCh3 errors.
Class responsible for handling of systematic error parameters with different types defined in the con...
CW: A reduced TF1 class only. Only saves parameters for each TF1 and how many parameters each paramet...
TF1 * ConstructTF1(const std::string &function, const int xmin, const int xmax)
KS: Make a TF1 from the reduced TF1.
void SetSize(M3::int_t nSpline)
Set the size.
void SetFunc(TF1 *&Func)
Set the function.
TF1_red()
Empty constructor.
void Print() override
Print detailed info.
M3::int_t GetNp() override
DL: Get number of points.
TF1_red(TF1 *&Function)
The TF1 constructor with deep copy.
M3::float_t * Par
The parameters.
virtual ~TF1_red()
Empty destructor.
double Eval(const double var) override
Evaluate a variation.
TF1_red(M3::int_t nSize, M3::float_t *Array)
The useful constructor with deep copy.
double GetParameter(M3::int_t Parameter)
Get a parameter value.
void SetParameter(M3::int_t Parameter, M3::float_t Value)
Set a parameter to a value.
int GetSize()
Get the size.
KS: A reduced ResponseFunction Generic function used for evaluating weight.
virtual ~TResponseFunction_red()
Empty destructor.
TResponseFunction_red()
Empty constructor.
virtual M3::int_t GetNp()=0
DL: Get number of points.
virtual double Eval(const double var)=0
Evaluate a variation.
virtual void Print()=0
KS: Printer.
CW: Reduced TSpline3 class.
double Eval(double var) override
CW: Evaluate the weight from a variation.
void GetKnot(int i, M3::float_t &xtmp, M3::float_t &ytmp)
int FindX(double x)
Find the segment relevant to this variation in x See root/hist/hist/src/TSpline3FindX(double) or samp...
virtual ~TSpline3_red()
Empty destructor.
M3::float_t ** Par
Always uses a third order polynomial, so hard-code the number of coefficients in implementation.
TSpline3_red(M3::float_t *X, M3::float_t *Y, M3::int_t N, M3::float_t **P)
constructor taking parameters
TSpline3 * ConstructTSpline3()
CW: Make a TSpline3 from the reduced splines.
M3::int_t nPoints
Number of points/knot in TSpline3.
TSpline3_red(TSpline3 *&spline, SplineInterpolation InterPolation=kTSpline3)
The constructor that takes a TSpline3 pointer and copies in to memory.
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.
void Print() override
Print detailed info.
void SetFunc(TSpline3 *&spline, SplineInterpolation InterPolation=kTSpline3)
Set the function.
M3::float_t * XPos
Positions of each x for each knot.
TSpline3_red()
Empty constructor.
M3::float_t * YResp
y-value for each knot
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.
static constexpr const double DefSplineKnotLowBound
Default value for spline knot capping, default mean not capping is being applied.
CW: Add a struct to hold info about the splinified xsec parameters and help with FindSplineSegment.
virtual ~FastSplineInfo()
Destructor.
std::vector< M3::float_t > xPts
Array of the knots positions.
const double * splineParsPointer
Array of the knots positions.
M3::int_t CurrSegment
Array of what segment of spline we're currently interested in. Gets updated once per MCMC iteration.
FastSplineInfo()
Constructor.
M3::int_t nPts
Number of points in spline.