MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
SplineBase.cpp
Go to the documentation of this file.
2
3#pragma GCC diagnostic ignored "-Wuseless-cast"
4#pragma GCC diagnostic ignored "-Wfloat-conversion"
5
6// *****************************************
8// *****************************************
9 nParams = 0;
10 SplineSegments = nullptr;
11 ParamValues = nullptr;
12}
13
14
15// *****************************************
17// *****************************************
18}
19
20// *************************
21// CW: Only need to do the binary search once per parameter, not once per event!
22// Takes down the number of binary searches from 1.2M to 17, ha!
23// For ROOT version see root/hist/hist/src/TSpline3.cxx TSpline3::FindX(double)
25// *************************
26 // Loop over the splines
27 //KS: Tried multithreading here with 48 splines and it is faster with one thread, maybe in future multithreading will be worth revisiting
28 for (M3::int_t i = 0; i < nParams; ++i)
29 {
30 const M3::int_t nPoints = SplineInfoArray[i].nPts;
31 const std::vector<M3::float_t>& xArray = SplineInfoArray[i].xPts;
32
33 // Get the variation for this reconfigure for the ith parameter
34 const float xvar = float(*SplineInfoArray[i].splineParsPointer);
35 ParamValues[i] = xvar;
36
37 // EM: if we have a parameter that has no response for any event (i.e. all splines have just one knot), then skip it and avoid a seg fault here
38 // In principle, such parameters shouldn't really be included in the first place, but with new det syst splines this
39 // could happen if say you were just running on one FHC run, then all RHC parameters would be flat and the code below would break.
40 if(xArray.size() == 0) continue;
41
42 // The segment we're interested in (klow in ROOT code)
43 M3::int_t segment = 0;
44 M3::int_t kHigh = nPoints-1;
45 //KS: We expect new segment is very close to previous
46 const M3::int_t PreviousSegment = SplineInfoArray[i].CurrSegment;
47
48 // If the variation is below the lowest saved spline point
49 if (xvar <= xArray[0]) {
50 segment = 0;
51 // If the variation is above the highest saved spline point
52 } else if (xvar >= xArray[nPoints-1]) {
53 //CW: Yes, the -2 is indeed correct, see TSpline.cxx:814 and //see: https://savannah.cern.ch/bugs/?71651
54 segment = kHigh;
55 //KS: It is quite probable the new segment is same as in previous step so try to avoid binary search, first we have to check if it is in bounds to avoid seg fault
56 } else if( xArray[PreviousSegment+1] > xvar && xvar >= xArray[PreviousSegment] ) {
57 segment = PreviousSegment;
58 // If the variation is between the maximum and minimum, perform a binary search
59 } else {
60 // The top point we've got
61 M3::int_t kHalf = 0;
62 // While there is still a difference in the points (we haven't yet found the segment)
63 // This is a binary search, incrementing segment and decrementing kHalf until we've found the segment
64 while (kHigh - segment > 1) {
65 // Increment the half-step
66 kHalf = M3::int_t((segment + kHigh)/2);
67 // If our variation is above the kHalf, set the segment to kHalf
68 if (xvar > xArray[kHalf]) {
69 segment = kHalf;
70 // Else move kHigh down
71 } else {
72 kHigh = kHalf;
73 }
74 } // End the while: we've now done our binary search
75 } // End the else: we've now found our point
76
77 if (segment >= nPoints-1 && nPoints > 1) segment = nPoints-2;
78
79 //CW: This way we avoid doing 1.2M+ binary searches on the GPU
80 // and literally just multiply lots of numbers together on the GPU without any algorithm
81 // Update the values and which segment it belongs to
82 SplineInfoArray[i].CurrSegment = segment;
83 SplineSegments[i] = short(SplineInfoArray[i].CurrSegment);
84
85 #ifdef DEBUG
86 if (SplineInfoArray[i].xPts[segment] > xvar && segment != 0) {
87 MACH3LOG_ERROR("Found a segment which is _ABOVE_ the variation!");
88 MACH3LOG_ERROR("IT SHOULD ALWAYS BE BELOW! (except when segment 0)");
89 MACH3LOG_ERROR("Spline: {}", i);
90 MACH3LOG_ERROR("Found segment = {}", segment);
91 MACH3LOG_ERROR("Doing variation = {}", xvar);
92 MACH3LOG_ERROR("x in spline = {}", SplineInfoArray[i].xPts[segment]);
93 for (M3::int_t j = 0; j < SplineInfoArray[j].nPts; ++j) {
94 MACH3LOG_ERROR(" {} = {}", j, SplineInfoArray[i].xPts[j]);
95 }
96 throw MaCh3Exception(__FILE__ , __LINE__ );
97 }
98 #endif
99 } //end loop over params
100}
101
102// *****************************************
103// Get the spline coefficients from the TSpline3 so that we can load ONLY these onto the GPU, not the whole TSpline3 object
104// This loads up coefficients into two arrays: one x array and one yabcd array
105// This should maximize our cache hits!
106void SplineBase::getTF1Coeff(TF1_red* &spl, int &nPoints, float *& coeffs) {
107// *****************************************
108 // Initialise all arrays to 1.0
109 for (int i = 0; i < _nTF1Coeff_; ++i) {
110 coeffs[i] = 0.0;
111 }
112
113 // Get number of points in spline
114 nPoints = spl->GetSize();
115
116 if(nPoints > _nTF1Coeff_)
117 {
118 MACH3LOG_ERROR("Too big number of points for TF1");
119 throw MaCh3Exception(__FILE__ , __LINE__ );
120 }
121 // TSpline3 can only take doubles, not floats
122 // But our GPU is slow with doubles, so need to cast to float
123 for (int i = 0; i < nPoints; i++) {
124 coeffs[i] = float(spl->GetParameter(M3::int_t(i)));
125 }
126 // The structure is now coeffs = {a,b,c,d,e}
127}
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:25
#define _nTF1Coeff_
KS: For TF1 we store at most 5 coefficients, we could make it more flexible but for now define it her...
Definition: SplineCommon.h:15
Custom exception class for MaCh3 errors.
virtual ~SplineBase()
Destructor.
Definition: SplineBase.cpp:16
short int nParams
Number of parameters that have splines.
Definition: SplineBase.h:64
void FindSplineSegment()
CW:Code used in step by step reweighting, Find Spline Segment for each param.
Definition: SplineBase.cpp:24
short int * SplineSegments
Definition: SplineBase.h:60
std::vector< FastSplineInfo > SplineInfoArray
Definition: SplineBase.h:57
float * ParamValues
Store parameter values they are not in FastSplineInfo as in case of GPU we need to copy paste it to G...
Definition: SplineBase.h:62
SplineBase()
Constructor.
Definition: SplineBase.cpp:7
void getTF1Coeff(TF1_red *&spl, int &nPoints, float *&coeffs)
CW: Gets the polynomial coefficients for TF1.
Definition: SplineBase.cpp:106
CW: A reduced TF1 class only. Only saves parameters for each TF1 and how many parameters each paramet...
double GetParameter(M3::int_t Parameter)
Get a parameter value.
int GetSize()
Get the size.
int int_t
Definition: Core.h:29