Line data Source code
1 : #include "triggeralgs/SWIFT/TAMakerSWIFTAlgorithm.hpp"
2 :
3 : #include "TRACE/trace.h"
4 : #define TRACE_NAME "TriggerActivityMakerSWIFTPlugin"
5 :
6 : #include <cassert>
7 : #include <iostream>
8 : #include <cmath>
9 :
10 : namespace triggeralgs {
11 :
12 : // configuration
13 0 : void TriggerActivityMakerSWIFT::configure(const nlohmann::json& config)
14 : {
15 : //window settings
16 0 : m_window_length = config.value("window_length", 32000); // in DTS ticks (32 * 1000 readout ticks @ 500ns/tick)
17 0 : m_inspect_energy_threshold_sadc = config.value("inspect_energy_threshold_sadc", 15000);
18 0 : m_accept_energy_threshold_sadc = config.value("accept_energy_threshold_sadc", 55000);
19 :
20 : //tp filtering settings
21 0 : m_min_adc_peak = config.value("min_adc_peak", 80); //ADC
22 0 : m_min_samples_over_threshold = config.value("min_samples_over_threshold", 256); // 8 * 32 DTS ticks while still using TPv1 FIXME
23 :
24 : //clustering
25 : //these depend on detector properties. default values for HD
26 0 : m_cm_per_tick = config.value("cm_per_tick", 0.016 * 0.16); // sampling rate [us/tick] * drift velocity [cm/us]
27 0 : m_wire_pitch = config.value("wire_pitch", 0.48); // cm
28 0 : m_db_min_samples = config.value("min_samples", 2); //min. number of TPs for valid cluster
29 0 : m_db_eps = config.value("epsilon", 2); //dbscan search radius in cm
30 0 : m_cluster_energy_cut_sadc = config.value("cluster_energy_cut_sadc", 22000); // min energy of dominant cluster eng. in window for acceptance
31 :
32 0 : assert(m_window_length > 0);
33 0 : }
34 :
35 : //TP refinement
36 0 : bool TriggerActivityMakerSWIFT::preprocess( const TriggerPrimitive& input_tp){
37 0 : if ((input_tp.adc_peak < m_min_adc_peak) && (input_tp.samples_over_threshold < m_min_samples_over_threshold )) {
38 0 : return false;
39 : }
40 : return true;
41 : }
42 :
43 : // Reset window state
44 0 : void TriggerActivityMakerSWIFT::reset_window_state(uint64_t new_window_start) {
45 0 : m_window_start = new_window_start;
46 0 : m_window_energy_sadc = 0;
47 0 : m_tp_count = 0;
48 0 : m_current_ta = TriggerActivity();
49 0 : m_current_ta.time_start = m_window_start;
50 0 : }
51 :
52 : //main function for binning TPs into fixed-size time windows
53 : //assumes TPs are strictly time-ordered
54 0 : void TriggerActivityMakerSWIFT::process(const TriggerPrimitive& input_tp, std::vector<TriggerActivity>& output_tas)
55 : {
56 :
57 : // Apply TP filtering
58 0 : if (!preprocess(input_tp)) {
59 : return;
60 : }
61 :
62 : // If TP is valid, determine which window this TP belongs to
63 0 : const uint64_t tp_window_start = (input_tp.time_start / m_window_length) * m_window_length;
64 :
65 : // Initialise on first TP
66 0 : if (!m_initialised) {
67 :
68 : //create TA instance once at run start
69 0 : m_current_ta = TriggerActivity();
70 :
71 : //reset window state
72 0 : reset_window_state(tp_window_start);
73 0 : m_initialised = true;
74 : }
75 :
76 : //if TP belongs to past window (already closed), reject it (not ideal)
77 0 : if (tp_window_start < m_window_start){
78 : return;
79 : }
80 :
81 : //if TP belongs to future window, close existing TA and start a new one at current time
82 0 : if (tp_window_start > m_window_start){
83 0 : close_window(output_tas);
84 0 : reset_window_state(tp_window_start);
85 : }
86 :
87 : //If we got here, the TP belongs to the current window
88 0 : m_current_ta.inputs.push_back(input_tp);
89 0 : m_window_energy_sadc += input_tp.adc_integral;
90 0 : ++m_tp_count;
91 : }
92 :
93 : //function which gets called when the window is closed.
94 : //main window categorisation & TA generation happens here
95 0 : void TriggerActivityMakerSWIFT::close_window(std::vector<TriggerActivity>& output_tas) {
96 :
97 0 : if (m_current_ta.inputs.empty()) return;
98 :
99 : //Prompt window categorisatoin : immidiate accept, inspect, reject based on local energy in window
100 0 : WindowDecision decision;
101 0 : if (m_window_energy_sadc >= m_accept_energy_threshold_sadc) decision = WindowDecision::kAccept;
102 0 : else if (m_window_energy_sadc >= m_inspect_energy_threshold_sadc) decision = WindowDecision::kInspect;
103 : else return; // Reject
104 :
105 : // cluster inspect cases
106 0 : if (decision == WindowDecision::kInspect) {
107 0 : const uint64_t max_cluster_energy = extract_dominant_cluster_energy(m_current_ta.inputs, m_db_eps, m_db_min_samples);
108 0 : if (max_cluster_energy <= m_cluster_energy_cut_sadc) return; // reject window if it didn't pass inspection
109 : }
110 :
111 : // Emit TA: should only reach this step if dealing with Accept, or Inspect windows that passed clustering
112 0 : set_ta_attributes();
113 0 : output_tas.push_back(m_current_ta);
114 : }
115 :
116 :
117 :
118 : //Clustering function: density-based clustering in t-z, where both coordinates are expressed in cm.
119 : //returns only max eng. for now, as that's the param. based on which decision is made.
120 : //FIXME eventually want this step to return nclusrers, mean cluster eng. etc.
121 0 : uint64_t TriggerActivityMakerSWIFT::extract_dominant_cluster_energy(const std::vector<TriggerPrimitive>& tps, float eps, int min_samples){
122 :
123 0 : struct Point {
124 : float z;
125 : float t;
126 : uint64_t adc;
127 : };
128 :
129 : //since time and channel are in different units, express TPs as points in unified coordinate system (z,t)
130 0 : std::vector<Point> points;
131 0 : points.reserve(tps.size());
132 0 : for (const auto& tp : tps) {
133 0 : points.push_back({tp.channel * m_wire_pitch, (tp.time_start - m_window_start) * m_cm_per_tick, tp.adc_integral});
134 : }
135 :
136 0 : const int N = points.size(); //number of elements to cluster
137 0 : int cluster_id = 0;
138 0 : std::vector<int> labels(N, -1); // initialise all labels to noise for now
139 0 : std::vector<uint8_t> visited(N, 0);
140 0 : float eps2 = eps * eps; //clustering radius
141 :
142 0 : for (int i = 0; i < N; ++i) {
143 0 : if (visited[i]) continue;
144 0 : visited[i] = 1;
145 :
146 0 : std::vector<int> neigh;
147 0 : const float zi = points[i].z;
148 0 : const float ti = points[i].t;
149 0 : for (int j = 0; j < N; ++j) {
150 0 : float dz = points[j].z - zi;
151 0 : float dt = points[j].t - ti;
152 0 : if (dz*dz + dt*dt <= eps2) neigh.push_back(j);
153 : }
154 :
155 0 : if (neigh.size() < static_cast<size_t>(min_samples)) {
156 0 : labels[i] = -1;
157 0 : continue;
158 : }
159 :
160 0 : labels[i] = cluster_id;
161 :
162 0 : size_t k = 0;
163 0 : while (k < neigh.size()) {
164 0 : int j = neigh[k];
165 :
166 0 : if (!visited[j]) {
167 0 : visited[j] = 1;
168 :
169 0 : std::vector<int> neigh2;
170 0 : const float zj = points[j].z;
171 0 : const float tj = points[j].t;
172 0 : for (int m = 0; m < N; ++m) {
173 0 : float dz = points[m].z - zj;
174 0 : float dt = points[m].t - tj;
175 0 : if (dz*dz + dt*dt <= eps2) neigh2.push_back(m);
176 : }
177 0 : if (neigh2.size() >= static_cast<size_t>(min_samples)) {
178 0 : for (int m : neigh2) {
179 0 : if (std::find(neigh.begin(), neigh.end(), m) == neigh.end())
180 0 : neigh.push_back(m);
181 : }
182 : }
183 0 : }
184 :
185 0 : if (labels[j] == -1) labels[j] = cluster_id;
186 0 : ++k;
187 : }
188 0 : ++cluster_id;
189 0 : }
190 :
191 : // once clusters are formed, calculate the energies
192 0 : std::vector<uint64_t> cluster_sums(cluster_id, 0);
193 0 : for (int i = 0; i < N; ++i) {
194 0 : if (labels[i] >= 0) cluster_sums[labels[i]] += points[i].adc;
195 : }
196 : //return dominant cluster energy in window
197 0 : if (cluster_sums.empty()) return 0;
198 0 : return *std::max_element(cluster_sums.begin(), cluster_sums.end());
199 0 : }
200 :
201 :
202 : //TA feature extraction - FIXME most of these fields are somewhat useless
203 0 : void TriggerActivityMakerSWIFT::set_ta_attributes()
204 : {
205 0 : m_current_ta.time_start = m_window_start;
206 0 : m_current_ta.time_end = m_window_start + m_window_length;
207 0 : m_current_ta.adc_integral = m_window_energy_sadc;
208 :
209 0 : const TriggerPrimitive& first_tp = m_current_ta.inputs.front();
210 0 : m_current_ta.detid = first_tp.detid;
211 0 : m_current_ta.type = TriggerActivity::Type::kTPC;
212 0 : m_current_ta.algorithm = TriggerActivity::Algorithm::kSWIFT;
213 :
214 :
215 0 : dunedaq::trgdataformats::channel_t min_ch = first_tp.channel;
216 0 : dunedaq::trgdataformats::channel_t max_ch = first_tp.channel;
217 :
218 : // Peak quantities
219 0 : m_current_ta.adc_peak = 0;
220 0 : for (const auto& tp : m_current_ta.inputs) {
221 0 : if (tp.channel < min_ch) min_ch = tp.channel;
222 0 : if (tp.channel > max_ch) max_ch = tp.channel;
223 :
224 0 : if (tp.adc_peak > m_current_ta.adc_peak) {
225 0 : m_current_ta.adc_peak = tp.adc_peak;
226 0 : m_current_ta.channel_peak = tp.channel;
227 0 : m_current_ta.time_peak = tp.time_start;
228 : }
229 : }
230 :
231 0 : m_current_ta.channel_start = min_ch;
232 0 : m_current_ta.channel_end = max_ch;
233 0 : m_current_ta.time_activity = m_current_ta.time_peak;
234 :
235 0 : }
236 :
237 :
238 0 : void TriggerActivityMakerSWIFT::flush(timestamp_t /*until*/, std::vector<TriggerActivity>& output_tas)
239 : {
240 0 : if (!m_initialised) return;
241 0 : close_window(output_tas);
242 : }
243 :
244 12 : REGISTER_TRIGGER_ACTIVITY_MAKER(TRACE_NAME, TriggerActivityMakerSWIFT)
245 :
246 : } // namespace triggeralgs
|
