CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
InterpolatedPulse.h
Go to the documentation of this file.
1 #ifndef CondFormats_HcalObjects_InterpolatedPulse_h_
2 #define CondFormats_HcalObjects_InterpolatedPulse_h_
3 
4 #include <cassert>
6 
7 #include "boost/serialization/access.hpp"
8 #include "boost/serialization/version.hpp"
9 
10 template <unsigned MaxLen>
12  template <unsigned Len2>
13  friend class InterpolatedPulse;
14 
15 public:
16  // Would normally do "static const" but genreflex has problems for it
17  enum { maxlen = MaxLen };
18 
19  // Default constructor creates a pulse which is zero everywhere
20  inline InterpolatedPulse() : tmin_(0.0), width_(1.0), length_(2U) { zeroOut(); }
21 
22  // Constructor from a single integer creates a pulse with the given
23  // number of discrete steps which is zero everywhere
24  inline explicit InterpolatedPulse(const unsigned len) : tmin_(0.0), width_(1.0), length_(len) {
25  if (length_ < 2 || length_ > MaxLen)
26  throw cms::Exception("In InterpolatedPulse constructor: invalid length");
27  zeroOut();
28  }
29 
30  inline InterpolatedPulse(const double tmin, const double tmax, const unsigned len)
31  : tmin_(tmin), width_(tmax - tmin), length_(len) {
32  if (length_ < 2 || length_ > MaxLen)
33  throw cms::Exception("In InterpolatedPulse constructor: invalid length");
34  if (width_ <= 0.0)
35  throw cms::Exception("In InterpolatedPulse constructor: invalid pulse width");
36  zeroOut();
37  }
38 
39  template <typename Real>
40  inline InterpolatedPulse(const double tmin, const double tmax, const Real* values, const unsigned len)
41  : tmin_(tmin), width_(tmax - tmin) {
42  if (width_ <= 0.0)
43  throw cms::Exception("In InterpolatedPulse constructor: invalid pulse width");
44  setShape(values, len);
45  }
46 
47  // Efficient copy constructor. Do not copy undefined values.
49  double* buf = &pulse_[0];
50  const double* rbuf = &r.pulse_[0];
51  for (unsigned i = 0; i < length_; ++i)
52  *buf++ = *rbuf++;
53  }
54 
55  // Converting copy constructor
56  template <unsigned Len2>
58  if (length_ > MaxLen)
59  throw cms::Exception("In InterpolatedPulse copy constructor: buffer is not long enough");
60  double* buf = &pulse_[0];
61  const double* rbuf = &r.pulse_[0];
62  for (unsigned i = 0; i < length_; ++i)
63  *buf++ = *rbuf++;
64  }
65 
66  // Efficient assignment operator
68  if (this != &r) {
69  tmin_ = r.tmin_;
70  width_ = r.width_;
71  length_ = r.length_;
72  double* buf = &pulse_[0];
73  const double* rbuf = &r.pulse_[0];
74  for (unsigned i = 0; i < length_; ++i)
75  *buf++ = *rbuf++;
76  }
77  return *this;
78  }
79 
80  // Converting assignment operator
81  template <unsigned Len2>
83  if (r.length_ > MaxLen)
84  throw cms::Exception("In InterpolatedPulse::operator=: buffer is not long enough");
85  tmin_ = r.tmin_;
86  width_ = r.width_;
87  length_ = r.length_;
88  double* buf = &pulse_[0];
89  const double* rbuf = &r.pulse_[0];
90  for (unsigned i = 0; i < length_; ++i)
91  *buf++ = *rbuf++;
92  return *this;
93  }
94 
95  template <typename Real>
96  inline void setShape(const Real* values, const unsigned len) {
97  if (len < 2 || len > MaxLen)
98  throw cms::Exception("In InterpolatedPulse::setShape: invalid length");
99  assert(values);
100  length_ = len;
101  double* buf = &pulse_[0];
102  for (unsigned i = 0; i < len; ++i)
103  *buf++ = *values++;
104  }
105 
106  // Zero out the signal
107  inline void zeroOut() {
108  for (unsigned i = 0; i < length_; ++i)
109  pulse_[i] = 0.0;
110  }
111 
112  // Simple inspectors
113  inline const double* getPulse() const { return &pulse_[0]; }
114  inline unsigned getLength() const { return length_; }
115  inline double getStartTime() const { return tmin_; }
116  inline double getStopTime() const { return tmin_ + width_; }
117  inline double getPulseWidth() const { return width_; }
118  inline double getTimeStep() const { return width_ / (length_ - 1U); }
119 
120  // Simple modifiers
121  inline void setStartTime(const double newStartTime) { tmin_ = newStartTime; }
122 
123  inline void setPulseWidth(const double newWidth) {
124  if (newWidth <= 0.0)
125  throw cms::Exception("In InterpolatedPulse::setPulseWidth: invalid pulse width");
126  width_ = newWidth;
127  }
128 
129  // Get the pulse value at the given time using linear interpolation
130  inline double operator()(const double t) const {
131  const volatile double tmax = tmin_ + width_;
132  if (t < tmin_ || t > tmax)
133  return 0.0;
134  const unsigned lm1 = length_ - 1U;
135  const double step = width_ / lm1;
136  const double nSteps = (t - tmin_) / step;
137  unsigned nbelow = nSteps;
138  unsigned nabove = nbelow + 1;
139  if (nabove > lm1) {
140  nabove = lm1;
141  nbelow = nabove - 1U;
142  }
143  const double delta = nSteps - nbelow;
144  return pulse_[nbelow] * (1.0 - delta) + pulse_[nabove] * delta;
145  }
146 
147  inline double derivative(const double t) const {
148  const volatile double tmax = tmin_ + width_;
149  if (t < tmin_ || t > tmax)
150  return 0.0;
151  const unsigned lm1 = length_ - 1U;
152  const double step = width_ / lm1;
153  const double nSteps = (t - tmin_) / step;
154  unsigned nbelow = nSteps;
155  unsigned nabove = nbelow + 1;
156  if (nabove > lm1) {
157  nabove = lm1;
158  nbelow = nabove - 1U;
159  }
160  const double delta = nSteps - nbelow;
161  if ((nbelow == 0U && delta <= 0.5) || (nabove == lm1 && delta >= 0.5))
162  return (pulse_[nabove] - pulse_[nbelow]) / step;
163  else if (delta >= 0.5) {
164  const double lower = pulse_[nabove] - pulse_[nbelow];
165  const double upper = pulse_[nabove + 1U] - pulse_[nabove];
166  return (upper * (delta - 0.5) + lower * (1.5 - delta)) / step;
167  } else {
168  const double lower = pulse_[nbelow] - pulse_[nbelow - 1U];
169  const double upper = pulse_[nabove] - pulse_[nbelow];
170  return (lower * (0.5 - delta) + upper * (0.5 + delta)) / step;
171  }
172  }
173 
174  inline double secondDerivative(const double t) const {
175  const volatile double tmax = tmin_ + width_;
176  if (t < tmin_ || t > tmax || length_ < 3U)
177  return 0.0;
178  const unsigned lm1 = length_ - 1U;
179  const double step = width_ / lm1;
180  const double stepSq = step * step;
181  const double nSteps = (t - tmin_) / step;
182  unsigned nbelow = nSteps;
183  unsigned nabove = nbelow + 1;
184  if (nabove > lm1) {
185  nabove = lm1;
186  nbelow = nabove - 1U;
187  }
188 
189  if (nbelow == 0U) {
190  // The first interval
191  return (pulse_[2] - 2.0 * pulse_[1] + pulse_[0]) / stepSq;
192  } else if (nabove == lm1) {
193  // The last interval
194  return (pulse_[lm1] - 2.0 * pulse_[lm1 - 1U] + pulse_[lm1 - 2U]) / stepSq;
195  } else {
196  // One of the middle intervals
197  const double lower = pulse_[nbelow - 1U] - 2.0 * pulse_[nbelow] + pulse_[nabove];
198  const double upper = pulse_[nbelow] - 2.0 * pulse_[nabove] + pulse_[nabove + 1U];
199  const double delta = nSteps - nbelow;
200  return (lower * (1.0 - delta) + upper * delta) / stepSq;
201  }
202  }
203 
204  inline InterpolatedPulse& operator*=(const double scale) {
205  if (scale != 1.0) {
206  double* buf = &pulse_[0];
207  for (unsigned i = 0; i < length_; ++i)
208  *buf++ *= scale;
209  }
210  return *this;
211  }
212 
213  // Add another pulse to this one. Note that addition of another pulse
214  // will not change the start time or the width of _this_ pulse. The
215  // added pulse will be truncated as needed.
216  template <unsigned Len2>
218  const double step = width_ / (length_ - 1U);
219  for (unsigned i = 0; i < length_; ++i)
220  pulse_[i] += r(tmin_ + i * step);
221  return *this;
222  }
223 
224  template <unsigned Len2>
225  inline bool operator==(const InterpolatedPulse<Len2>& r) const {
226  if (!(tmin_ == r.tmin_ && width_ == r.width_ && length_ == r.length_))
227  return false;
228  const double* buf = &pulse_[0];
229  const double* rbuf = &r.pulse_[0];
230  for (unsigned i = 0; i < length_; ++i)
231  if (*buf++ != *rbuf++)
232  return false;
233  return true;
234  }
235 
236  template <unsigned Len2>
237  inline bool operator!=(const InterpolatedPulse<Len2>& r) const {
238  return !(*this == r);
239  }
240 
241  // Simple trapezoidal integration
242  inline double getIntegral() const {
243  const double* buf = &pulse_[0];
244  long double sum = buf[0] / 2.0;
245  const unsigned nIntervals = length_ - 1U;
246  for (unsigned i = 1U; i < nIntervals; ++i)
247  sum += buf[i];
248  sum += buf[nIntervals] / 2.0;
249  return sum * width_ / nIntervals;
250  }
251 
252  inline void setIntegral(const double newValue) {
253  const double integ = this->getIntegral();
254  if (integ == 0.0)
255  throw cms::Exception("In InterpolatedPulse::setIntegral division by zero");
256  *this *= (newValue / integ);
257  }
258 
259  inline double getPeakValue() const {
260  const double* buf = &pulse_[0];
261  double peak = buf[0];
262  for (unsigned i = 1U; i < length_; ++i)
263  if (buf[i] > peak)
264  peak = buf[i];
265  return peak;
266  }
267 
268  inline void setPeakValue(const double newValue) {
269  const double peak = this->getPeakValue();
270  if (peak == 0.0)
271  throw cms::Exception("In InterpolatedPulse::setPeakValue: division by zero");
272  *this *= (newValue / peak);
273  }
274 
275 private:
276  double pulse_[MaxLen];
277  double tmin_;
278  double width_;
279  unsigned length_;
280 
281  friend class boost::serialization::access;
282 
283  template <class Archive>
284  inline void serialize(Archive& ar, unsigned /* version */) {
285  ar& tmin_& width_& length_;
286 
287  // In case we are reading, it may be useful to verify
288  // that the length is reasonable
289  if (length_ > MaxLen)
290  throw cms::Exception("In InterpolatedPulse::serialize: buffer is not long enough");
291 
292  for (unsigned i = 0; i < length_; ++i)
293  ar& pulse_[i];
294  }
295 };
296 
297 // boost serialization version number for this template
298 namespace boost {
299  namespace serialization {
300  template <unsigned MaxLen>
301  struct version<InterpolatedPulse<MaxLen> > {
302  BOOST_STATIC_CONSTANT(int, value = 1);
303  };
304  } // namespace serialization
305 } // namespace boost
306 
307 #endif // CondFormats_HcalObjects_InterpolatedPulse_h_
Definition: CLHEP.h:16
void setStartTime(const double newStartTime)
void serialize(Archive &ar, unsigned)
double pulse_[MaxLen]
bool operator==(const InterpolatedPulse< Len2 > &r) const
double getStartTime() const
InterpolatedPulse & operator+=(const InterpolatedPulse< Len2 > &r)
void setShape(const Real *values, const unsigned len)
assert(be >=bs)
double getIntegral() const
double getTimeStep() const
InterpolatedPulse(const double tmin, const double tmax, const unsigned len)
double getPulseWidth() const
InterpolatedPulse(const double tmin, const double tmax, const Real *values, const unsigned len)
InterpolatedPulse & operator=(const InterpolatedPulse< Len2 > &r)
Definition: value.py:1
const double * getPulse() const
double operator()(const double t) const
InterpolatedPulse(const InterpolatedPulse< Len2 > &r)
static const double tmax[3]
double Real
double secondDerivative(const double t) const
InterpolatedPulse & operator=(const InterpolatedPulse &r)
double getPeakValue() const
InterpolatedPulse(const InterpolatedPulse &r)
InterpolatedPulse & operator*=(const double scale)
void setPeakValue(const double newValue)
bool operator!=(const InterpolatedPulse< Len2 > &r) const
double getStopTime() const
step
Definition: StallMonitor.cc:83
double derivative(const double t) const
void setPulseWidth(const double newWidth)
void setIntegral(const double newValue)
unsigned getLength() const
InterpolatedPulse(const unsigned len)