Local/src/TotemTimingRecHitProducerAlgorithm.cc

0001 /****************************************************************************
0002  *
0003  * This is a part of CTPPS offline software.
0004  * Authors:
0005  *   Laurent Forthomme (laurent.forthomme@cern.ch)
0006  *   Nicola Minafra
0007  *   Christopher Misan (krzysztof.misan@cern.ch)
0008  *
0009  ****************************************************************************/
0010
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012
0013 #include "RecoPPS/Local/interface/TotemTimingRecHitProducerAlgorithm.h"
0014 #include "RecoPPS/Local/interface/TotemTimingConversions.h"
0015
0016 #include <numeric>
0017
0018 //----------------------------------------------------------------------------------------------------
0019
0020 TotemTimingRecHitProducerAlgorithm::TotemTimingRecHitProducerAlgorithm(const edm::ParameterSet& iConfig)
0021     : TimingRecHitProducerAlgorithm(iConfig),
0022       mergeTimePeaks_(iConfig.getParameter<bool>("mergeTimePeaks")),
0023       baselinePoints_(iConfig.getParameter<int>("baselinePoints")),
0024       saturationLimit_(iConfig.getParameter<double>("saturationLimit")),
0025       cfdFraction_(iConfig.getParameter<double>("cfdFraction")),
0026       smoothingPoints_(iConfig.getParameter<int>("smoothingPoints")),
0027       lowPassFrequency_(iConfig.getParameter<double>("lowPassFrequency")),
0028       hysteresis_(iConfig.getParameter<double>("hysteresis")),
0029       sampicOffset_(iConfig.getParameter<double>("sampicOffset")),
0030       sampicSamplingPeriodNs_(iConfig.getParameter<double>("sampicSamplingPeriodNs")) {}
0031
0032 //----------------------------------------------------------------------------------------------------
0033
0034 void TotemTimingRecHitProducerAlgorithm::setCalibration(const PPSTimingCalibration& calib) {
0035   sampicConversions_ = std::make_unique<TotemTimingConversions>(sampicSamplingPeriodNs_, mergeTimePeaks_, calib);
0036 }
0037
0038 //----------------------------------------------------------------------------------------------------
0039
0040 void TotemTimingRecHitProducerAlgorithm::build(const CTPPSGeometry& geom,
0041                                                const edm::DetSetVector<TotemTimingDigi>& input,
0042                                                edm::DetSetVector<TotemTimingRecHit>& output) {
0043   for (const auto& vec : input) {
0044     const CTPPSDetId detid(vec.detId());
0045
0046     float x_pos = 0.f, y_pos = 0.f, z_pos = 0.f;
0047     float x_width = 0.f, y_width = 0.f, z_width = 0.f;
0048
0049     // retrieve the geometry element associated to this DetID ( if present )
0050     const DetGeomDesc* det = geom.sensorNoThrow(detid);
0051
0052     if (det) {
0053       x_pos = det->translation().x();
0054       y_pos = det->translation().y();
0055       z_pos = det->parentZPosition();  // retrieve the plane position;
0056
0057       const auto& diamondDimensions = det->getDiamondDimensions();
0058       x_width = 2.0 * diamondDimensions.xHalfWidth;  // parameters stand for half the size
0059       y_width = 2.0 * diamondDimensions.yHalfWidth;
0060       z_width = 2.0 * diamondDimensions.zHalfWidth;
0061     } else
0062       throw cms::Exception("TotemTimingRecHitProducerAlgorithm") << "Failed to retrieve a sensor for " << detid;
0063
0064     if (!sampicConversions_)
0065       throw cms::Exception("TotemTimingRecHitProducerAlgorithm") << "No timing conversion retrieved.";
0066
0067     edm::DetSet<TotemTimingRecHit>& rec_hits = output.find_or_insert(detid);
0068
0069     for (const auto& digi : vec) {
0070       const float triggerCellTimeInstant(sampicConversions_->triggerTime(digi));
0071       const float timePrecision(sampicConversions_->timePrecision(digi));
0072       const std::vector<float> time(sampicConversions_->timeSamples(digi));
0073       std::vector<float> data(sampicConversions_->voltSamples(digi));
0074
0075       auto [min_it, max_it] = std::minmax_element(data.begin(), data.end());
0076
0077       if (det->name() == "CTPPS_Diamond_Segment") {
0078         //flip the signal
0079         for (unsigned int i = 0; i < data.size(); ++i)
0080           data[i] = -data[i] + sampicOffset_;
0081       }
0082       RegressionResults baselineRegression = simplifiedLinearRegression(time, data, 0, baselinePoints_);
0083
0084       // remove baseline
0085       std::vector<float> dataCorrected(data.size());
0086       for (unsigned int i = 0; i < data.size(); ++i)
0087         dataCorrected[i] = data[i] - (baselineRegression.q + baselineRegression.m * time[i]);
0088
0089       auto max_corrected_it = std::max_element(dataCorrected.begin(), dataCorrected.end());
0090
0091       float t = TotemTimingRecHit::NO_T_AVAILABLE;
0092
0093       //if det==diamond then if (*min_it) > saturationLimit_), calculate the t, otherwise t=-100
0094       //if det==ufsd then if (*max_it) < saturationLimit_), calculate the t, otherwise t=-100
0095       //the difference between ufsd and sampic conditions comes from the signal flip
0096       if ((det->name() == "CTPPS_Diamond_Segment" && (*min_it) > saturationLimit_) ||
0097           (det->name() == "CTPPS_UFSD_Segment" && (*max_it) < saturationLimit_))
0098         t = constantFractionDiscriminator(time, dataCorrected);
0099       mode_ = TotemTimingRecHit::CFD;
0100       rec_hits.emplace_back(x_pos,
0101                             x_width,
0102                             y_pos,
0103                             y_width,
0104                             z_pos,
0105                             z_width,
0106                             t,
0107                             triggerCellTimeInstant,
0108                             timePrecision,
0109                             *max_corrected_it,
0110                             baselineRegression.rms,
0111                             mode_);
0112     }
0113   }
0114 }
0115
0116 //----------------------------------------------------------------------------------------------------
0117
0118 TotemTimingRecHitProducerAlgorithm::RegressionResults TotemTimingRecHitProducerAlgorithm::simplifiedLinearRegression(
0119     const std::vector<float>& time,
0120     const std::vector<float>& data,
0121     const unsigned int start_at,
0122     const unsigned int points) const {
0123   RegressionResults results;
0124   if (time.size() != data.size())
0125     return results;
0126   if (start_at > time.size())
0127     return results;
0128   unsigned int stop_at = std::min((unsigned int)time.size(), start_at + points);
0129   unsigned int realPoints = stop_at - start_at;
0130   auto t_begin = std::next(time.begin(), start_at);
0131   auto t_end = std::next(time.begin(), stop_at);
0132   auto d_begin = std::next(data.begin(), start_at);
0133   auto d_end = std::next(data.begin(), stop_at);
0134
0135   const float sx = std::accumulate(t_begin, t_end, 0.);
0136   const float sxx = std::inner_product(t_begin, t_end, t_begin, 0.);  // sum(t**2)
0137   const float sy = std::accumulate(d_begin, d_end, 0.);
0138
0139   float sxy = 0.;
0140   for (unsigned int i = 0; i < realPoints; ++i)
0141     sxy += time[i] * data[i];
0142
0143   // y = mx + q
0144   results.m = (sxy * realPoints - sx * sy) / (sxx * realPoints - sx * sx);
0145   results.q = sy / realPoints - results.m * sx / realPoints;
0146
0147   float correctedSyy = 0.;
0148   for (unsigned int i = 0; i < realPoints; ++i)
0149     correctedSyy += pow(data[i] - (results.m * time[i] + results.q), 2);
0150   results.rms = sqrt(correctedSyy / realPoints);
0151
0152   return results;
0153 }
0154
0155 //----------------------------------------------------------------------------------------------------
0156
0157 int TotemTimingRecHitProducerAlgorithm::fastDiscriminator(const std::vector<float>& data, float threshold) const {
0158   int threholdCrossingIndex = -1;
0159   bool above = false;
0160   bool lockForHysteresis = false;
0161
0162   for (unsigned int i = 0; i < data.size(); ++i) {
0163     // Look for first edge
0164     if (!above && !lockForHysteresis && data[i] > threshold) {
0165       threholdCrossingIndex = i;
0166       above = true;
0167       lockForHysteresis = true;
0168     }
0169     if (above && lockForHysteresis)  // NOTE: not else if!, "above" can be set in
0170                                      // the previous if
0171     {
0172       // Lock until above threshold_+hysteresis
0173       if (lockForHysteresis && data[i] > threshold + hysteresis_)
0174         lockForHysteresis = false;
0175       // Ignore noise peaks
0176       if (lockForHysteresis && data[i] < threshold) {
0177         above = false;
0178         lockForHysteresis = false;
0179         threholdCrossingIndex = -1;  // assigned because of noise
0180       }
0181     }
0182   }
0183
0184   return threholdCrossingIndex;
0185 }
0186
0187 float TotemTimingRecHitProducerAlgorithm::constantFractionDiscriminator(const std::vector<float>& time,
0188                                                                         const std::vector<float>& data) {
0189   std::vector<float> dataProcessed(data);
0190   if (lowPassFrequency_ != 0)  // Smoothing
0191     for (unsigned short i = 0; i < data.size(); ++i)
0192       for (unsigned short j = -smoothingPoints_ / 2; j <= +smoothingPoints_ / 2; ++j)
0193         if ((i + j) >= 0 && (i + j) < data.size() && j != 0) {
0194           float x = SINC_COEFFICIENT * lowPassFrequency_ * j;
0195           dataProcessed[i] += data[i + j] * std::sin(x) / x;
0196         }
0197   auto max_it = std::max_element(dataProcessed.begin(), dataProcessed.end());
0198   float max = *max_it;
0199
0200   float threshold = cfdFraction_ * max;
0201   int indexOfThresholdCrossing = fastDiscriminator(dataProcessed, threshold);
0202
0203   //--- necessary sanity checks before proceeding with the extrapolation
0204   return (indexOfThresholdCrossing >= 1 && indexOfThresholdCrossing >= baselinePoints_ &&
0205           indexOfThresholdCrossing < (int)time.size())
0206              ? (time[indexOfThresholdCrossing - 1] - time[indexOfThresholdCrossing]) /
0207                        (dataProcessed[indexOfThresholdCrossing - 1] - dataProcessed[indexOfThresholdCrossing]) *
0208                        (threshold - dataProcessed[indexOfThresholdCrossing]) +
0209                    time[indexOfThresholdCrossing]
0210              : (float)TotemTimingRecHit::NO_T_AVAILABLE;
0211 }