MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
SplineStructs.h
Go to the documentation of this file.
1#pragma once
2
3// MaCh3 includes
6
7#include <cmath>
8
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"
13
18
19// *******************
22// *******************
25 nPts = -999;
26 CurrSegment = 0;
27 splineParsPointer = nullptr;
28 }
29
31 virtual ~FastSplineInfo() {
32
33 }
34
37
39 std::vector<M3::float_t> xPts;
40
43
45 const double* splineParsPointer;
46};
47
48
49// ***************************************************************************
51inline void ApplyKnotWeightCap(TGraph* xsecgraph, const int splineParsIndex, ParameterHandlerGeneric* XsecCov) {
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}
97
98// ***************************************************************************
100inline void ApplyKnotWeightCapTSpline3(TSpline3* &Spline, const int splineParsIndex, ParameterHandlerGeneric* XsecCov) {
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}
130
131// ************************
134// ************************
135public:
141 virtual double Eval(const double var)=0;
143 virtual void Print()=0;
145 virtual M3::int_t GetNp()=0;
146};
147
148// ************************
151// ************************
152public:
155 length = 0;
156 Par = NULL;
157 }
158
160 virtual ~TF1_red() {
161 if (Par != NULL) {
162 delete[] Par;
163 Par = NULL;
164 }
165 }
166
169 length = nSize;
170 for (int i = 0; i < length; ++i) {
171 Par[i] = Array[i];
172 }
173 }
174
176 TF1_red(TF1* &Function) : TResponseFunction_red() {
177 Par = NULL;
178 SetFunc(Function);
179 }
180
182 inline void SetFunc(TF1* &Func) {
183 length = M3::int_t(Func->GetNpar());
184 if (Par != NULL) delete[] Par;
185 Par = new M3::float_t[length];
186 for (int i = 0; i < length; ++i) {
187 Par[i] = M3::float_t(Func->GetParameter(i));
188 }
189 delete Func;
190 Func = NULL;
191 }
192
194 inline double Eval(const double var) override {
195 return Par[1]+Par[0]*var;
196
197 /* FIXME in future we might introduce more TF1
198 //If we have 5 parameters we're using a fifth order polynomial
199 if (Type == kFifthOrderPolynomial) {
200 return 1+Par[0]*var+Par[1]*var*var+Par[2]*var*var*var+Par[3]*var*var*var*var+Par[4]*var*var*var*var*var;
201 } else if (Type == kTwoLinears) {
202 return (var <= 0)*(Par[2]+Par[0]*var)+(var > 0)*(Par[2]+Par[1]*var);
203 } else if (Type == kLinear) {
204 return (Par[1]+Par[0]*var);
205 } else if (Type == kPseudoHeaviside) {
206 return (var <= 0)*(1+Par[0]*var) + (1 >= var)*(var > 0)*(1+Par[1]*var) + (var > 1)*(Par[3]+Par[2]*var);
207 }else {
208 MACH3LOG_ERROR(" Class only knows about 5th order polynomial, two superposed linear functions, linear function, or pseudo Heaviside.");
209 MACH3LOG_ERROR(" You have tried something else than this, which remains unimplemented.");
210 MACH3LOG_ERROR("{}: {}", __FILE__, __LINE__);
211 throw MaCh3Exception(__FILE__ , __LINE__ );
212 }
213 */
214 }
215
217 inline void SetParameter(M3::int_t Parameter, M3::float_t Value) {
218 Par[Parameter] = Value;
219 }
220
222 double GetParameter(M3::int_t Parameter) {
223 if (Parameter > length) {
224 MACH3LOG_ERROR("You requested parameter number {} but length is {} parameters", Parameter, length);
225 throw MaCh3Exception(__FILE__ , __LINE__ );
226 return -999.999;
227 }
228 return Par[Parameter];
229 }
230
232 inline void SetSize(M3::int_t nSpline) {
233 length = nSpline;
234 Par = new M3::float_t[length];
235 }
237 inline int GetSize() { return length; }
239 inline void Print() override {
240 MACH3LOG_INFO("Printing TF1_red:");
241 MACH3LOG_INFO(" Length = {}", length);
242 for (int i = 0; i < length; i++) {
243 MACH3LOG_INFO(" Coeff {} = {}", i, Par[i]);
244 }
245 }
246
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]);
252 }
253 return func;
254 }
255
257 inline M3::int_t GetNp() override { return length; }
258
259private:
263};
264
265// ************************
268// ************************
269public:
272 nPoints = 0;
273 Par = NULL;
274 XPos = NULL;
275 YResp = NULL;
276 }
277
279 TSpline3_red(TSpline3* &spline, SplineInterpolation InterPolation = kTSpline3) : TResponseFunction_red() {
280 Par = NULL;
281 XPos = NULL;
282 YResp = NULL;
283 SetFunc(spline, InterPolation);
284 }
285
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
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 }
312 inline void SetFunc(TSpline3* &spline, SplineInterpolation InterPolation = kTSpline3) {
313 nPoints = M3::int_t(spline->GetNp());
314 if (Par != NULL) {
315 for (int i = 0; i < nPoints; ++i) {
316 delete[] Par[i];
317 Par[i] = NULL;
318 }
319 delete[] Par;
320 Par = NULL;
321 }
322 if (XPos != NULL) delete[] XPos;
323 if (YResp != NULL) delete[] YResp;
324 // Save the parameters for each knot
325 Par = new M3::float_t*[nPoints];
326 // Save the positions of the knots
327 XPos = new M3::float_t[nPoints];
328 // Save the y response at each knot
330
331 //KS: Default TSpline3 ROOT implementation
332 if(InterPolation == kTSpline3)
333 {
334 for (int i = 0; i < nPoints; ++i) {
335 // 3 is the size of the TSpline3 coefficients
336 Par[i] = new M3::float_t[3];
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);
339 XPos[i] = M3::float_t(x);
340 YResp[i] = M3::float_t(y);
341 Par[i][0] = M3::float_t(b);
342 Par[i][1] = M3::float_t(c);
343 Par[i][2] = M3::float_t(d);
344 }
345 }
346 //CW: Reduce to use linear spline interpolation for certain parameters
347 // 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
348 // Get which parameters should be linear from the fit manager
349 // Convert the spline number to global xsec parameter
350 // Loop over the splines points
351 // KS: kLinearFunc should be used with TF1, this is just as safety
352 else if(InterPolation == kLinear || InterPolation == kLinearFunc)
353 {
354 for (int k = 0; k < nPoints; ++k) {
355 // 3 is the size of the TSpline3 coefficients
356 Par[k] = new M3::float_t[3];
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);
361
362 XPos[k] = M3::float_t(x1);
363 YResp[k] = M3::float_t(y1);
364 Par[k][0] = M3::float_t(tempb);
365 Par[k][1] = M3::float_t(0);
366 Par[k][2] = M3::float_t(0);
367 }
368 }
369 //EM: Akima spline is similar to regular cubic spline but is allowed to be discontinuous in 2nd derivative and coefficients in any segment
370 // only depend on th 2 nearest points on either side
371 else if(InterPolation == kAkima)
372 {
373 // get the knot values for the spline
374 for (int i = 0; i < nPoints; ++i) {
375 // 3 is the size of the TSpline3 coefficients
376 Par[i] = new M3::float_t[3];
377
378 double x = -999.99, y = -999.99;
379 spline->GetKnot(i, x, y);
380
381 XPos[i] = M3::float_t(x);
382 YResp[i] = M3::float_t(y);
383 }
384
385 M3::float_t* mvals = new M3::float_t[nPoints + 3];
386 M3::float_t* svals = new M3::float_t[nPoints + 1];
387
388 for (int i = -2; i <= nPoints; ++i) {
389 // if segment is first or last or 2nd to first or last, needs to be dealt with slightly differently;
390 // need to estimate the values for additinal points which would lie outside of the spline
391 if(i ==-2){
392 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]));
393 }
394 else if(i==-1){
395 mvals[i+2] = M3::float_t(2.0 * (YResp[1] - YResp[0]) / (XPos[1] - XPos[0]) - (YResp[2] - YResp[1]) / (XPos[2] - XPos[1]));
396 }
397 else if(i==nPoints){
398 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]));
399 }
400 else if(i == nPoints - 1){
401 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]));
402 }
403 //standard internal segment
404 else{
405 mvals[i+2] = (YResp[i+1] - YResp[i])/ (XPos[i+1] - XPos[i]);
406 }
407 }
408
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]));
412 }
413 else{svals[i-2] = mvals[i];}
414 }
415
416 // calculate the coefficients for the spline
417 for(int i = 0; i <nPoints; i++){
418 M3::float_t b, c, d = M3::float_t(-999.999);
419
420 b = svals[i];
421 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]);
422 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]));
423
424 Par[i][0] = b;
425 Par[i][1] = c;
426 Par[i][2] = d;
427 }
428
429 // check the input spline for linear segments, if there are any then overwrite the calculated coefficients
430 // 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
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);
434
435 if((c == 0.0 && d == 0.0)){
436 Par[i][0] = M3::float_t(b);
437 Par[i][1] = M3::float_t(0.0);
438 Par[i][2] = M3::float_t(0.0);
439 }
440 }
441 delete[] mvals;
442 delete[] svals;
443 }
444 //EM: Monotone spline is similar to regular cubic spline but enforce the condition that the interpolated value at any point
445 // must be between its two nearest knots, DOES NOT make the entire spline monotonic, only the segments
446 else if(InterPolation == kMonotonic)
447 {
448 // values of the secants at each point (for calculating monotone spline)
449 M3::float_t * Secants = new M3::float_t[nPoints -1];
450 // values of the tangens at each point (for calculating monotone spline)
451 M3::float_t * Tangents = new M3::float_t[nPoints];
452
453 // get the knot values for the spline
454 for (int i = 0; i < nPoints; ++i) {
455 // 3 is the size of the TSpline3 coefficients
456 Par[i] = new M3::float_t[3];
457
458 double x = -999.99, y = -999.99;
459 spline->GetKnot(i, x, y);
460
461 XPos[i] = M3::float_t(x);
462 YResp[i] = M3::float_t(y);
463
464 Tangents[i] = 0.0;
465 }
466
467 // deal with the case of two points (just do linear interpolation between them)
468 if (nPoints == 2){
469 Par[0][0] = (YResp[1] - YResp[0]) / (XPos[1] - XPos[0]);
470 Par[0][1] = 0.0;
471 Par[0][2] = 0.0;
472 // extra "virtual" segment at end to make Par array shape fit with knot arrays shapes
473 Par[1][1] = 0.0;
474 Par[1][2] = 0.0;
475
476 return;
477 } // if nPoints !=2 do full monotonic spline treatment:
478 else
479 {
480 // first pass over knots to calculate the secants
481 for (int i = 0; i < nPoints-1; ++i) {
482 Secants[i] = (YResp[i+1] - YResp[i]) / (XPos[i+1] - XPos[i]);
483 MACH3LOG_TRACE("Secant {}: {}", i, Secants[i]);
484 }
485
486 Tangents[0] = Secants[0];
487 Tangents[nPoints-1] = Secants[nPoints -2];
488
489 M3::float_t alpha;
490 M3::float_t beta;
491
492 // second pass over knots to calculate tangents
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)){ //check for same sign
495 Tangents[i] = M3::float_t((Secants[i-1] + Secants[i]) /2.0);
496 }
497 }
498
499 // third pass over knots to rescale tangents
500 for (int i = 0; i < nPoints-1; ++i) {
501 if (Secants[i] == 0.0){
502 Tangents[i] = 0.0;
503 Tangents[i+1] = 0.0;
504 }
505
506 else{
507 alpha = Tangents[i] / Secants[i];
508 beta = Tangents[i+1] / Secants[i];
509
510 if (alpha <0.0){
511 Tangents[i] = 0.0;
512 }
513 if (beta < 0.0){
514 Tangents[i+1] = 0.0;
515 }
516
517 if (alpha * alpha + beta * beta >9.0){
518 M3::float_t tau = M3::float_t(3.0 / std::sqrt(alpha * alpha + beta * beta));
519 Tangents[i] = tau * alpha * Secants[i];
520 Tangents[i+1] = tau * beta * Secants[i];
521 }
522 }
523 } // finished rescaling tangents
524 // fourth pass over knots to calculate the coefficients for the spline
525 M3::float_t dx;
526 for(int i = 0; i <nPoints-1; i++){
527 M3::float_t b, c, d = M3::float_t(-999.999);
528 dx = XPos[i+1] - XPos[i];
529
530 b = Tangents[i] * dx;
531 c = M3::float_t(3.0* (YResp[i+1] - YResp[i]) -2.0 *dx * Tangents[i] - dx * Tangents[i +1]);
532 d = M3::float_t(2.0* (YResp[i] - YResp[i+1]) + dx * (Tangents[i] + Tangents[i+1]));
533
534 Par[i][0] = b / dx;
535 Par[i][1] = c / (dx * dx);
536 Par[i][2] = d / (dx * dx * dx);
537
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.",
540 i, Par[i][0], Par[i][1], Par[i][2]);
541 }
542 MACH3LOG_TRACE("b: {}", b);
543 MACH3LOG_TRACE("dx: {}, x_0: {}, x_1: {}", dx, XPos[i], XPos[i+1]);
544 MACH3LOG_TRACE(" y_0: {}, y_1: {}", YResp[i], YResp[i+1]);
545 }
546
547 // include params for final "segment" outside of the spline so that par array fits with x and y arrays,
548 // should never actually get used but if not set then the GPU code gets very angry
549 Par[nPoints-1][0] = 0.0;
550 Par[nPoints-1][1] = 0.0;
551 Par[nPoints-1][2] = 0.0;
552
553 // check the input spline for linear segments, if there are any then overwrite the calculated coefficients
554 // 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
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);
558
559 if((c == 0.0 && d == 0.0)){
560 Par[i][0] = M3::float_t(b);
561 Par[i][1] = 0.0;
562 Par[i][2] = 0.0;
563 }
564 }
565 delete[] Secants;
566 delete[] Tangents;
567 } // end of if(nPoints !=2)
568 }
569 else
570 {
571 MACH3LOG_ERROR("Unsupported interpolation type {}", static_cast<int>(InterPolation));
572 throw MaCh3Exception(__FILE__ , __LINE__ );
573 }
574
575 delete spline;
576 spline = NULL;
577 }
578
580 virtual ~TSpline3_red() {
581 if(Par != NULL) {
582 for (int i = 0; i < nPoints; ++i) {
583 if (Par[i] != NULL) {
584 delete[] Par[i];
585 }
586 }
587 delete[] Par;
588 }
589 if(XPos != NULL) delete[] XPos;
590 if(YResp != NULL) delete[] YResp;
591 Par = NULL;
592 XPos = YResp = NULL;
593 }
594
597 inline int FindX(double x) {
598 // The segment we're interested in (klow in ROOT code)
599 int segment = 0;
600 int kHigh = nPoints-1;
601 // If the variation is below the lowest saved spline point
602 if (x <= XPos[0]){
603 segment = 0;
604 // If the variation is above the highest saved spline point
605 } else if (x >= XPos[nPoints-1]) {
606 // Yes, the -2 is indeed correct, see TSpline.cxx:814 and //see: https://savannah.cern.ch/bugs/?71651
607 segment = kHigh;
608 // If the variation is between the maximum and minimum, perform a binary search
609 } else {
610 // The top point we've got
611 int kHalf = 0;
612 // While there is still a difference in the points (we haven't yet found the segment)
613 // This is a binary search, incrementing segment and decrementing kHalf until we've found the segment
614 while (kHigh - segment > 1) {
615 // Increment the half-step
616 kHalf = (segment + kHigh)/2;
617 // If our variation is above the kHalf, set the segment to kHalf
618 if (x > XPos[kHalf]) {
619 segment = kHalf;
620 // Else move kHigh down
621 } else {
622 kHigh = kHalf;
623 }
624 } // End the while: we've now done our binary search
625 } // End the else: we've now found our point
626 if (segment >= nPoints-1 && nPoints > 1) segment = nPoints-2;
627 return segment;
628 }
629
631 inline double Eval(double var) override {
632 // Get the segment for this variation
633 int segment = FindX(var);
634 // The get the coefficients for this variation
635 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);
636 GetCoeff(segment, x, y, b, c, d);
637 double dx = var - x;
638 // Evaluate the third order polynomial
639 double weight = y+dx*(b+dx*(c+d*dx));
640 return weight;
641 }
642
644 inline M3::int_t GetNp() override { return nPoints; }
645 // Get the ith knot's x and y position
646 inline void GetKnot(int i, M3::float_t &xtmp, M3::float_t &ytmp) {
647 xtmp = XPos[i];
648 ytmp = YResp[i];
649 }
650
652 inline void GetCoeff(int segment, M3::float_t &x, M3::float_t &y, M3::float_t &b, M3::float_t &c, M3::float_t &d) {
653 b = Par[segment][0];
654 c = Par[segment][1];
655 d = Par[segment][2];
656 x = XPos[segment];
657 y = YResp[segment];
658 }
659
661 inline TSpline3* ConstructTSpline3() {
662 // KS: Sadly ROOT only accepts double...
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]); // Convert float to double
668 yPosDoubles[i] = static_cast<Double_t>(YResp[i]); // Convert float to double
669 }
670 TSpline3 *spline = new TSpline3("Spline", xPosDoubles.data(), yPosDoubles.data(), static_cast<int>(nPoints));
671 #else
672 TSpline3 *spline = new TSpline3("Spline", XPos, YResp, nPoints);
673 #endif
674 return spline;
675 }
676
678 inline void Print() override {
679 MACH3LOG_INFO("Printing TSpline_red:");
680 MACH3LOG_INFO(" Nknots = {}", nPoints);
681 for (int i = 0; i < nPoints; ++i) {
682 MACH3LOG_INFO(" i = {} x = {} y = {} b = {} c = {} d = {}",
683 i, XPos[i], YResp[i], Par[i][0], Par[i][1], Par[i][2]);
684 }
685 }
686
687 protected: //changed to protected from private so can be accessed by derived classes
696};
697
698// *****************************************
701inline bool isFlat(TSpline3_red* &spl) {
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}
716
717// *********************************
720inline std::vector<std::vector<TSpline3_red*> > ReduceTSpline3(std::vector<std::vector<TSpline3*> > &MasterSpline) {
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}
754
755// *********************************
758inline std::vector<std::vector<TF1_red*> > ReduceTF1(std::vector<std::vector<TF1*> > &MasterSpline) {
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}
792
793// *********************************
800 const RespFuncType SplineRespFuncType,
801 const SplineInterpolation SplineInterpolationType,
802 const std::string& Title) {
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}
863
864#pragma GCC diagnostic pop
865
std::vector< bool > isFlat
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:25
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:23
#define MACH3LOG_TRACE
Definition: MaCh3Logger.h:21
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...
Definition: SplineStructs.h:51
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.
M3::int_t length
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.
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
double float_t
Definition: Core.h:28
int int_t
Definition: Core.h:29
CW: Add a struct to hold info about the splinified xsec parameters and help with FindSplineSegment.
Definition: SplineStructs.h:21
virtual ~FastSplineInfo()
Destructor.
Definition: SplineStructs.h:31
std::vector< M3::float_t > xPts
Array of the knots positions.
Definition: SplineStructs.h:39
const double * splineParsPointer
Array of the knots positions.
Definition: SplineStructs.h:45
M3::int_t CurrSegment
Array of what segment of spline we're currently interested in. Gets updated once per MCMC iteration.
Definition: SplineStructs.h:42
FastSplineInfo()
Constructor.
Definition: SplineStructs.h:24
M3::int_t nPts
Number of points in spline.
Definition: SplineStructs.h:36