FastTimingCommon/src/MTDShapeBase.cc

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "SimFastTiming/FastTimingCommon/interface/MTDShapeBase.h"
0003
0004 MTDShapeBase::~MTDShapeBase() {}
0005
0006 MTDShapeBase::MTDShapeBase()
0007     : qNSecPerBin_(1. / kNBinsPerNSec), indexOfMax_(0), timeOfMax_(0.), shape_(DVec(k1NSecBinsTotal, 0.0)) {}
0008
0009 std::array<float, 3> MTDShapeBase::timeAtThr(const float scale, const float threshold1, const float threshold2) const {
0010   std::array<float, 3> times_tmp = {{0., 0., 0.}};
0011
0012   // --- Check if the pulse amplitude is greater than threshold 2
0013   if (shape_[indexOfMax_] * scale < threshold2)
0014     return times_tmp;
0015
0016   // --- Find the times corresponding to thresholds 1 and 2 on the pulse leading edge
0017   //     NB: To speed up the search we loop only on the rising edge
0018   unsigned int index_LT1 = 0;
0019   unsigned int index_LT2 = 0;
0020
0021   for (unsigned int i = 0; i < indexOfMax_; ++i) {
0022     float amplitude = shape_[i] * scale;
0023
0024     if (amplitude > threshold1 && index_LT1 == 0)
0025       index_LT1 = i;
0026
0027     if (amplitude > threshold2 && index_LT2 == 0) {
0028       index_LT2 = i;
0029       break;
0030     }
0031   }
0032
0033   // --- Find the time corresponding to thresholds 1 on the pulse falling edge
0034   unsigned int index_FT1 = 0;
0035
0036   for (unsigned int i = shape_.size() - 1; i > indexOfMax_; i--) {
0037     float amplitude = shape_[i] * scale;
0038
0039     if (amplitude > threshold1 && index_FT1 == 0) {
0040       index_FT1 = i + 1;
0041       break;
0042     }
0043   }
0044
0045   if (index_LT1 != 0)
0046     times_tmp[0] = linear_interpolation(threshold1,
0047                                         (index_LT1 - 1) * qNSecPerBin_,
0048                                         index_LT1 * qNSecPerBin_,
0049                                         shape_[index_LT1 - 1] * scale,
0050                                         shape_[index_LT1] * scale);
0051
0052   if (index_LT2 != 0)
0053     times_tmp[1] = linear_interpolation(threshold2,
0054                                         (index_LT2 - 1) * qNSecPerBin_,
0055                                         index_LT2 * qNSecPerBin_,
0056                                         shape_[index_LT2 - 1] * scale,
0057                                         shape_[index_LT2] * scale);
0058
0059   if (index_FT1 != 0)
0060     times_tmp[2] = linear_interpolation(threshold1,
0061                                         (index_FT1 - 1) * qNSecPerBin_,
0062                                         index_FT1 * qNSecPerBin_,
0063                                         shape_[index_FT1 - 1] * scale,
0064                                         shape_[index_FT1] * scale);
0065
0066   return times_tmp;
0067 }
0068
0069 unsigned int MTDShapeBase::indexOfMax() const { return indexOfMax_; }
0070
0071 double MTDShapeBase::timeOfMax() const { return timeOfMax_; }
0072
0073 void MTDShapeBase::buildMe() {
0074   // --- Fill the vector with the pulse shape
0075   fillShape(shape_);
0076
0077   // --- Find the index of maximum
0078   for (unsigned int i = 0; i < shape_.size(); ++i) {
0079     if (shape_[indexOfMax_] < shape_[i])
0080       indexOfMax_ = i;
0081   }
0082
0083   if (indexOfMax_ != 0)
0084     timeOfMax_ = indexOfMax_ * qNSecPerBin_;
0085 }
0086
0087 unsigned int MTDShapeBase::timeIndex(double aTime) const {
0088   const unsigned int index = aTime * kNBinsPerNSec;
0089
0090   const bool bad = (k1NSecBinsTotal <= index);
0091
0092   if (bad)
0093     LogDebug("MTDShapeBase") << " MTD pulse shape requested for out of range time " << aTime;
0094
0095   return (bad ? k1NSecBinsTotal : index);
0096 }
0097
0098 double MTDShapeBase::operator()(double aTime) const {
0099   // return pulse amplitude for request time in ns
0100   const unsigned int index(timeIndex(aTime));
0101   return (k1NSecBinsTotal == index ? 0 : shape_[index]);
0102 }
0103
0104 double MTDShapeBase::linear_interpolation(
0105     const double& y, const double& x1, const double& x2, const double& y1, const double& y2) const {
0106   if (x1 == x2)
0107     throw cms::Exception("BadValue") << " MTDShapeBase: Trying to interpolate two points with the same x coordinate!";
0108
0109   double a = (y2 - y1) / (x2 - x1);
0110   double b = y1 - a * x1;
0111
0112   return (y - b) / a;
0113 }