File indexing completed on 2023-03-17 11:24:07
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
0013 if (shape_[indexOfMax_] * scale < threshold2)
0014 return times_tmp;
0015
0016
0017
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
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
0075 fillShape(shape_);
0076
0077
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
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 }