Bsm trigger dev/refactor PR for update to ProtoDUNEBSMWindow TAMaker

CMakeLists.txt

@@ -62,6 +62,7 @@ add_library(triggeralgs SHARED
   src/ProtoDUNEBSMWindow/ProtoDUNEBSMWindow.cpp
   src/ProtoDUNEBSMWindow/treelitemodel.cpp
   src/ProtoDUNEBSMWindow/models/treelitePDHDmodel.cpp
+  src/ProtoDUNEBSMWindow/models/treelitePDVDmodel.cpp
   src/ProtoDUNEBSMWindow/CompiledModelInterface.cpp
   src/TAWindow.cpp
   src/TPWindow.cpp

include/triggeralgs/ProtoDUNEBSMWindow/CompiledModelInterface.hpp

@@ -17,20 +17,18 @@ class TreeliteModelBase;
 class CompiledModelInterface {
   public:
 
-    CompiledModelInterface(int nbatch);
+    CompiledModelInterface(int nbatch, bool is_pdvd);
 
     ~CompiledModelInterface();
 
     // Get number of features in model 
     int GetNumFeatures();
 
-    void ModelWarmUp(Entry *input);
-
     // Run prediction with GBDT
     void Predict(Entry *input, float *result);
 
     // Is it a neutrino or cosmic according to GBDT?
-    bool Classify(const float *result, float &bdt_threshold);
+    bool Classify(const float *result, float bdt_threshold);
 
   protected:
 

include/triggeralgs/ProtoDUNEBSMWindow/ProtoDUNEBSMWindow.hpp

@@ -27,23 +27,18 @@ class ProtoDUNEBSMWindow {
 
     // Bins the TPs as a function of channel and time
     // If running in PD-VD then effective channel ID is used
-    void bin_window(std::vector<float> &input, timestamp_t time_bin_width, channel_t chan_bin_width, 
-                    int num_time_bins, int num_chan_bins, channel_t first_channel,
-                    std::unique_ptr<PDVDEffectiveChannelMap> const &effective_channel_mapper, bool use_pdvd_map);
+    void bin_window(std::vector<float> &input, 
+                    int num_time_bins, timestamp_t time_bin_width,
+                    int num_chan_bins, channel_t n_channels_on_plane, channel_t first_channel,
+                    std::unique_ptr<PDVDEffectiveChannelMap> const &effective_channel_mapper,
+                    bool use_pdvd_map);
 
     void fill_entry_window(std::vector<Entry> &entry_input, std::vector<float> &input);
 
-    // Calculate average properties of TPs in a time window
-    float mean_sadc();
-    float mean_adc_peak();
-    float mean_tot();
-
     friend std::ostream& operator<<(std::ostream& os, const ProtoDUNEBSMWindow& window);
 
     timestamp_t time_start;
     uint32_t adc_integral;
-    uint64_t adc_peak_sum;
-    uint64_t tot_sum;
     std::vector<TriggerPrimitive> tp_list;
   };
 }

include/triggeralgs/ProtoDUNEBSMWindow/TAMakerProtoDUNEBSMWindowAlgorithm.hpp

@@ -28,7 +28,7 @@ class TAMakerProtoDUNEBSMWindowAlgorithm : public TriggerActivityMaker
   void process(const TriggerPrimitive& input_tp, std::vector<TriggerActivity>& output_ta); 
   void configure(const nlohmann::json &config);
 
-  ~TAMakerProtoDUNEBSMWindowAlgorithm() override;
+  ~TAMakerProtoDUNEBSMWindowAlgorithm() override = default;
 
 private:
   // Function to handle XGBoost classification
@@ -40,7 +40,7 @@ class TAMakerProtoDUNEBSMWindowAlgorithm : public TriggerActivityMaker
   // The current time window of TPs
   ProtoDUNEBSMWindow m_current_window;
 
-  timestamp_t m_last_pred_time;
+  timestamp_t m_last_pred_time = 0;
   uint64_t m_primitive_count = 0;
 
   // Possible to do batch predictions with XGBoost
@@ -52,19 +52,24 @@ class TAMakerProtoDUNEBSMWindowAlgorithm : public TriggerActivityMaker
   std::vector<Entry> flat_batched_Entries;
 
   // Configurable parameters.
-  uint32_t m_adc_threshold = 200000;
-  float m_ratio_threshold = 0.65;
-  float m_bdt_threshold = 0.99;
-  timestamp_t m_window_length = 20000;
+  uint32_t m_adc_threshold_induction = 12000000;
+  float m_bdt_threshold = 0.999;
   std::string m_channel_map_name = "PD2VDTPCChannelMap";
   // End of configurable parameters
+
+  // Constant parameters defined by XGBoost model training
+  // Currently both PD-HD and PD-VD XGBoost models have been trained
+  // on data filtered to have an ADC sum > 200k ADC
+  const uint32_t m_adc_threshold_collection = 200000;
+  const timestamp_t m_window_length = 20000;
+  // End constant parameters
 
   // Define time binning
-  timestamp_t m_bin_length = 4000;
-  int m_num_timebins = 5;
-  // Define channel binning
-  channel_t m_chan_bin_length = 100;
-  int m_num_chanbins = 5;
+  timestamp_t m_bin_length = 2000;
+  // Number of bins fixed by model training
+  const int m_num_timebins = 10;
+  // Number of bins fixed by model training
+  const int m_num_chanbins = 10;
 
   // Geometry information for binning
   std::shared_ptr<dunedaq::detchannelmaps::TPCChannelMap> channelMap; 
@@ -78,9 +83,12 @@ class TAMakerProtoDUNEBSMWindowAlgorithm : public TriggerActivityMaker
   std::unique_ptr<PDVDEffectiveChannelMap> m_pdvd_eff_channel_mapper = nullptr;
   // If in NP02 and using a PD-VD channel map, set this to true
   bool m_pdvd_map = true;
+  // Only use the XGBoost model if we are looking at collection plane TPs
+  bool m_collection_plane = false;
   // first and last channel on the plane
   channel_t m_first_channel;
-  channel_t m_last_channel;
+  // Number of channels on the plane (will be no. effective channels for PD-VD)
+  channel_t m_n_channels_on_plane;
 
   // Compiled treelite model interface
   std::unique_ptr<CompiledModelInterface> m_compiled_model_interface;

include/triggeralgs/ProtoDUNEBSMWindow/treelitemodel.hpp

@@ -57,8 +57,14 @@ namespace triggeralgs {
     virtual int32_t get_num_feature() const override;
     void predict(union Entry* data, int pred_margin, float* result) const override;
   };
-
-  // Add model later for PD-VD signal findings
+
+  /// Model for signal finding in PD-VD — trained on PD-VD data
+  class TreelitePDVDModel : public TreeliteModelBase {
+  public:
+    using TreeliteModelBase::TreeliteModelBase;
+    virtual int32_t get_num_feature() const override;
+    void predict(union Entry* data, int pred_margin, float* result) const override;
+  };
 
 } // namespace triggeralgs
 #endif // TRIGGERALGS_TREELITEMODEL_HPP_

src/ProtoDUNEBSMWindow/CompiledModelInterface.cpp

@@ -4,8 +4,12 @@
 
 namespace triggeralgs {
 
-CompiledModelInterface::CompiledModelInterface(int nbatch) : num_batch(nbatch) {
-  model_ptr = std::make_unique<TreelitePDHDModel>();
+CompiledModelInterface::CompiledModelInterface(int nbatch, bool is_pdvd) : num_batch(nbatch) {
+  if (is_pdvd) {  
+    model_ptr = std::make_unique<TreelitePDVDModel>();
+  } else {
+    model_ptr = std::make_unique<TreelitePDHDModel>();
+  }
 }
 
 CompiledModelInterface::~CompiledModelInterface() {}
@@ -14,23 +18,13 @@ int CompiledModelInterface::GetNumFeatures() {
   return model_ptr->get_num_feature();
 }
 
-void CompiledModelInterface::ModelWarmUp(Entry *input) {
-  // Warm the BDT up here
-  float result[num_batch];
-  for (int rid = 0; rid < num_batch; ++rid) {
-    for (int i = 0; i < 100; i++) {
-      model_ptr->predict(input, 0, result);
-    }
-  }
-}
-
 void CompiledModelInterface::Predict(Entry *input, float *result) {
   for (int rid = 0; rid < num_batch; ++rid) {
     model_ptr->predict(input, 0, result);
   }
 }
 
-bool CompiledModelInterface::Classify(const float *result, float &bdt_threshold) {
+bool CompiledModelInterface::Classify(const float *result, float bdt_threshold) {
   for (uint64_t rid = 0; rid < num_batch; rid++) {
     if (result[rid] > bdt_threshold) {
       return true;

src/ProtoDUNEBSMWindow/DetectorPlaneMap.cpp

@@ -7,8 +7,10 @@ PlaneInfo DetectorPlaneMap::get_plane_info(std::string channel_map_name, int det
 
   const std::map<std::pair<int,int>, PlaneInfo> *plane_map;
   if (channel_map_name == "PD2HDTPCChannelMap") plane_map = &pdhd_plane_map;
-  else if (channel_map_name == "PD2VDTPCChannelMap") plane_map = &pdvd_plane_map;
-  else plane_map = &pdvd_plane_map;
+  else if (channel_map_name == "PD2VDTPCChannelMap" ||
+           channel_map_name == "PD2VDBottomTPCChannelMap" ||
+           channel_map_name == "PD2VDTopTPCChannelMap") plane_map = &pdvd_plane_map;
+  else throw std::invalid_argument("Unknown channel map: " + channel_map_name);
 
   auto it = plane_map->find({detelement, plane});
   if (it == plane_map->end()) {

src/ProtoDUNEBSMWindow/ProtoDUNEBSMWindow.cpp

@@ -11,12 +11,9 @@ bool ProtoDUNEBSMWindow::is_empty() const{
 };
 
 void ProtoDUNEBSMWindow::add(TriggerPrimitive const &input_tp){
-  // Add the input TP's contribution to the total ADC and add it to
-  // the TP list. Also keep running sum of all the samples over threshold
-  // and the peak ADC. These are used for samples/peak ratio cut
+  // Add the input TP's contribution to the total 
+  // ADC and add it to the TP list.
   adc_integral += input_tp.adc_integral;
-  adc_peak_sum += input_tp.adc_peak;
-  tot_sum += input_tp.samples_over_threshold;
   tp_list.push_back(input_tp);
 };
 
@@ -34,16 +31,14 @@ void ProtoDUNEBSMWindow::move(TriggerPrimitive const &input_tp, timestamp_t cons
     if(!(input_tp.time_start-tp.time_start < window_length)){
       n_tps_to_erase++;
       adc_integral -= tp.adc_integral;
-      adc_peak_sum -= tp.adc_peak;
-      tot_sum -= tp.samples_over_threshold;
     }
     else break;
   }
   // Erase the TPs from the window.
   tp_list.erase(tp_list.begin(), tp_list.begin()+n_tps_to_erase);
   // Make the window start time the start time of what is now the
   // first TP.
-  if(tp_list.size()!=0){
+  if(!tp_list.empty()){
     time_start = tp_list.front().time_start;
     add(input_tp);
   }
@@ -63,15 +58,15 @@ void ProtoDUNEBSMWindow::reset(TriggerPrimitive const &input_tp){
 };
 
 void ProtoDUNEBSMWindow::bin_window(
-    std::vector<float> &input, timestamp_t time_bin_width, 
-    channel_t chan_bin_width, int num_time_bins, 
-    int num_chan_bins, channel_t first_channel,
-    std::unique_ptr<PDVDEffectiveChannelMap> const &effective_channel_mapper, 
+    std::vector<float> &input, 
+    int num_time_bins, timestamp_t time_bin_width,
+    int num_chan_bins, channel_t n_channels_on_plane, channel_t first_channel,
+    std::unique_ptr<PDVDEffectiveChannelMap> const &effective_channel_mapper,
     bool use_pdvd_map) {
+
   std::fill(input.begin(), input.end(), 0.0f);
 
   const float inv_time_bin_width = 1.0f / time_bin_width;
-  const float inv_chan_bin_width = 1.0f / chan_bin_width;
 
   for (const TriggerPrimitive& tp : tp_list) {
     channel_t temp_tp_channel = tp.channel;
@@ -80,7 +75,9 @@ void ProtoDUNEBSMWindow::bin_window(
       temp_tp_channel = effective_channel_mapper->remapCollectionPlaneChannel(temp_tp_channel);
     }
     size_t time_bin = static_cast<size_t>((tp.time_start - time_start) * inv_time_bin_width);
-    size_t channel_bin = static_cast<size_t>((temp_tp_channel - first_channel) * inv_chan_bin_width);
+    // Channel bin calculation matches model training - remainder channels are appended to 10th channel bin 
+    size_t channel_bin = static_cast<size_t>(((temp_tp_channel - first_channel) * num_chan_bins) / n_channels_on_plane);
+
     if (time_bin < num_time_bins && channel_bin < num_chan_bins) {
       size_t index = channel_bin * num_time_bins + time_bin;
       input[index] += tp.adc_integral;
@@ -95,16 +92,6 @@ void ProtoDUNEBSMWindow::fill_entry_window(std::vector<Entry> &entry_input, std:
   }
 }
 
-float ProtoDUNEBSMWindow::mean_sadc() {
-  return static_cast<float>(adc_integral / tp_list.size());;
-}
-float ProtoDUNEBSMWindow::mean_adc_peak() {
-  return static_cast<float>(adc_peak_sum / tp_list.size());
-}
-float ProtoDUNEBSMWindow::mean_tot() {
-  return static_cast<float>(tot_sum / tp_list.size());
-}
-
 std::ostream& operator<<(std::ostream& os, const ProtoDUNEBSMWindow& window){
   if(window.is_empty()) os << "Window is empty!\n";
   else{