Line data Source code
1 : /**
2 : * @file TAMakerProtoDUNEBSMWindowAlgorithm.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/ProtoDUNEBSMWindow/TAMakerProtoDUNEBSMWindowAlgorithm.hpp"
10 :
11 : #include "TRACE/trace.h"
12 : #define TRACE_NAME "TAMakerProtoDUNEBSMWindowAlgorithm"
13 :
14 : #include <vector>
15 : #include <chrono>
16 :
17 : using namespace triggeralgs;
18 : using Logging::TLVL_DEBUG_ALL;
19 : using Logging::TLVL_DEBUG_HIGH;
20 : using Logging::TLVL_DEBUG_LOW;
21 : using Logging::TLVL_IMPORTANT;
22 :
23 : void
24 0 : TAMakerProtoDUNEBSMWindowAlgorithm::process(const TriggerPrimitive& input_tp, std::vector<TriggerActivity>& output_ta)
25 : {
26 :
27 0 : if(m_current_window.is_empty()){
28 : // Reset window with new TP
29 0 : m_current_window.reset(input_tp);
30 : // Initialise last time an XGBoost prediction was made
31 0 : m_last_pred_time = input_tp.time_start;
32 : // Iterate number of TPs in the window
33 0 : m_primitive_count++;
34 : // First time operator is called set ROP first and last channel
35 0 : unsigned int detelement = channelMap->get_element_id_from_offline_channel(input_tp.channel);
36 0 : unsigned int plane = channelMap->get_plane_from_offline_channel(input_tp.channel);
37 :
38 : // Are we on the collection plane? Use XGBoost model for collection plane TPs
39 : // evaluate the sum of the TP charge if on induction planes
40 : // Induction plane IDs = 0, 1
41 : // Collection plane ID = 2
42 0 : if (plane > 1) m_collection_plane = true;
43 :
44 : // Use PlaneInfo object to get the first and last channels on plane
45 0 : PlaneInfo plane_info = m_det_plane_map.get_plane_info(m_channel_map_name, detelement, plane);
46 0 : m_first_channel = static_cast<channel_t>(plane_info.min_channel);
47 0 : m_n_channels_on_plane = static_cast<channel_t>(plane_info.n_channels);
48 :
49 : // If we are in PD-VD use 'effective' channel mapping for CRPs
50 : // (but only for collection plane)
51 0 : if (m_pdvd_map && m_collection_plane) {
52 0 : m_pdvd_eff_channel_mapper = std::make_unique<PDVDEffectiveChannelMap>(plane_info.min_channel, plane_info.n_channels);
53 : // Get the first effective channel and number of effective channels on the plane
54 0 : m_first_channel = m_pdvd_eff_channel_mapper->remapCollectionPlaneChannel(m_first_channel);
55 0 : m_n_channels_on_plane = m_pdvd_eff_channel_mapper->getNEffectiveChannels();
56 : }
57 :
58 0 : TLOG_DEBUG(TLVL_DEBUG_ALL) << "[TAM:BSMW] 1st Chan = " << m_first_channel << ", N channels on plane = " << m_n_channels_on_plane << std::endl
59 0 : << "Number of channel bins = " << m_num_chanbins;
60 0 : return;
61 : }
62 :
63 : // If the difference between the current TP's start time and the start of the window
64 : // is less than the specified window size, add the TP to the window.
65 0 : if((input_tp.time_start - m_current_window.time_start) < m_window_length){
66 0 : TLOG_DEBUG(TLVL_DEBUG_HIGH) << "[TAM:BSMW] Window not yet complete, adding the input_tp to the window.";
67 0 : m_current_window.add(input_tp);
68 : }
69 :
70 : // If the addition of the current TP to the window would make it longer
71 : // than the specified window length, don't add it
72 : // First, if these are not collection plane TPs, just evaluate the total charge
73 : // If the total charge on the induction plane crosses a threshold, create a TA
74 0 : else if(!m_collection_plane && m_current_window.adc_integral > m_adc_threshold_induction){
75 0 : TLOG_DEBUG(TLVL_DEBUG_LOW) << "[TAM:BSMW] ADC integral in window is greater than specified threshold.";
76 0 : output_ta.push_back(construct_ta());
77 0 : TLOG_DEBUG(TLVL_DEBUG_HIGH) << "[TAM:BSMW] Resetting window with input_tp.";
78 0 : m_current_window.reset(input_tp);
79 0 : }
80 :
81 : // If the addition of the current TP to the window would make it longer
82 : // than the specified window length, don't add it
83 : // Check the TPs are on the collection plane - if they are we can use XGBoost
84 : // Instead check whether it has been long enough since the last XGBoost prediction
85 : // then run the model to determine whether to create a TA
86 0 : else if (m_collection_plane &&
87 0 : (m_current_window.time_start - m_last_pred_time) > m_bin_length && // check enough time has passed since last window
88 0 : m_current_window.adc_integral > m_adc_threshold_collection && // set a low minimum threshold for the ADC integral sum
89 0 : compute_treelite_classification() // XGBoost classifier
90 : )
91 : {
92 0 : TLOG_DEBUG(TLVL_DEBUG_LOW) << "[TAM:BSMW] XGBoost neutrino prob. is greater than specified threshold.";
93 0 : output_ta.push_back(construct_ta());
94 0 : TLOG_DEBUG(TLVL_DEBUG_HIGH) << "[TAM:BSMW] Resetting window with input_tp.";
95 0 : m_current_window.reset(input_tp);
96 : }
97 : // If it is not, move the window along.
98 : else{
99 0 : TLOG_DEBUG(TLVL_DEBUG_ALL) << "[TAM:BSMW] Window is at required length but adc/bdt threshold not met, shifting window along.";
100 0 : m_current_window.move(input_tp, m_window_length);
101 : }
102 :
103 0 : TLOG_DEBUG(TLVL_DEBUG_ALL) << "[TAM:BSMW] " << m_current_window;
104 :
105 0 : m_primitive_count++;
106 :
107 0 : return;
108 :
109 : }
110 :
111 : void
112 0 : TAMakerProtoDUNEBSMWindowAlgorithm::configure(const nlohmann::json &config)
113 : {
114 0 : if (config.is_object()){
115 0 : if (config.contains("channel_map_name")) m_channel_map_name = config["channel_map_name"];
116 0 : if (config.contains("adc_threshold_induction")) m_adc_threshold_induction = config["adc_threshold_induction"];
117 0 : if (config.contains("bdt_threshold")) m_bdt_threshold = config["bdt_threshold"];
118 : }
119 : else{
120 0 : TLOG_DEBUG(TLVL_IMPORTANT) << "[TAM:BSMW] Use DEFAULT values of channel_map_name, adc_threshold_induction and bdt_threshold.";
121 : }
122 :
123 0 : TLOG_DEBUG(TLVL_DEBUG_ALL) << "[TAM:BSMW] Channel map name is " << m_channel_map_name <<
124 0 : ". ADC threshold for the induction planes set to " << m_adc_threshold_induction <<
125 0 : ". BDT threshold for collection plane set to " << m_bdt_threshold;
126 :
127 0 : channelMap = dunedaq::detchannelmaps::make_tpc_map(m_channel_map_name);
128 :
129 : // If we are in PD-VD, set boolean to true to enable effective channel mapping
130 0 : if (m_channel_map_name == "PD2VDTPCChannelMap" || m_channel_map_name == "PD2VDBottomTPCChannelMap" ||
131 0 : m_channel_map_name == "PD2VDTopTPCChannelMap") {
132 0 : m_pdvd_map = true;
133 : } else { // else we are in PD-HD and we use true channel mapping
134 0 : m_pdvd_map = false;
135 : }
136 :
137 0 : m_bin_length = static_cast<timestamp_t>(m_window_length / m_num_timebins);
138 :
139 0 : m_compiled_model_interface = std::make_unique<CompiledModelInterface>(nbatch, m_pdvd_map);
140 :
141 0 : const size_t num_feature = m_compiled_model_interface->GetNumFeatures();
142 :
143 0 : flat_batched_inputs.resize(num_feature);
144 :
145 0 : flat_batched_Entries.clear();
146 0 : for (size_t i = 0; i < num_feature; ++i) {
147 0 : union Entry zero;
148 0 : zero.fvalue = 0.0;
149 0 : flat_batched_Entries.emplace_back(zero);
150 : }
151 0 : }
152 :
153 : TriggerActivity
154 0 : TAMakerProtoDUNEBSMWindowAlgorithm::construct_ta() const
155 : {
156 0 : TLOG_DEBUG(TLVL_DEBUG_LOW) << "[TAM:BSMW] I am constructing a trigger activity!";
157 :
158 0 : TriggerPrimitive latest_tp_in_window = m_current_window.tp_list.back();
159 : // The time_peak, time_activity, channel_* and adc_peak fields of this TA are irrelevent
160 : // for the purpose of this trigger alg.
161 0 : TriggerActivity ta;
162 0 : ta.time_start = m_current_window.time_start;
163 0 : ta.time_end = latest_tp_in_window.time_start + latest_tp_in_window.samples_over_threshold * 32;
164 0 : ta.time_peak = latest_tp_in_window.samples_to_peak * 32 + latest_tp_in_window.time_start;
165 0 : ta.time_activity = ta.time_peak;
166 0 : ta.channel_start = latest_tp_in_window.channel;
167 0 : ta.channel_end = latest_tp_in_window.channel;
168 0 : ta.channel_peak = latest_tp_in_window.channel;
169 0 : ta.adc_integral = m_current_window.adc_integral;
170 0 : ta.adc_peak = latest_tp_in_window.adc_peak;
171 0 : ta.detid = latest_tp_in_window.detid;
172 0 : ta.type = TriggerActivity::Type::kTPC;
173 0 : ta.algorithm = TriggerActivity::Algorithm::kProtoDUNEBSMWindow;
174 0 : ta.inputs = m_current_window.tp_list;
175 0 : return ta;
176 0 : }
177 :
178 0 : bool TAMakerProtoDUNEBSMWindowAlgorithm::compute_treelite_classification() {
179 :
180 0 : m_last_pred_time = m_current_window.time_start;
181 :
182 0 : m_current_window.bin_window(
183 0 : flat_batched_inputs,
184 0 : m_num_timebins, m_bin_length,
185 0 : m_num_chanbins, m_n_channels_on_plane, m_first_channel,
186 0 : m_pdvd_eff_channel_mapper, m_pdvd_map
187 : );
188 :
189 0 : m_current_window.fill_entry_window(flat_batched_Entries, flat_batched_inputs);
190 :
191 0 : std::vector<float> result(nbatch, 0.0f);
192 :
193 0 : m_compiled_model_interface->Predict(flat_batched_Entries.data(), result.data());
194 :
195 0 : return m_compiled_model_interface->Classify(result.data(), m_bdt_threshold);
196 :
197 0 : }
198 :
199 : // Register algo in TA Factory
200 11 : REGISTER_TRIGGER_ACTIVITY_MAKER(TRACE_NAME, TAMakerProtoDUNEBSMWindowAlgorithm)
|
