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.),
      fallTime_(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_; }

float MTDShapeBase::maximum() const { return shape_[indexOfMax_]; }

float MTDShapeBase::fallTime() const { return fallTime_; }

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;
  }

  for (unsigned int i = indexOfMax_; i < shape_.size(); ++i) {
    if (shape_[indexOfMax_] * 0.01 > shape_[i]) {
      fallTime_ = i * qNSecPerBin_;
      break;
    }
  }

  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;
}