HcalObjects/interface/InterpolatedPulse.h

0001 #ifndef CondFormats_HcalObjects_InterpolatedPulse_h_
0002 #define CondFormats_HcalObjects_InterpolatedPulse_h_
0003
0004 #include <cassert>
0005 #include "FWCore/Utilities/interface/Exception.h"
0006
0007 #include "boost/serialization/access.hpp"
0008 #include "boost/serialization/version.hpp"
0009
0010 template <unsigned MaxLen>
0011 class InterpolatedPulse {
0012   template <unsigned Len2>
0013   friend class InterpolatedPulse;
0014
0015 public:
0016   // Would normally do "static const" but genreflex has problems for it
0017   enum { maxlen = MaxLen };
0018
0019   // Default constructor creates a pulse which is zero everywhere
0020   inline InterpolatedPulse() : tmin_(0.0), width_(1.0), length_(2U) { zeroOut(); }
0021
0022   // Constructor from a single integer creates a pulse with the given
0023   // number of discrete steps which is zero everywhere
0024   inline explicit InterpolatedPulse(const unsigned len) : tmin_(0.0), width_(1.0), length_(len) {
0025     if (length_ < 2 || length_ > MaxLen)
0026       throw cms::Exception("In InterpolatedPulse constructor: invalid length");
0027     zeroOut();
0028   }
0029
0030   inline InterpolatedPulse(const double tmin, const double tmax, const unsigned len)
0031       : tmin_(tmin), width_(tmax - tmin), length_(len) {
0032     if (length_ < 2 || length_ > MaxLen)
0033       throw cms::Exception("In InterpolatedPulse constructor: invalid length");
0034     if (width_ <= 0.0)
0035       throw cms::Exception("In InterpolatedPulse constructor: invalid pulse width");
0036     zeroOut();
0037   }
0038
0039   template <typename Real>
0040   inline InterpolatedPulse(const double tmin, const double tmax, const Real* values, const unsigned len)
0041       : tmin_(tmin), width_(tmax - tmin) {
0042     if (width_ <= 0.0)
0043       throw cms::Exception("In InterpolatedPulse constructor: invalid pulse width");
0044     setShape(values, len);
0045   }
0046
0047   // Efficient copy constructor. Do not copy undefined values.
0048   inline InterpolatedPulse(const InterpolatedPulse& r) : tmin_(r.tmin_), width_(r.width_), length_(r.length_) {
0049     double* buf = &pulse_[0];
0050     const double* rbuf = &r.pulse_[0];
0051     for (unsigned i = 0; i < length_; ++i)
0052       *buf++ = *rbuf++;
0053   }
0054
0055   // Converting copy constructor
0056   template <unsigned Len2>
0057   inline InterpolatedPulse(const InterpolatedPulse<Len2>& r) : tmin_(r.tmin_), width_(r.width_), length_(r.length_) {
0058     if (length_ > MaxLen)
0059       throw cms::Exception("In InterpolatedPulse copy constructor: buffer is not long enough");
0060     double* buf = &pulse_[0];
0061     const double* rbuf = &r.pulse_[0];
0062     for (unsigned i = 0; i < length_; ++i)
0063       *buf++ = *rbuf++;
0064   }
0065
0066   // Efficient assignment operator
0067   inline InterpolatedPulse& operator=(const InterpolatedPulse& r) {
0068     if (this != &r) {
0069       tmin_ = r.tmin_;
0070       width_ = r.width_;
0071       length_ = r.length_;
0072       double* buf = &pulse_[0];
0073       const double* rbuf = &r.pulse_[0];
0074       for (unsigned i = 0; i < length_; ++i)
0075         *buf++ = *rbuf++;
0076     }
0077     return *this;
0078   }
0079
0080   // Converting assignment operator
0081   template <unsigned Len2>
0082   inline InterpolatedPulse& operator=(const InterpolatedPulse<Len2>& r) {
0083     if (r.length_ > MaxLen)
0084       throw cms::Exception("In InterpolatedPulse::operator=: buffer is not long enough");
0085     tmin_ = r.tmin_;
0086     width_ = r.width_;
0087     length_ = r.length_;
0088     double* buf = &pulse_[0];
0089     const double* rbuf = &r.pulse_[0];
0090     for (unsigned i = 0; i < length_; ++i)
0091       *buf++ = *rbuf++;
0092     return *this;
0093   }
0094
0095   template <typename Real>
0096   inline void setShape(const Real* values, const unsigned len) {
0097     if (len < 2 || len > MaxLen)
0098       throw cms::Exception("In InterpolatedPulse::setShape: invalid length");
0099     assert(values);
0100     length_ = len;
0101     double* buf = &pulse_[0];
0102     for (unsigned i = 0; i < len; ++i)
0103       *buf++ = *values++;
0104   }
0105
0106   // Zero out the signal
0107   inline void zeroOut() {
0108     for (unsigned i = 0; i < length_; ++i)
0109       pulse_[i] = 0.0;
0110   }
0111
0112   // Simple inspectors
0113   inline const double* getPulse() const { return &pulse_[0]; }
0114   inline unsigned getLength() const { return length_; }
0115   inline double getStartTime() const { return tmin_; }
0116   inline double getStopTime() const { return tmin_ + width_; }
0117   inline double getPulseWidth() const { return width_; }
0118   inline double getTimeStep() const { return width_ / (length_ - 1U); }
0119
0120   // Simple modifiers
0121   inline void setStartTime(const double newStartTime) { tmin_ = newStartTime; }
0122
0123   inline void setPulseWidth(const double newWidth) {
0124     if (newWidth <= 0.0)
0125       throw cms::Exception("In InterpolatedPulse::setPulseWidth: invalid pulse width");
0126     width_ = newWidth;
0127   }
0128
0129   // Get the pulse value at the given time using linear interpolation
0130   inline double operator()(const double t) const {
0131     const volatile double tmax = tmin_ + width_;
0132     if (t < tmin_ || t > tmax)
0133       return 0.0;
0134     const unsigned lm1 = length_ - 1U;
0135     const double step = width_ / lm1;
0136     const double nSteps = (t - tmin_) / step;
0137     unsigned nbelow = nSteps;
0138     unsigned nabove = nbelow + 1;
0139     if (nabove > lm1) {
0140       nabove = lm1;
0141       nbelow = nabove - 1U;
0142     }
0143     const double delta = nSteps - nbelow;
0144     return pulse_[nbelow] * (1.0 - delta) + pulse_[nabove] * delta;
0145   }
0146
0147   inline double derivative(const double t) const {
0148     const volatile double tmax = tmin_ + width_;
0149     if (t < tmin_ || t > tmax)
0150       return 0.0;
0151     const unsigned lm1 = length_ - 1U;
0152     const double step = width_ / lm1;
0153     const double nSteps = (t - tmin_) / step;
0154     unsigned nbelow = nSteps;
0155     unsigned nabove = nbelow + 1;
0156     if (nabove > lm1) {
0157       nabove = lm1;
0158       nbelow = nabove - 1U;
0159     }
0160     const double delta = nSteps - nbelow;
0161     if ((nbelow == 0U && delta <= 0.5) || (nabove == lm1 && delta >= 0.5))
0162       return (pulse_[nabove] - pulse_[nbelow]) / step;
0163     else if (delta >= 0.5) {
0164       const double lower = pulse_[nabove] - pulse_[nbelow];
0165       const double upper = pulse_[nabove + 1U] - pulse_[nabove];
0166       return (upper * (delta - 0.5) + lower * (1.5 - delta)) / step;
0167     } else {
0168       const double lower = pulse_[nbelow] - pulse_[nbelow - 1U];
0169       const double upper = pulse_[nabove] - pulse_[nbelow];
0170       return (lower * (0.5 - delta) + upper * (0.5 + delta)) / step;
0171     }
0172   }
0173
0174   inline double secondDerivative(const double t) const {
0175     const volatile double tmax = tmin_ + width_;
0176     if (t < tmin_ || t > tmax || length_ < 3U)
0177       return 0.0;
0178     const unsigned lm1 = length_ - 1U;
0179     const double step = width_ / lm1;
0180     const double stepSq = step * step;
0181     const double nSteps = (t - tmin_) / step;
0182     unsigned nbelow = nSteps;
0183     unsigned nabove = nbelow + 1;
0184     if (nabove > lm1) {
0185       nabove = lm1;
0186       nbelow = nabove - 1U;
0187     }
0188
0189     if (nbelow == 0U) {
0190       // The first interval
0191       return (pulse_[2] - 2.0 * pulse_[1] + pulse_[0]) / stepSq;
0192     } else if (nabove == lm1) {
0193       // The last interval
0194       return (pulse_[lm1] - 2.0 * pulse_[lm1 - 1U] + pulse_[lm1 - 2U]) / stepSq;
0195     } else {
0196       // One of the middle intervals
0197       const double lower = pulse_[nbelow - 1U] - 2.0 * pulse_[nbelow] + pulse_[nabove];
0198       const double upper = pulse_[nbelow] - 2.0 * pulse_[nabove] + pulse_[nabove + 1U];
0199       const double delta = nSteps - nbelow;
0200       return (lower * (1.0 - delta) + upper * delta) / stepSq;
0201     }
0202   }
0203
0204   inline InterpolatedPulse& operator*=(const double scale) {
0205     if (scale != 1.0) {
0206       double* buf = &pulse_[0];
0207       for (unsigned i = 0; i < length_; ++i)
0208         *buf++ *= scale;
0209     }
0210     return *this;
0211   }
0212
0213   // Add another pulse to this one. Note that addition of another pulse
0214   // will not change the start time or the width of _this_ pulse. The
0215   // added pulse will be truncated as needed.
0216   template <unsigned Len2>
0217   inline InterpolatedPulse& operator+=(const InterpolatedPulse<Len2>& r) {
0218     const double step = width_ / (length_ - 1U);
0219     for (unsigned i = 0; i < length_; ++i)
0220       pulse_[i] += r(tmin_ + i * step);
0221     return *this;
0222   }
0223
0224   template <unsigned Len2>
0225   inline bool operator==(const InterpolatedPulse<Len2>& r) const {
0226     if (!(tmin_ == r.tmin_ && width_ == r.width_ && length_ == r.length_))
0227       return false;
0228     const double* buf = &pulse_[0];
0229     const double* rbuf = &r.pulse_[0];
0230     for (unsigned i = 0; i < length_; ++i)
0231       if (*buf++ != *rbuf++)
0232         return false;
0233     return true;
0234   }
0235
0236   template <unsigned Len2>
0237   inline bool operator!=(const InterpolatedPulse<Len2>& r) const {
0238     return !(*this == r);
0239   }
0240
0241   // Simple trapezoidal integration
0242   inline double getIntegral() const {
0243     const double* buf = &pulse_[0];
0244     long double sum = buf[0] / 2.0;
0245     const unsigned nIntervals = length_ - 1U;
0246     for (unsigned i = 1U; i < nIntervals; ++i)
0247       sum += buf[i];
0248     sum += buf[nIntervals] / 2.0;
0249     return sum * width_ / nIntervals;
0250   }
0251
0252   inline void setIntegral(const double newValue) {
0253     const double integ = this->getIntegral();
0254     if (integ == 0.0)
0255       throw cms::Exception("In InterpolatedPulse::setIntegral division by zero");
0256     *this *= (newValue / integ);
0257   }
0258
0259   inline double getPeakValue() const {
0260     const double* buf = &pulse_[0];
0261     double peak = buf[0];
0262     for (unsigned i = 1U; i < length_; ++i)
0263       if (buf[i] > peak)
0264         peak = buf[i];
0265     return peak;
0266   }
0267
0268   inline void setPeakValue(const double newValue) {
0269     const double peak = this->getPeakValue();
0270     if (peak == 0.0)
0271       throw cms::Exception("In InterpolatedPulse::setPeakValue: division by zero");
0272     *this *= (newValue / peak);
0273   }
0274
0275 private:
0276   double pulse_[MaxLen];
0277   double tmin_;
0278   double width_;
0279   unsigned length_;
0280
0281   friend class boost::serialization::access;
0282
0283   template <class Archive>
0284   inline void serialize(Archive& ar, unsigned /* version */) {
0285     ar & tmin_ & width_ & length_;
0286
0287     // In case we are reading, it may be useful to verify
0288     // that the length is reasonable
0289     if (length_ > MaxLen)
0290       throw cms::Exception("In InterpolatedPulse::serialize: buffer is not long enough");
0291
0292     for (unsigned i = 0; i < length_; ++i)
0293       ar& pulse_[i];
0294   }
0295 };
0296
0297 // boost serialization version number for this template
0298 namespace boost {
0299   namespace serialization {
0300     template <unsigned MaxLen>
0301     struct version<InterpolatedPulse<MaxLen> > {
0302       BOOST_STATIC_CONSTANT(int, value = 1);
0303     };
0304   }  // namespace serialization
0305 }  // namespace boost
0306
0307 #endif  // CondFormats_HcalObjects_InterpolatedPulse_h_