forked from cms-sw/cmssw
/
TotemTimingRecHitProducerAlgorithm.cc
195 lines (159 loc) · 8.34 KB
/
TotemTimingRecHitProducerAlgorithm.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
/****************************************************************************
*
* This is a part of CTPPS offline software.
* Authors:
* Laurent Forthomme (laurent.forthomme@cern.ch)
* Nicola Minafra
*
****************************************************************************/
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "RecoCTPPS/TotemRPLocal/interface/TotemTimingRecHitProducerAlgorithm.h"
#include <numeric>
//----------------------------------------------------------------------------------------------------
TotemTimingRecHitProducerAlgorithm::TotemTimingRecHitProducerAlgorithm(const edm::ParameterSet& iConfig)
: mergeTimePeaks_(iConfig.getParameter<bool>("mergeTimePeaks")),
baselinePoints_(iConfig.getParameter<int>("baselinePoints")),
saturationLimit_(iConfig.getParameter<double>("saturationLimit")),
cfdFraction_(iConfig.getParameter<double>("cfdFraction")),
smoothingPoints_(iConfig.getParameter<int>("smoothingPoints")),
lowPassFrequency_(iConfig.getParameter<double>("lowPassFrequency")),
hysteresis_(iConfig.getParameter<double>("hysteresis")) {}
//----------------------------------------------------------------------------------------------------
void TotemTimingRecHitProducerAlgorithm::setCalibration(const PPSTimingCalibration& calib) {
sampicConversions_.reset(new TotemTimingConversions(mergeTimePeaks_, calib));
}
//----------------------------------------------------------------------------------------------------
void TotemTimingRecHitProducerAlgorithm::build(const CTPPSGeometry& geom,
const edm::DetSetVector<TotemTimingDigi>& input,
edm::DetSetVector<TotemTimingRecHit>& output) {
for (const auto& vec : input) {
const TotemTimingDetId detid(vec.detId());
float x_pos = 0.f, y_pos = 0.f, z_pos = 0.f;
float x_width = 0.f, y_width = 0.f, z_width = 0.f;
// retrieve the geometry element associated to this DetID ( if present )
const DetGeomDesc* det = geom.getSensorNoThrow(detid);
if (det) {
x_pos = det->translation().x();
y_pos = det->translation().y();
z_pos = det->parentZPosition(); // retrieve the plane position;
x_width = 2.0 * det->params()[0], // parameters stand for half the size
y_width = 2.0 * det->params()[1], z_width = 2.0 * det->params()[2];
} else
throw cms::Exception("TotemTimingRecHitProducerAlgorithm") << "Failed to retrieve a sensor for " << detid;
if (!sampicConversions_)
throw cms::Exception("TotemTimingRecHitProducerAlgorithm") << "No timing conversion retrieved.";
edm::DetSet<TotemTimingRecHit>& rec_hits = output.find_or_insert(detid);
for (const auto& digi : vec) {
const float triggerCellTimeInstant(sampicConversions_->triggerTime(digi));
const float timePrecision(sampicConversions_->timePrecision(digi));
const std::vector<float> time(sampicConversions_->timeSamples(digi));
std::vector<float> data(sampicConversions_->voltSamples(digi));
auto max_it = std::max_element(data.begin(), data.end());
RegressionResults baselineRegression = simplifiedLinearRegression(time, data, 0, baselinePoints_);
// remove baseline
std::vector<float> dataCorrected(data.size());
for (unsigned int i = 0; i < data.size(); ++i)
dataCorrected[i] = data[i] - (baselineRegression.q + baselineRegression.m * time[i]);
auto max_corrected_it = std::max_element(dataCorrected.begin(), dataCorrected.end());
float t = TotemTimingRecHit::NO_T_AVAILABLE;
if (*max_it < saturationLimit_)
t = constantFractionDiscriminator(time, dataCorrected);
mode_ = TotemTimingRecHit::CFD;
rec_hits.push_back(TotemTimingRecHit(x_pos,
x_width,
y_pos,
y_width,
z_pos,
z_width, // spatial information
t,
triggerCellTimeInstant,
timePrecision,
*max_corrected_it,
baselineRegression.rms,
mode_));
}
}
}
//----------------------------------------------------------------------------------------------------
TotemTimingRecHitProducerAlgorithm::RegressionResults TotemTimingRecHitProducerAlgorithm::simplifiedLinearRegression(
const std::vector<float>& time,
const std::vector<float>& data,
const unsigned int start_at,
const unsigned int points) const {
RegressionResults results;
if (time.size() != data.size())
return results;
if (start_at > time.size())
return results;
unsigned int stop_at = std::min((unsigned int)time.size(), start_at + points);
unsigned int realPoints = stop_at - start_at;
auto t_begin = std::next(time.begin(), start_at);
auto t_end = std::next(time.begin(), stop_at);
auto d_begin = std::next(data.begin(), start_at);
auto d_end = std::next(data.begin(), stop_at);
const float sx = std::accumulate(t_begin, t_end, 0.);
const float sxx = std::inner_product(t_begin, t_end, t_begin, 0.); // sum(t**2)
const float sy = std::accumulate(d_begin, d_end, 0.);
float sxy = 0.;
for (unsigned int i = 0; i < realPoints; ++i)
sxy += time[i] * data[i];
// y = mx + q
results.m = (sxy * realPoints - sx * sy) / (sxx * realPoints - sx * sx);
results.q = sy / realPoints - results.m * sx / realPoints;
float correctedSyy = 0.;
for (unsigned int i = 0; i < realPoints; ++i)
correctedSyy += pow(data[i] - (results.m * time[i] + results.q), 2);
results.rms = sqrt(correctedSyy / realPoints);
return results;
}
//----------------------------------------------------------------------------------------------------
int TotemTimingRecHitProducerAlgorithm::fastDiscriminator(const std::vector<float>& data, float threshold) const {
int threholdCrossingIndex = -1;
bool above = false;
bool lockForHysteresis = false;
for (unsigned int i = 0; i < data.size(); ++i) {
// Look for first edge
if (!above && !lockForHysteresis && data[i] > threshold) {
threholdCrossingIndex = i;
above = true;
lockForHysteresis = true;
}
if (above && lockForHysteresis) // NOTE: not else if!, "above" can be set in
// the previous if
{
// Lock until above threshold_+hysteresis
if (lockForHysteresis && data[i] > threshold + hysteresis_)
lockForHysteresis = false;
// Ignore noise peaks
if (lockForHysteresis && data[i] < threshold) {
above = false;
lockForHysteresis = false;
threholdCrossingIndex = -1; // assigned because of noise
}
}
}
return threholdCrossingIndex;
}
float TotemTimingRecHitProducerAlgorithm::constantFractionDiscriminator(const std::vector<float>& time,
const std::vector<float>& data) {
std::vector<float> dataProcessed(data);
if (lowPassFrequency_ != 0) // Smoothing
for (unsigned short i = 0; i < data.size(); ++i)
for (unsigned short j = -smoothingPoints_ / 2; j <= +smoothingPoints_ / 2; ++j)
if ((i + j) >= 0 && (i + j) < data.size() && j != 0) {
float x = SINC_COEFFICIENT * lowPassFrequency_ * j;
dataProcessed[i] += data[i + j] * std::sin(x) / x;
}
auto max_it = std::max_element(dataProcessed.begin(), dataProcessed.end());
float max = *max_it;
float threshold = cfdFraction_ * max;
int indexOfThresholdCrossing = fastDiscriminator(dataProcessed, threshold);
//--- necessary sanity checks before proceeding with the extrapolation
return (indexOfThresholdCrossing >= 1 && indexOfThresholdCrossing >= baselinePoints_ &&
indexOfThresholdCrossing < (int)time.size())
? (time[indexOfThresholdCrossing - 1] - time[indexOfThresholdCrossing]) /
(dataProcessed[indexOfThresholdCrossing - 1] - dataProcessed[indexOfThresholdCrossing]) *
(threshold - dataProcessed[indexOfThresholdCrossing]) +
time[indexOfThresholdCrossing]
: (float)TotemTimingRecHit::NO_T_AVAILABLE;
}