File indexing completed on 2023-03-17 11:16:02
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #include <cstring>
0015 #include <cmath>
0016
0017 #include "PhysicsTools/MVAComputer/interface/Spline.h"
0018
0019 namespace PhysicsTools {
0020
0021 float Spline::Segment::eval(float x) const {
0022 float tmp;
0023 float y = 0.0;
0024 y += coeffs[0];
0025 tmp = x;
0026 y += coeffs[1] * tmp;
0027 tmp *= x;
0028 y += coeffs[2] * tmp;
0029 tmp *= x;
0030 y += coeffs[3] * tmp;
0031 return y;
0032 }
0033
0034 float Spline::Segment::deriv(float x) const {
0035 float tmp;
0036 float d = 0.0;
0037 d += coeffs[1];
0038 tmp = x;
0039 d += coeffs[2] * tmp * 2.0;
0040 tmp *= x;
0041 d += coeffs[3] * tmp * 3.0;
0042 return d;
0043 }
0044
0045 float Spline::Segment::integral(float x) const {
0046 float tmp = x;
0047 float area = this->area;
0048 area += coeffs[0] * tmp;
0049 tmp *= x;
0050 area += coeffs[1] * tmp * (1.0 / 2.0);
0051 tmp *= x;
0052 area += coeffs[2] * tmp * (1.0 / 3.0);
0053 tmp *= x;
0054 area += coeffs[3] * tmp * (1.0 / 4.0);
0055 return area;
0056 }
0057
0058 Spline::Spline() : n(0), segments(nullptr), area(0.0) {}
0059
0060 Spline::Spline(const Spline &orig) : n(orig.n), area(orig.area) {
0061 segments = new Segment[n];
0062 std::memcpy(segments, orig.segments, sizeof(Segment) * n);
0063 }
0064
0065 Spline::Spline(unsigned int n_, const double *vals) : n(0), segments(nullptr), area(0.0) { set(n_, vals); }
0066
0067 void Spline::set(unsigned int n_, const double *vals) {
0068 n = n_ - 1;
0069 area = 0.0;
0070
0071 delete[] segments;
0072 segments = new Segment[n];
0073
0074 if (n == 1) {
0075 Segment *seg = &segments[0];
0076 seg->coeffs[0] = vals[0];
0077 seg->coeffs[1] = vals[1] - vals[0];
0078 seg->coeffs[2] = 0.0;
0079 seg->coeffs[3] = 0.0;
0080 seg->area = 0.0;
0081 area = seg->integral(1.0);
0082 return;
0083 }
0084
0085 float m0, m1;
0086 Segment *seg = &segments[0];
0087 m0 = 0.0, m1 = 0.5 * (vals[2] - vals[0]);
0088 seg->coeffs[0] = vals[0];
0089 seg->coeffs[1] = -2.0 * vals[0] + 2.0 * vals[1] - m1;
0090 seg->coeffs[2] = vals[0] - vals[1] + m1;
0091 seg->coeffs[3] = 0.0;
0092 seg->area = 0.0;
0093 area = seg->integral(1.0);
0094 m0 = m1;
0095 seg++, vals++;
0096
0097 for (unsigned int i = 1; i < n - 1; i++, seg++, vals++) {
0098 m1 = 0.5 * (vals[2] - vals[0]);
0099 seg->coeffs[0] = vals[0];
0100 seg->coeffs[1] = m0;
0101 seg->coeffs[2] = -3.0 * vals[0] - 2.0 * m0 + 3.0 * vals[1] - m1;
0102 seg->coeffs[3] = 2.0 * vals[0] + m0 - 2.0 * vals[1] + m1;
0103 seg->area = area;
0104 area = seg->integral(1.0);
0105 m0 = m1;
0106 }
0107
0108 seg->coeffs[0] = vals[0];
0109 seg->coeffs[1] = m0;
0110 seg->coeffs[2] = -vals[0] - m0 + vals[1];
0111 seg->coeffs[3] = 0.0;
0112 seg->area = area;
0113 area = seg->integral(1.0);
0114 }
0115
0116 Spline::~Spline() { delete[] segments; }
0117
0118 Spline &Spline::operator=(const Spline &orig) {
0119 delete[] segments;
0120 n = orig.n;
0121 segments = new Segment[n];
0122 std::memcpy(segments, orig.segments, sizeof(Segment) * n);
0123 area = orig.area;
0124 return *this;
0125 }
0126
0127 float Spline::eval(float x) const {
0128 if (x <= 0.0)
0129 return segments[0].eval(0.0);
0130 if (x >= 1.0)
0131 return segments[n - 1].eval(1.0);
0132
0133 float total;
0134 float rest = std::modf(x * n, &total);
0135
0136 return segments[(unsigned int)total].eval(rest);
0137 }
0138
0139 float Spline::deriv(float x) const {
0140 if (x < 0.0 || x > 1.0)
0141 return 0.0;
0142 else if (x == 0.0)
0143 return segments[0].deriv(0.0);
0144 else if (x == 1.0)
0145 return segments[n - 1].deriv(1.0);
0146
0147 float total;
0148 float rest = std::modf(x * n, &total);
0149
0150 return segments[(unsigned int)total].deriv(rest);
0151 }
0152
0153 float Spline::integral(float x) const {
0154 if (x <= 0.0)
0155 return 0.0;
0156 if (x >= 1.0)
0157 return 1.0;
0158
0159 if (area < 1.0e-9)
0160 return 0.0;
0161
0162 float total;
0163 float rest = std::modf(x * n, &total);
0164
0165 return segments[(unsigned int)total].integral(rest) / area;
0166 }
0167
0168 }
