File indexing completed on 2024-04-06 12:27:40
0001
0002
0003
0004
0005
0006
0007
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
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();
0056
0057 const auto& diamondDimensions = det->getDiamondDimensions();
0058 x_width = 2.0 * diamondDimensions.xHalfWidth;
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
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
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
0094
0095
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.);
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
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
0164 if (!above && !lockForHysteresis && data[i] > threshold) {
0165 threholdCrossingIndex = i;
0166 above = true;
0167 lockForHysteresis = true;
0168 }
0169 if (above && lockForHysteresis)
0170
0171 {
0172
0173 if (lockForHysteresis && data[i] > threshold + hysteresis_)
0174 lockForHysteresis = false;
0175
0176 if (lockForHysteresis && data[i] < threshold) {
0177 above = false;
0178 lockForHysteresis = false;
0179 threholdCrossingIndex = -1;
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)
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
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 }