Line data Source code
1 : /**
2 : * @file TAMakerHorizontalMuonAlgorithm.cpp
3 : *
4 : * This is part of the DUNE DAQ Application Framework, copyright 2021.
5 : * Licensing/copyright details are in the COPYING file that you should have
6 : * received with this code.
7 : */
8 :
9 : #include "triggeralgs/HorizontalMuon/TAMakerHorizontalMuonAlgorithm.hpp"
10 : #include "TRACE/trace.h"
11 : #define TRACE_NAME "TAMakerHorizontalMuonAlgorithm"
12 : #include <vector>
13 : #include <math.h>
14 : #include <vector>
15 :
16 : using namespace triggeralgs;
17 :
18 : using Logging::TLVL_DEBUG_ALL;
19 : using Logging::TLVL_DEBUG_MEDIUM;
20 :
21 : void
22 0 : TAMakerHorizontalMuonAlgorithm::process(const TriggerPrimitive& input_tp,
23 : std::vector<TriggerActivity>& output_ta)
24 : {
25 :
26 0 : uint16_t adjacency;
27 :
28 : // Add useful info about recived TPs here for FW and SW TPG guys.
29 0 : if (m_print_tp_info) {
30 0 : TLOG_DEBUG(TLVL_DEBUG_ALL) << "[TAM:HM] TP Start Time: " << input_tp.time_start
31 0 : << ", TP ADC Sum: " << input_tp.adc_integral
32 0 : << ", TP SOT: " << input_tp.samples_over_threshold << ", TP ADC Peak: " << input_tp.adc_peak
33 0 : << ", TP Offline Channel ID: " << input_tp.channel;
34 0 : TLOG_DEBUG(TLVL_DEBUG_ALL) << "[TAM:HM] Adjacency of current window is: " << check_adjacency();
35 : }
36 :
37 : // 0) FIRST TP =====================================================================
38 : // The first time process() is called, reset the window object.
39 0 : if (m_current_window.is_empty()) {
40 0 : m_current_window.reset(input_tp);
41 0 : return;
42 : }
43 :
44 : // If the difference between the current TP's start time and the start of the window
45 : // is less than the specified window size, add the TP to the window.
46 0 : if ((input_tp.time_start - m_current_window.time_start) < m_window_length) {
47 0 : m_current_window.add(input_tp);
48 : }
49 :
50 : // 1) ADC THRESHOLD EXCEEDED =======================================================
51 : // If the addition of the current TP to the window would make it longer specified
52 : // window length, don't add it but check whether the ADC integral if the existing
53 : // window is above the configured threshold. If it is, and we are triggering on ADC,
54 : // make a TA and start a fresh window with the current TP.
55 0 : else if (m_current_window.adc_integral > m_adc_threshold && m_trigger_on_adc) {
56 :
57 0 : auto ta = construct_ta();
58 :
59 0 : TLOG_DEBUG(TLVL_DEBUG_MEDIUM) << "[TA:HM]: Emitting ADC threshold trigger with " << m_current_window.adc_integral
60 0 : << " window ADC integral. ta.time_start=" << ta.time_start
61 0 : << " ta.time_end=" << ta.time_end;
62 :
63 0 : output_ta.push_back(ta);
64 0 : m_current_window.reset(input_tp);
65 0 : }
66 :
67 : // 2) MULTIPLICITY - N UNQIUE CHANNELS EXCEEDED =====================================
68 : // If the addition of the current TP to the window would make it longer than the
69 : // specified window length, don't add it but check whether the number of hit channels
70 : // in the existing window is above the specified threshold. If it is, and we are triggering
71 : // on channel multiplicity, make a TA and start a fresh window with the current TP.
72 0 : else if (m_current_window.n_channels_hit() > m_n_channels_threshold && m_trigger_on_n_channels) {
73 :
74 0 : TLOG_DEBUG(TLVL_DEBUG_MEDIUM) << "[TAM:HM] Emitting multiplicity trigger with "
75 0 : << m_current_window.n_channels_hit() << " unique channels hit.";
76 :
77 0 : output_ta.push_back(construct_ta());
78 0 : m_current_window.reset(input_tp);
79 : }
80 :
81 : // 3) ADJACENCY THRESHOLD EXCEEDED ==================================================
82 : // If the addition of the current TP to the window would make it longer than the
83 : // specified window length, don't add it but check whether the adjacency of the
84 : // current window exceeds the configured threshold. If it does, and we are triggering
85 : // on adjacency, then create a TA and reset the window with the new/current TP.
86 0 : else if ((adjacency = check_adjacency()) > m_adjacency_threshold && m_trigger_on_adjacency) {
87 :
88 : //for (auto tp : m_current_window.inputs){ dump_tp(tp); }
89 :
90 : // Check for a new maximum, display the largest seen adjacency in the log.
91 : // uint16_t adjacency = check_adjacency();
92 0 : if (adjacency > m_max_adjacency) {
93 0 : m_max_adjacency = adjacency;
94 : }
95 0 : TLOG_DEBUG(TLVL_DEBUG_MEDIUM) << "[TAM:HM] Emitting track and multiplicity TA with adjacency "
96 0 : << check_adjacency() << " and multiplicity " << m_current_window.n_channels_hit()
97 0 : << ". The ADC integral of this TA is " << m_current_window.adc_integral
98 0 : << " and the largest longest track seen so far is " << m_max_adjacency;
99 :
100 0 : output_ta.push_back(construct_ta());
101 0 : m_current_window.reset(input_tp);
102 : }
103 :
104 : // Temporary triggering logic for Adam's large SOT TPs. Trigger on very large SOT TPs.
105 0 : else if (m_trigger_on_sot && input_tp.samples_over_threshold > m_sot_threshold) {
106 :
107 : // If the incoming TP has a large time over threshold, we might have a cluster of
108 : // interesting physics activity surrounding it. Trigger on that.
109 0 : TLOG_DEBUG(TLVL_DEBUG_MEDIUM) << "[TAM:HM] Emitting a TA due to a TP with a very large samples over threshold: "
110 0 : << input_tp.samples_over_threshold << " ticks and offline channel: " << input_tp.channel
111 0 : << ", where the ADC integral of that TP is " << input_tp.adc_integral;
112 0 : output_ta.push_back(construct_ta());
113 0 : m_current_window.reset(input_tp);
114 0 : }
115 :
116 : // 4) Otherwise, slide the window along using the current TP.
117 : else {
118 0 : m_current_window.move(input_tp, m_window_length);
119 : }
120 :
121 : return;
122 : }
123 :
124 : void
125 0 : TAMakerHorizontalMuonAlgorithm::configure(const nlohmann::json& config)
126 : {
127 0 : TriggerActivityMaker::configure(config);
128 :
129 0 : if (config.is_object()) {
130 0 : if (config.contains("trigger_on_adc"))
131 0 : m_trigger_on_adc = config["trigger_on_adc"];
132 0 : if (config.contains("trigger_on_n_channels"))
133 0 : m_trigger_on_n_channels = config["trigger_on_n_channels"];
134 0 : if (config.contains("adc_threshold"))
135 0 : m_adc_threshold = config["adc_threshold"];
136 0 : if (config.contains("n_channels_threshold"))
137 0 : m_n_channels_threshold = config["n_channels_threshold"];
138 0 : if (config.contains("window_length"))
139 0 : m_window_length = config["window_length"];
140 0 : if (config.contains("trigger_on_adjacency"))
141 0 : m_trigger_on_adjacency = config["trigger_on_adjacency"];
142 0 : if (config.contains("adjacency_tolerance"))
143 0 : m_adj_tolerance = config["adjacency_tolerance"];
144 0 : if (config.contains("adjacency_threshold"))
145 0 : m_adjacency_threshold = config["adjacency_threshold"];
146 0 : if (config.contains("print_tp_info"))
147 0 : m_print_tp_info = config["print_tp_info"];
148 0 : if (config.contains("trigger_on_sot"))
149 0 : m_trigger_on_sot = config["trigger_on_sot"];
150 0 : if (config.contains("sot_threshold"))
151 0 : m_sot_threshold = config["sot_threshold"];
152 : }
153 :
154 0 : }
155 :
156 : TriggerActivity
157 0 : TAMakerHorizontalMuonAlgorithm::construct_ta() const
158 : {
159 :
160 0 : TriggerActivity ta;
161 :
162 0 : TriggerPrimitive last_tp = m_current_window.inputs.back();
163 :
164 0 : ta.time_start = last_tp.time_start;
165 0 : ta.time_end = last_tp.time_start + last_tp.samples_over_threshold * 32; // FIXME: Replace the hard-coded SOT to TOT scaling.
166 0 : ta.time_peak = last_tp.samples_to_peak * 32 + last_tp.time_start; // FIXME: Replace STP to `time_peak` conversion.
167 0 : ta.time_activity = ta.time_peak;
168 0 : ta.channel_start = last_tp.channel;
169 0 : ta.channel_end = last_tp.channel;
170 0 : ta.channel_peak = last_tp.channel;
171 0 : ta.adc_integral = m_current_window.adc_integral;
172 0 : ta.adc_peak = last_tp.adc_integral;
173 0 : ta.detid = last_tp.detid;
174 0 : ta.type = TriggerActivity::Type::kTPC;
175 0 : ta.algorithm = TriggerActivity::Algorithm::kHorizontalMuon;
176 0 : ta.inputs = m_current_window.inputs;
177 :
178 0 : for( const auto& tp : ta.inputs ) {
179 0 : ta.time_start = std::min(ta.time_start, tp.time_start);
180 0 : ta.time_end = std::max(ta.time_end, tp.time_start + tp.samples_over_threshold * 32); // FIXME: Replace the hard-coded SOT to TOT scaling.
181 0 : ta.channel_start = std::min(ta.channel_start, channel_t(tp.channel));
182 0 : ta.channel_end = std::max(ta.channel_end, channel_t(tp.channel));
183 0 : if (tp.adc_peak > ta.adc_peak) {
184 0 : ta.time_peak = tp.samples_to_peak * 32 + tp.time_start; // FIXME: Replace STP to `time_peak` conversion.
185 0 : ta.adc_peak = tp.adc_peak;
186 0 : ta.channel_peak = tp.channel;
187 : }
188 : }
189 :
190 0 : return ta;
191 0 : }
192 :
193 : uint16_t
194 0 : TAMakerHorizontalMuonAlgorithm::check_adjacency() const
195 : {
196 : // This function returns the adjacency value for the current window, where adjacency
197 : // is defined as the maximum number of consecutive wires containing hits. It accepts
198 : // a configurable tolerance paramter, which allows up to adj_tolerance missing hits
199 : // on adjacent wires before restarting the adjacency count. The maximum gap is 4 which
200 : // comes from tuning on December 2021 coldbox data, and June 2022 coldbox runs.
201 :
202 0 : uint16_t adj = 1; // Initialise adjacency, 1 for the first wire.
203 0 : uint16_t max = 0; // Maximum adjacency of window, which this function returns
204 0 : unsigned int channel = 0; // Current channel ID
205 0 : unsigned int next_channel = 0; // Next channel ID
206 0 : unsigned int next = 0; // The next position in the hit channels vector
207 0 : unsigned int tol_count = 0; // Tolerance count, should not pass adj_tolerance
208 :
209 : // Generate a channelID ordered list of hit channels for this window
210 0 : std::vector<int> chanList;
211 0 : for (auto tp : m_current_window.inputs) {
212 0 : chanList.push_back(channel_t(tp.channel));
213 : }
214 0 : std::sort(chanList.begin(), chanList.end());
215 :
216 : // ADAJACENCY LOGIC ====================================================================
217 : // =====================================================================================
218 : // Adjcancency Tolerance = Number of times prepared to skip missed hits before resetting
219 : // the adjacency count. This accounts for things like dead channels / missed TPs. The
220 : // maximum gap is 4 which comes from tuning on December 2021 coldbox data, and June 2022
221 : // coldbox runs.
222 0 : for (int i = 0; i < chanList.size(); ++i) {
223 :
224 0 : next = (i + 1) % chanList.size(); // Loops back when outside of channel list range
225 0 : channel = chanList.at(i);
226 0 : next_channel = chanList.at(next); // Next channel with a hit
227 :
228 : // End of vector condition.
229 0 : if (next_channel == 0) {
230 0 : next_channel = channel - 1;
231 : }
232 :
233 : // Skip same channel hits.
234 0 : if (next_channel == channel) {
235 0 : continue;
236 : }
237 :
238 : // If next hit is on next channel, increment the adjacency count.
239 0 : else if (next_channel == channel + 1) {
240 0 : ++adj;
241 : }
242 :
243 : // If next channel is not on the next hit, but the 'second next', increase adjacency
244 : // but also tally up with the tolerance counter.
245 0 : else if (((next_channel == channel + 2) || (next_channel == channel + 3) || (next_channel == channel + 4) ||
246 0 : (next_channel == channel + 5)) &&
247 0 : (tol_count < m_adj_tolerance)) {
248 0 : ++adj;
249 0 : for (int i = 0; i < next_channel - channel; ++i) {
250 0 : ++tol_count;
251 : }
252 : }
253 :
254 : // If next hit isn't within reach, end the adjacency count and check for a new max.
255 : // Reset variables for next iteration.
256 : else {
257 0 : if (adj > max) {
258 : max = adj;
259 : }
260 : adj = 1;
261 : tol_count = 0;
262 : }
263 : }
264 :
265 0 : return max;
266 0 : }
267 :
268 : // =====================================================================================
269 : // Functions below this line are for debugging purposes.
270 : // =====================================================================================
271 : void
272 0 : TAMakerHorizontalMuonAlgorithm::add_window_to_record(TPWindow window)
273 : {
274 0 : m_window_record.push_back(window);
275 0 : return;
276 : }
277 :
278 : // Function to dump the details of the TA window currently on record
279 : void
280 0 : TAMakerHorizontalMuonAlgorithm::dump_window_record()
281 : {
282 0 : std::ofstream outfile;
283 0 : outfile.open("window_record_tam.csv", std::ios_base::app);
284 :
285 0 : for (auto window : m_window_record) {
286 0 : outfile << window.time_start << ",";
287 0 : outfile << window.inputs.back().time_start << ",";
288 0 : outfile << window.inputs.back().time_start - window.time_start << ",";
289 0 : outfile << window.adc_integral << ",";
290 0 : outfile << window.n_channels_hit() << ","; // Number of unique channels with hits
291 0 : outfile << window.inputs.size() << ","; // Number of TPs in TPWindow
292 0 : outfile << window.inputs.back().channel << ","; // Last TP Channel ID
293 0 : outfile << window.inputs.front().channel << ","; // First TP Channel ID
294 0 : outfile << check_adjacency() << ","; // New adjacency value for the window
295 0 : outfile << check_sot() << std::endl; // Summed window SOT
296 0 : }
297 :
298 0 : outfile.close();
299 0 : m_window_record.clear();
300 :
301 0 : return;
302 0 : }
303 :
304 : // Function to add current TP details to a text file for testing and debugging.
305 : void
306 0 : TAMakerHorizontalMuonAlgorithm::dump_tp(TriggerPrimitive const& input_tp)
307 : {
308 0 : std::ofstream outfile;
309 0 : outfile.open("coldbox_tps.txt", std::ios_base::app);
310 :
311 : // Output relevant TP information to file
312 0 : outfile << input_tp.time_start << " ";
313 0 : outfile << input_tp.samples_over_threshold << " "; // 50MHz ticks
314 0 : outfile << input_tp.samples_to_peak << " ";
315 0 : outfile << input_tp.channel << " "; // Offline channel ID
316 0 : outfile << input_tp.adc_integral << " ";
317 0 : outfile << input_tp.adc_peak << " ";
318 0 : outfile << input_tp.detid << " "; // Det ID - Identifies detector element
319 0 : outfile.close();
320 :
321 0 : return;
322 0 : }
323 :
324 : int
325 0 : TAMakerHorizontalMuonAlgorithm::check_sot() const
326 : {
327 : // Here, we just want to sum up all the sot values for each TP within window,
328 : // and return this sot of the window.
329 0 : int window_sot = 0;
330 0 : for (auto tp : m_current_window.inputs) {
331 0 : window_sot += tp.samples_over_threshold;
332 : }
333 :
334 0 : return window_sot;
335 : }
336 :
337 : // Register algo in TA Factory
338 12 : REGISTER_TRIGGER_ACTIVITY_MAKER(TRACE_NAME, TAMakerHorizontalMuonAlgorithm)
|
