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