Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 11:16:02

0001 // -*- C++ -*-
0002 //
0003 // Package:     MVAComputer
0004 // Class  :     Spline
0005 //
0006 
0007 // Implementation:
0008 //     Simple cubic spline implementation for equidistant points in x.
0009 //
0010 // Author:      Christophe Saout
0011 // Created:     Sat Apr 24 15:18 CEST 2007
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 }  // namespace PhysicsTools