HcalAlgos/src/LowPassFilterTiming.cc

0001 #include <cmath>
0002 #include <cassert>
0003 #include "FWCore/Utilities/interface/Exception.h"
0004
0005 #include "CalibCalorimetry/HcalAlgos/interface/LowPassFilterTiming.h"
0006
0007 #define NPARAMETERS 6U
0008
0009 //
0010 // The formula below join_location - h is aleft*(x_in - join_location) + bleft.
0011 // The formula above join_location + h is aright*(x_in - join_location) + bright.
0012 // Smooth cubic interpolation is performed within +-h of join_location.
0013 //
0014 static double two_joined_lines(const double x_in,
0015                                const double aleft,
0016                                const double bleft,
0017                                const double aright,
0018                                const double bright,
0019                                const double join_location,
0020                                const double h) {
0021   assert(h >= 0.0);
0022   const double x = x_in - join_location;
0023   if (x <= -h)
0024     return aleft * x + bleft;
0025   else if (x >= h)
0026     return aright * x + bright;
0027   else {
0028     const double vleft = -h * aleft + bleft;
0029     const double vright = aright * h + bright;
0030     const double b = (aright - aleft) / 4.0 / h;
0031     const double d = (vright + vleft - 2 * b * h * h) / 2.0;
0032     const double a = (d + aright * h - b * h * h - vright) / (2.0 * h * h * h);
0033     const double c = -(3 * d + aright * h + b * h * h - 3 * vright) / (2.0 * h);
0034     return ((a * x + b) * x + c) * x + d;
0035   }
0036 }
0037
0038 unsigned LowPassFilterTiming::nParameters() const { return NPARAMETERS; }
0039
0040 //
0041 // The time constant decreases linearly in the space of log(V + Vbias),
0042 // from tauMin+tauDelta when V = 0 to tauMin when V = V1.
0043 // Around log(V1 + Vbias), the curve is joined by a third order polynomial.
0044 // The width of the join is dLogV on both sides of log(V1 + Vbias).
0045 //
0046 // Value "tauDelta" = 0 can be used to create a constant time filter.
0047 //
0048 double LowPassFilterTiming::operator()(const double v, const double* pars, const unsigned nParams) const {
0049   assert(nParams == NPARAMETERS);
0050   assert(pars);
0051   unsigned ipar = 0;
0052   const double logVbias = pars[ipar++];
0053   const double logTauMin = pars[ipar++];
0054
0055   // The middle of the join region. Not a log actually,
0056   // it is simple in the log space.
0057   const double logV0 = pars[ipar++];
0058
0059   // Log of the width of the join region
0060   const double logDelta = pars[ipar++];
0061
0062   // Log of the minus negative slope
0063   const double slopeLog = pars[ipar++];
0064
0065   // Log of the maximum delay time (cutoff)
0066   const double tauMax = pars[ipar++];
0067   assert(ipar == NPARAMETERS);
0068
0069   // What happens for large (in magnitude) negative voltage inputs?
0070   const double Vbias = exp(logVbias);
0071   const double shiftedV = v + Vbias;
0072   if (shiftedV <= 0.0)
0073     return tauMax;
0074
0075   const double lg = log(shiftedV);
0076   const double delta = exp(logDelta);
0077   const double tauMin = exp(logTauMin);
0078   double result = two_joined_lines(lg, -exp(slopeLog), tauMin, 0.0, tauMin, logV0, delta);
0079   if (result > tauMax)
0080     result = tauMax;
0081   return result;
0082 }