Back to home page

Project CMSSW displayed by LXR

 
 

    


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