51 {
52
53
54
55
56 std::map<int, const ldmx::FittedSiStripHit*> channel_map;
57 for (const auto& h : hits) {
58 const int ch = h.getStripID();
59 auto it = channel_map.find(ch);
60 if (it == channel_map.end()) {
61 channel_map[ch] = &h;
62 } else {
63
64 if (std::abs(h.getT0() - mean_time_ns_) <
65 std::abs(it->second->getT0() - mean_time_ns_)) {
66 it->second = &h;
67 }
68 }
69 }
70
71
72
73
74
75 const double snr_neighbor = neighbor_threshold_ * noise_sigma_adc_;
76 const double snr_seed = seed_threshold_ * noise_sigma_adc_;
77
78 std::set<int> clusterable_set;
79 std::vector<int> seed_channels;
80
81 for (const auto& [ch, hp] : channel_map) {
82 const double amp = hp->getAmplitude();
83 if (amp >= snr_neighbor) {
84 clusterable_set.insert(ch);
85 }
86 if (amp >= snr_seed && passesSeedCuts(*hp)) {
87 seed_channels.push_back(ch);
88 }
89 }
90
91
92 std::sort(seed_channels.begin(), seed_channels.end(), [&](int a, int b) {
93 return channel_map.at(a)->getAmplitude() >
94 channel_map.at(b)->getAmplitude();
95 });
96
97
98
99
100 std::vector<ClusterCandidate> clusters;
101
102 for (int seed_ch : seed_channels) {
103
104 if (clusterable_set.find(seed_ch) == clusterable_set.end()) continue;
105
106 ClusterCandidate cand;
107 double cluster_weighted_t = 0.0;
108 double cluster_total_amp = 0.0;
109 double cluster_noise_sq = 0.0;
110
111 std::deque<int> unchecked;
112 unchecked.push_back(seed_ch);
113 clusterable_set.erase(seed_ch);
114
115 while (!unchecked.empty()) {
116 const int cur_ch = unchecked.front();
117 unchecked.pop_front();
118
121
122
123 cand.strip_ids.push_back(cur_ch);
124 cluster_total_amp += amp;
125 cluster_weighted_t += amp * hit.
getT0();
126 cluster_noise_sq += noise_sigma_adc_ * noise_sigma_adc_;
127
128
129 for (int delta : {-1, +1}) {
130 const int nb_ch = cur_ch + delta;
131 if (clusterable_set.find(nb_ch) == clusterable_set.end()) continue;
132
133
134 if (!passesNeighborCuts(*channel_map.at(nb_ch), cluster_weighted_t,
135 cluster_total_amp)) {
136 continue;
137 }
138
139 unchecked.push_back(nb_ch);
140 clusterable_set.erase(nb_ch);
141 }
142 }
143
144
145 if (cluster_noise_sq <= 0.0) continue;
146 if (cluster_total_amp / std::sqrt(cluster_noise_sq) < cluster_threshold_) {
147 continue;
148 }
149
150
151
152
153 double sum_amp_strip = 0.0;
154 for (int ch : cand.strip_ids) {
155 sum_amp_strip += channel_map.at(ch)->getAmplitude() * ch;
156 }
157
158 cand.centroid_strip = sum_amp_strip / cluster_total_amp;
159 cand.total_amplitude = cluster_total_amp;
160 cand.time_ns = cluster_weighted_t / cluster_total_amp;
161 cand.n_strips = static_cast<int>(cand.strip_ids.size());
162
163
164
165
166
167
168 double sum_amp_dsq = 0.0;
169 for (int ch : cand.strip_ids) {
170 double d = ch - cand.centroid_strip;
171 sum_amp_dsq += channel_map.at(ch)->getAmplitude() * d * d;
172 }
173 constexpr double k_single_strip_sigma = 1.0 / 3.4641;
174 const double rms = std::sqrt(sum_amp_dsq / cluster_total_amp);
175 cand.sigma_strip = (rms > 0.0) ? rms : k_single_strip_sigma;
176 cand.layer_id = channel_map.at(seed_ch)->getLayerID();
177
178 clusters.push_back(std::move(cand));
179 }
180
181 return clusters;
182}
Result of fitting a pulse shape to the ADC samples of a single readout strip.
float getAmplitude() const
Fitted pedestal-subtracted peak amplitude [ADC counts].
float getT0() const
Fitted hit arrival time [ns] in the sample-window reference frame.