MTDShapeBase.cc

Go to the documentation of this file.

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "SimFastTiming/FastTimingCommon/interface/MTDShapeBase.h"
 
MTDShapeBase::~MTDShapeBase() {}
 
MTDShapeBase::MTDShapeBase()
    : qNSecPerBin_(1. / kNBinsPerNSec), indexOfMax_(0), timeOfMax_(0.), shape_(DVec(k1NSecBinsTotal, 0.0)) {}
 
std::array<float, 3> MTDShapeBase::timeAtThr(const float scale, const float threshold1, const float threshold2) const {
  std::array<float, 3> times_tmp = {{0., 0., 0.}};
 
  // --- Check if the pulse amplitude is greater than threshold 2
  if (shape_[indexOfMax_] * scale < threshold2)
    return times_tmp;
 
  // --- Find the times corresponding to thresholds 1 and 2 on the pulse leading edge
  //     NB: To speed up the search we loop only on the rising edge
  unsigned int index_LT1 = 0;
  unsigned int index_LT2 = 0;
 
  for (unsigned int i = 0; i < indexOfMax_; ++i) {
    float amplitude = shape_[i] * scale;
 
    if (amplitude > threshold1 && index_LT1 == 0)
      index_LT1 = i;
 
    if (amplitude > threshold2 && index_LT2 == 0) {
      index_LT2 = i;
      break;
    }
  }
 
  // --- Find the time corresponding to thresholds 1 on the pulse falling edge
  unsigned int index_FT1 = 0;
 
  for (unsigned int i = shape_.size() - 1; i > indexOfMax_; i--) {
    float amplitude = shape_[i] * scale;
 
    if (amplitude > threshold1 && index_FT1 == 0) {
      index_FT1 = i + 1;
      break;
    }
  }
 
  if (index_LT1 != 0)
    times_tmp[0] = linear_interpolation(threshold1,
                                        (index_LT1 - 1) * qNSecPerBin_,
                                        index_LT1 * qNSecPerBin_,
                                        shape_[index_LT1 - 1] * scale,
                                        shape_[index_LT1] * scale);
 
  if (index_LT2 != 0)
    times_tmp[1] = linear_interpolation(threshold2,
                                        (index_LT2 - 1) * qNSecPerBin_,
                                        index_LT2 * qNSecPerBin_,
                                        shape_[index_LT2 - 1] * scale,
                                        shape_[index_LT2] * scale);
 
  if (index_FT1 != 0)
    times_tmp[2] = linear_interpolation(threshold1,
                                        (index_FT1 - 1) * qNSecPerBin_,
                                        index_FT1 * qNSecPerBin_,
                                        shape_[index_FT1 - 1] * scale,
                                        shape_[index_FT1] * scale);
 
  return times_tmp;
}
 
unsigned int MTDShapeBase::indexOfMax() const { return indexOfMax_; }
 
double MTDShapeBase::timeOfMax() const { return timeOfMax_; }
 
void MTDShapeBase::buildMe() {
  // --- Fill the vector with the pulse shape
  fillShape(shape_);
 
  // --- Find the index of maximum
  for (unsigned int i = 0; i < shape_.size(); ++i) {
    if (shape_[indexOfMax_] < shape_[i])
      indexOfMax_ = i;
  }
 
  if (indexOfMax_ != 0)
    timeOfMax_ = indexOfMax_ * qNSecPerBin_;
}
 
unsigned int MTDShapeBase::timeIndex(double aTime) const {
  const unsigned int index = aTime * kNBinsPerNSec;
 
  const bool bad = (k1NSecBinsTotal <= index);
 
  if (bad)
    LogDebug("MTDShapeBase") << " MTD pulse shape requested for out of range time " << aTime;
 
  return (bad ? k1NSecBinsTotal : index);
}
 
double MTDShapeBase::operator()(double aTime) const {
  // return pulse amplitude for request time in ns
  const unsigned int index(timeIndex(aTime));
  return (k1NSecBinsTotal == index ? 0 : shape_[index]);
}
 
double MTDShapeBase::linear_interpolation(
    const double& y, const double& x1, const double& x2, const double& y1, const double& y2) const {
  if (x1 == x2)
    throw cms::Exception("BadValue") << " MTDShapeBase: Trying to interpolate two points with the same x coordinate!";
 
  double a = (y2 - y1) / (x2 - x1);
  double b = y1 - a * x1;
 
  return (y - b) / a;
}