From 2f52c5690a557e041ef1c07d11785a5178a8f65e Mon Sep 17 00:00:00 2001 From: kwawrows Date: Mon, 13 Apr 2026 07:55:24 +0200 Subject: [PATCH 1/3] Adding SWIFT TAMaker --- CMakeLists.txt | 1 + .../SWIFT/TAMakerSWIFTAlgorithm.hpp | 68 +++++ src/TAMakerSWIFTAlgorithm.cpp | 263 ++++++++++++++++++ 3 files changed, 332 insertions(+) create mode 100644 include/triggeralgs/SWIFT/TAMakerSWIFTAlgorithm.hpp create mode 100644 src/TAMakerSWIFTAlgorithm.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 875df3a..fc94fdd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -51,6 +51,7 @@ add_library(triggeralgs SHARED src/TCMakerDBSCANAlgorithm.cpp src/TAMakerChannelDistanceAlgorithm.cpp src/TCMakerChannelDistanceAlgorithm.cpp + src/TAMakerSWIFTAlgorithm.cpp src/TAMakerBundleNAlgorithm.cpp src/TCMakerBundleNAlgorithm.cpp src/TAMakerChannelAdjacencyAlgorithm.cpp diff --git a/include/triggeralgs/SWIFT/TAMakerSWIFTAlgorithm.hpp b/include/triggeralgs/SWIFT/TAMakerSWIFTAlgorithm.hpp new file mode 100644 index 0000000..b7d4c48 --- /dev/null +++ b/include/triggeralgs/SWIFT/TAMakerSWIFTAlgorithm.hpp @@ -0,0 +1,68 @@ +/* + * @File TAMakerSWIFTAlgorithm.hpp + * + * This is part of the DUNE DAQ Application Framework, copyright 2021. + * Licensing/copyright details are in the COPYING file that you should have + * received with this code. + */ + +#ifndef TRIGGERALGS_SWIFT_TRIGGERACTIVITYMAKERSWIFT_HPP_ +#define TRIGGERALGS_SWIFT_TRIGGERACTIVITYMAKERSWIFT_HPP_ + +#include "triggeralgs/TriggerActivityFactory.hpp" +#include + +namespace triggeralgs { + + class TriggerActivityMakerSWIFT : public TriggerActivityMaker { + public: + //void operator()(const TriggerPrimitive& input_tp, std::vector& output_tas); + void process(const TriggerPrimitive& input_tp, std::vector& output_tas) override; + void configure(const nlohmann::json& config); + void set_ta_attributes(); + void flush(timestamp_t until, std::vector& output_tas) override; + bool preprocess( const TriggerPrimitive& input_tp); + + private: + TriggerActivity m_current_ta; + + // -- Algortihm configuration -- + // static variables which get set once at run start & kept fixed afterwards + + //Window settings + uint64_t m_window_length = 32000; // fixed window size in DTS ticks (32 * 1k readout ticks @ 500ns/tick) + uint64_t m_inspect_energy_threshold = 15000; + uint64_t m_accept_energy_threshold = 55000; + + //TP pre-processing + uint16_t m_min_adc_peak = 80; + uint16_t m_min_samples_over_threshold = 8; + + //Clustering settings + //these depend on detector properties. default values for HD + float m_cm_per_tick = 0.016 * 0.16; // sampling rate [us/tick] * drift velocity [cm/us] + float m_wire_pitch = 0.48; //cm + int m_db_min_samples = 2; //min samples for valid cluster + float m_db_eps = 2; //nearest neighbour search radius in cm + uint64_t m_cluster_energy_cut = 22000; // min energy of dominant cluster for window to be accepted + + + // Dynamic vaiables which get updated during running to keep track of alg state. These are non-configurable + + //window state + bool m_initialised = false; + uint64_t m_window_start = 0; + uint64_t m_window_energy = 0; + uint16_t m_tp_count = 0; + + + //Additional SWIFT-specific functions and classes + enum class WindowDecision { Reject, Inspect, Accept }; + void close_window(std::vector& output_tas); + void reset_window_state(uint64_t new_window_start); + uint64_t extract_dominant_cluster_energy(const std::vector& tps, float eps, int min_samples); +}; + +} // namespace triggeralgs + +#endif // TRIGGERALGS_SWIFT_TRIGGERACTIVITYMAKERSWIFT_HPP_ diff --git a/src/TAMakerSWIFTAlgorithm.cpp b/src/TAMakerSWIFTAlgorithm.cpp new file mode 100644 index 0000000..e5c267b --- /dev/null +++ b/src/TAMakerSWIFTAlgorithm.cpp @@ -0,0 +1,263 @@ +#include "triggeralgs/SWIFT/TAMakerSWIFTAlgorithm.hpp" + +#include "TRACE/trace.h" +#define TRACE_NAME "TriggerActivityMakerSWIFTPlugin" + +#include +#include +#include + +namespace triggeralgs { + + // configuration + void TriggerActivityMakerSWIFT::configure(const nlohmann::json& config) + { + //window settings + if (config.contains("window_length")) + m_window_length = config["window_length"]; + + if (config.contains("inspect_energy_threshold")) + m_inspect_energy_threshold = config["inspect_energy_threshold"]; + + if (config.contains("accept_energy_threshold")) + m_accept_energy_threshold = config["accept_energy_threshold"]; + + //tp filtering settings + if (config.contains("min_adc_peak")) + m_min_adc_peak = config["min_adc_peak"]; + + if (config.contains("min_samples_over_threshold")) + m_min_samples_over_threshold = config["min_samples_over_threshold"]; + + //clustering + if (config.contains("cm_per_tick")) + m_cm_per_tick = config["cm_per_tick"]; + + if (config.contains("wire_pitch")) + m_wire_pitch = config["wire_pitch"]; + + if (config.contains("min_samples")) + m_db_min_samples = config["min_samples"]; + + if (config.contains("epsilon")) // NN search radius + m_db_eps = config["epsilon"]; + + if (config.contains("cluster_energy_cut")) + m_cluster_energy_cut = config["cluster_energy_cut"]; + + assert(m_window_length > 0); + } + + //TP refinement + bool TriggerActivityMakerSWIFT::preprocess( const TriggerPrimitive& input_tp){ + if ((input_tp.adc_peak < m_min_adc_peak) && (input_tp.samples_over_threshold < m_min_samples_over_threshold )) { + return false; + } + return true; + } + + // Reset window state + void TriggerActivityMakerSWIFT::reset_window_state(uint64_t new_window_start) { + m_window_start = new_window_start; + m_window_energy = 0; + m_tp_count = 0; + m_current_ta = TriggerActivity(); //.inputs.clear(); // reuse existing TA heap allocation + m_current_ta.time_start = m_window_start; + } + + //main function for binning TPs into fixed-size time windows + //assumes TPs are strictly time-ordered + void TriggerActivityMakerSWIFT::process(const TriggerPrimitive& input_tp, std::vector& output_tas) + { + + // Apply TP filtering + if (!preprocess(input_tp)) { + return; + } + + // If TP is valid, determine which window this TP belongs to + const uint64_t tp_window_start = (input_tp.time_start / m_window_length) * m_window_length; + + // Initialise on first TP + if (!m_initialised) { + + //create TA instance once at run start + m_current_ta = TriggerActivity(); + + //reset window state + reset_window_state(tp_window_start); + m_initialised = true; + } + + //if TP is tardy and belongs to past window - reject it for now + if (tp_window_start < m_window_start){ + return; + } + + //if TP belongs to future window, close current and jump ahead + if (tp_window_start > m_window_start){ + close_window(output_tas); + reset_window_state(tp_window_start); + } + + //If we got here, the TP belongs to the current window + m_current_ta.inputs.push_back(input_tp); + m_window_energy += input_tp.adc_integral; + ++m_tp_count; + } + + //function which gets called when the window is closed. + //main window categorisation & TA generation happens here + void TriggerActivityMakerSWIFT::close_window(std::vector& output_tas) { + + if (m_current_ta.inputs.empty()) return; + + //Prompt window categorisatoin : immidiate accept, inspect, reject based on local energy in window + WindowDecision decision; + if (m_window_energy >= m_accept_energy_threshold) decision = WindowDecision::Accept; + else if (m_window_energy >= m_inspect_energy_threshold) decision = WindowDecision::Inspect; + else return; // Reject + + // cluster inspect cases + if (decision == WindowDecision::Inspect) { + const uint64_t max_cluster_energy = extract_dominant_cluster_energy(m_current_ta.inputs, m_db_eps, m_db_min_samples); + if (max_cluster_energy <= m_cluster_energy_cut) return; // reject window if it didn't pass inspection + //for the time being, not updating the flag to keep track of what went through the inspect-->accept pipeline. FIXME? + } + + // Emit TA: should only reach this step if dealing with Accept, or Inspect windows that passed clustering + set_ta_attributes(); + output_tas.push_back(m_current_ta); + } + + + + //Clustering function: density-based clustering in t-z, where both coordinates are expressed in cm. + //returns only max eng. for now, as that's the param. based on which decision is made. + //FIXME eventually want this step to return nclusrers, mean cluster eng. etc. + uint64_t TriggerActivityMakerSWIFT::extract_dominant_cluster_energy(const std::vector& tps, float eps, int min_samples){ + + struct Point { + float z; + float t; + uint64_t adc; + }; + + //since time and channel are in different units, express TPs as points in unified coordinate system (z,t) + std::vector points; + points.reserve(tps.size()); + for (const auto& tp : tps) { + points.push_back({tp.channel * m_wire_pitch, (tp.time_start - m_window_start) * m_cm_per_tick, tp.adc_integral}); + } + + const int N = points.size(); //number of elements to cluster + int cluster_id = 0; + std::vector labels(N, -1); // initialise all labels to noise for now + std::vector visited(N, 0); + float eps2 = eps * eps; //clustering radius + + for (int i = 0; i < N; ++i) { + if (visited[i]) continue; + visited[i] = 1; + + std::vector neigh; + const float zi = points[i].z; + const float ti = points[i].t; + for (int j = 0; j < N; ++j) { + float dz = points[j].z - zi; + float dt = points[j].t - ti; + if (dz*dz + dt*dt <= eps2) neigh.push_back(j); + } + + if (neigh.size() < static_cast(min_samples)) { + labels[i] = -1; + continue; + } + + labels[i] = cluster_id; + + size_t k = 0; + while (k < neigh.size()) { + int j = neigh[k]; + + if (!visited[j]) { + visited[j] = 1; + + std::vector neigh2; + const float zj = points[j].z; + const float tj = points[j].t; + for (int m = 0; m < N; ++m) { + float dz = points[m].z - zj; + float dt = points[m].t - tj; + if (dz*dz + dt*dt <= eps2) neigh2.push_back(m); + } + if (neigh2.size() >= static_cast(min_samples)) { + for (int m : neigh2) { + if (std::find(neigh.begin(), neigh.end(), m) == neigh.end()) + neigh.push_back(m); + } + } + } + + if (labels[j] == -1) labels[j] = cluster_id; + ++k; + } + ++cluster_id; + } + + // once clusters are formed, calculate the energies + std::vector cluster_sums(cluster_id, 0); + for (int i = 0; i < N; ++i) { + if (labels[i] >= 0) cluster_sums[labels[i]] += points[i].adc; + } + //return dominant cluster energy in window + if (cluster_sums.empty()) return 0; + return *std::max_element(cluster_sums.begin(), cluster_sums.end()); + } + + + //TA feature extraction - FIXME most of these fields are somewhat useless + void TriggerActivityMakerSWIFT::set_ta_attributes() + { + m_current_ta.time_start = m_window_start; + m_current_ta.time_end = m_window_start + m_window_length; + m_current_ta.adc_integral = m_window_energy; + + const TriggerPrimitive& first_tp = m_current_ta.inputs.front(); + m_current_ta.detid = first_tp.detid; + m_current_ta.type = TriggerActivity::Type::kTPC; + m_current_ta.algorithm = TriggerActivity::Algorithm::kUnknown; + + + dunedaq::trgdataformats::channel_t min_ch = first_tp.channel; + dunedaq::trgdataformats::channel_t max_ch = first_tp.channel; + + // Peak quantities + m_current_ta.adc_peak = 0; + for (const auto& tp : m_current_ta.inputs) { + if (tp.channel < min_ch) min_ch = tp.channel; + if (tp.channel > max_ch) max_ch = tp.channel; + + if (tp.adc_peak > m_current_ta.adc_peak) { + m_current_ta.adc_peak = tp.adc_peak; + m_current_ta.channel_peak = tp.channel; + m_current_ta.time_peak = tp.time_start; + } + } + + m_current_ta.channel_start = min_ch; + m_current_ta.channel_end = max_ch; + m_current_ta.time_activity = m_current_ta.time_peak; + + } + + + void TriggerActivityMakerSWIFT::flush(timestamp_t /*until*/, std::vector& output_tas) + { + if (!m_initialised) return; + close_window(output_tas); + } + + REGISTER_TRIGGER_ACTIVITY_MAKER(TRACE_NAME, TriggerActivityMakerSWIFT) + +} // namespace triggeralgs From 99e51ea0d0af6ca7d3a65fe9af071e1257c7978b Mon Sep 17 00:00:00 2001 From: kwawrows Date: Fri, 17 Apr 2026 13:28:03 +0200 Subject: [PATCH 2/3] Improving configuration and updating alg. variable names --- .../SWIFT/TAMakerSWIFTAlgorithm.hpp | 28 ++++---- src/TAMakerSWIFTAlgorithm.cpp | 65 +++++++------------ 2 files changed, 37 insertions(+), 56 deletions(-) diff --git a/include/triggeralgs/SWIFT/TAMakerSWIFTAlgorithm.hpp b/include/triggeralgs/SWIFT/TAMakerSWIFTAlgorithm.hpp index b7d4c48..6874212 100644 --- a/include/triggeralgs/SWIFT/TAMakerSWIFTAlgorithm.hpp +++ b/include/triggeralgs/SWIFT/TAMakerSWIFTAlgorithm.hpp @@ -30,21 +30,20 @@ namespace triggeralgs { // static variables which get set once at run start & kept fixed afterwards //Window settings - uint64_t m_window_length = 32000; // fixed window size in DTS ticks (32 * 1k readout ticks @ 500ns/tick) - uint64_t m_inspect_energy_threshold = 15000; - uint64_t m_accept_energy_threshold = 55000; + uint64_t m_window_length; + uint64_t m_inspect_energy_threshold_sadc; + uint64_t m_accept_energy_threshold_sadc; //TP pre-processing - uint16_t m_min_adc_peak = 80; - uint16_t m_min_samples_over_threshold = 8; + uint16_t m_min_adc_peak; + uint16_t m_min_samples_over_threshold; //Clustering settings - //these depend on detector properties. default values for HD - float m_cm_per_tick = 0.016 * 0.16; // sampling rate [us/tick] * drift velocity [cm/us] - float m_wire_pitch = 0.48; //cm - int m_db_min_samples = 2; //min samples for valid cluster - float m_db_eps = 2; //nearest neighbour search radius in cm - uint64_t m_cluster_energy_cut = 22000; // min energy of dominant cluster for window to be accepted + float m_cm_per_tick; + float m_wire_pitch; + int m_db_min_samples; + float m_db_eps; + uint64_t m_cluster_energy_cut_sadc; // Dynamic vaiables which get updated during running to keep track of alg state. These are non-configurable @@ -52,16 +51,15 @@ namespace triggeralgs { //window state bool m_initialised = false; uint64_t m_window_start = 0; - uint64_t m_window_energy = 0; + uint64_t m_window_energy_sadc = 0; uint16_t m_tp_count = 0; - //Additional SWIFT-specific functions and classes - enum class WindowDecision { Reject, Inspect, Accept }; + enum class WindowDecision { kReject, kInspect, kAccept }; void close_window(std::vector& output_tas); void reset_window_state(uint64_t new_window_start); uint64_t extract_dominant_cluster_energy(const std::vector& tps, float eps, int min_samples); -}; + }; } // namespace triggeralgs diff --git a/src/TAMakerSWIFTAlgorithm.cpp b/src/TAMakerSWIFTAlgorithm.cpp index e5c267b..39a1992 100644 --- a/src/TAMakerSWIFTAlgorithm.cpp +++ b/src/TAMakerSWIFTAlgorithm.cpp @@ -13,37 +13,21 @@ namespace triggeralgs { void TriggerActivityMakerSWIFT::configure(const nlohmann::json& config) { //window settings - if (config.contains("window_length")) - m_window_length = config["window_length"]; - - if (config.contains("inspect_energy_threshold")) - m_inspect_energy_threshold = config["inspect_energy_threshold"]; - - if (config.contains("accept_energy_threshold")) - m_accept_energy_threshold = config["accept_energy_threshold"]; + m_window_length = config.value("window_length", 32000); // in DTS ticks (32 * 1000 readout ticks @ 500ns/tick) + m_inspect_energy_threshold_sadc = config.value("inspect_energy_threshold_sadc", 15000); + m_accept_energy_threshold_sadc = config.value("accept_energy_threshold_sadc", 55000); //tp filtering settings - if (config.contains("min_adc_peak")) - m_min_adc_peak = config["min_adc_peak"]; - - if (config.contains("min_samples_over_threshold")) - m_min_samples_over_threshold = config["min_samples_over_threshold"]; + m_min_adc_peak = config.value("min_adc_peak", 80); //ADC + m_min_samples_over_threshold = config.value("min_samples_over_threshold", 256); // 8 * 32 DTS ticks while still using TPv1 FIXME //clustering - if (config.contains("cm_per_tick")) - m_cm_per_tick = config["cm_per_tick"]; - - if (config.contains("wire_pitch")) - m_wire_pitch = config["wire_pitch"]; - - if (config.contains("min_samples")) - m_db_min_samples = config["min_samples"]; - - if (config.contains("epsilon")) // NN search radius - m_db_eps = config["epsilon"]; - - if (config.contains("cluster_energy_cut")) - m_cluster_energy_cut = config["cluster_energy_cut"]; + //these depend on detector properties. default values for HD + m_cm_per_tick = config.value("cm_per_tick", 0.016 * 0.16); // sampling rate [us/tick] * drift velocity [cm/us] + m_wire_pitch = config.value("wire_pitch", 0.48); // cm + m_db_min_samples = config.value("min_samples", 2); //min. number of TPs for valid cluster + m_db_eps = config.value("epsilon", 2); //dbscan search radius in cm + m_cluster_energy_cut_sadc = config.value("cluster_energy_cut_sadc", 22000); // min energy of dominant cluster eng. in window for acceptance assert(m_window_length > 0); } @@ -58,10 +42,10 @@ namespace triggeralgs { // Reset window state void TriggerActivityMakerSWIFT::reset_window_state(uint64_t new_window_start) { - m_window_start = new_window_start; - m_window_energy = 0; - m_tp_count = 0; - m_current_ta = TriggerActivity(); //.inputs.clear(); // reuse existing TA heap allocation + m_window_start = new_window_start; + m_window_energy_sadc = 0; + m_tp_count = 0; + m_current_ta = TriggerActivity(); m_current_ta.time_start = m_window_start; } @@ -89,12 +73,12 @@ namespace triggeralgs { m_initialised = true; } - //if TP is tardy and belongs to past window - reject it for now + //if TP belongs to past window (already closed), reject it (not ideal) if (tp_window_start < m_window_start){ return; } - //if TP belongs to future window, close current and jump ahead + //if TP belongs to future window, close existing TA and start a new one at current time if (tp_window_start > m_window_start){ close_window(output_tas); reset_window_state(tp_window_start); @@ -102,7 +86,7 @@ namespace triggeralgs { //If we got here, the TP belongs to the current window m_current_ta.inputs.push_back(input_tp); - m_window_energy += input_tp.adc_integral; + m_window_energy_sadc += input_tp.adc_integral; ++m_tp_count; } @@ -114,15 +98,14 @@ namespace triggeralgs { //Prompt window categorisatoin : immidiate accept, inspect, reject based on local energy in window WindowDecision decision; - if (m_window_energy >= m_accept_energy_threshold) decision = WindowDecision::Accept; - else if (m_window_energy >= m_inspect_energy_threshold) decision = WindowDecision::Inspect; + if (m_window_energy_sadc >= m_accept_energy_threshold_sadc) decision = WindowDecision::kAccept; + else if (m_window_energy_sadc >= m_inspect_energy_threshold_sadc) decision = WindowDecision::kInspect; else return; // Reject // cluster inspect cases - if (decision == WindowDecision::Inspect) { + if (decision == WindowDecision::kInspect) { const uint64_t max_cluster_energy = extract_dominant_cluster_energy(m_current_ta.inputs, m_db_eps, m_db_min_samples); - if (max_cluster_energy <= m_cluster_energy_cut) return; // reject window if it didn't pass inspection - //for the time being, not updating the flag to keep track of what went through the inspect-->accept pipeline. FIXME? + if (max_cluster_energy <= m_cluster_energy_cut_sadc) return; // reject window if it didn't pass inspection } // Emit TA: should only reach this step if dealing with Accept, or Inspect windows that passed clustering @@ -221,12 +204,12 @@ namespace triggeralgs { { m_current_ta.time_start = m_window_start; m_current_ta.time_end = m_window_start + m_window_length; - m_current_ta.adc_integral = m_window_energy; + m_current_ta.adc_integral = m_window_energy_sadc; const TriggerPrimitive& first_tp = m_current_ta.inputs.front(); m_current_ta.detid = first_tp.detid; m_current_ta.type = TriggerActivity::Type::kTPC; - m_current_ta.algorithm = TriggerActivity::Algorithm::kUnknown; + m_current_ta.algorithm = TriggerActivity::Algorithm::kUnknown; //FIXME dunedaq::trgdataformats::channel_t min_ch = first_tp.channel; From 4332320dbaa227ade215e4667ba4180c54ddd8b4 Mon Sep 17 00:00:00 2001 From: kwawrows Date: Thu, 30 Apr 2026 11:46:11 +0200 Subject: [PATCH 3/3] Replace default kUnknown with kSWIFT after alg. registration in trgdataformats --- src/TAMakerSWIFTAlgorithm.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TAMakerSWIFTAlgorithm.cpp b/src/TAMakerSWIFTAlgorithm.cpp index 39a1992..eaa358a 100644 --- a/src/TAMakerSWIFTAlgorithm.cpp +++ b/src/TAMakerSWIFTAlgorithm.cpp @@ -209,7 +209,7 @@ namespace triggeralgs { const TriggerPrimitive& first_tp = m_current_ta.inputs.front(); m_current_ta.detid = first_tp.detid; m_current_ta.type = TriggerActivity::Type::kTPC; - m_current_ta.algorithm = TriggerActivity::Algorithm::kUnknown; //FIXME + m_current_ta.algorithm = TriggerActivity::Algorithm::kSWIFT; dunedaq::trgdataformats::channel_t min_ch = first_tp.channel;