Merge pull request #396 from WireCell/feature/wgu_filt_resp

Add filter on field response

sigproc/inc/WireCellSigProc/FilterResponse.h

@@ -0,0 +1,40 @@
+/** This component provides per-wire filtering of
+ * field responses based on a configuration data file. */
+
+#ifndef WIRECELLSIGPROC_FILTERRESPONSE
+#define WIRECELLSIGPROC_FILTERRESPONSE
+
+#include "WireCellIface/IChannelResponse.h"
+#include "WireCellIface/IConfigurable.h"
+#include "WireCellUtil/Units.h"
+
+#include <string>
+#include <unordered_map>
+
+namespace WireCell {
+    namespace SigProc {
+        class FilterResponse : public IChannelResponse, public IConfigurable {
+           public:
+            FilterResponse(const char* filename = "", const int planeid = 0);
+
+            virtual ~FilterResponse();
+
+            // IChannelResponse
+            virtual const Waveform::realseq_t& channel_response(int channel_ident) const;
+            virtual Binning channel_response_binning() const;
+
+            // IConfigurable
+            virtual void configure(const WireCell::Configuration& config);
+            virtual WireCell::Configuration default_configuration() const;
+
+           private:
+            std::string m_filename;
+            int m_planeid;
+            std::unordered_map<int, Waveform::realseq_t> m_cr;
+            Binning m_bins;
+        };
+
+    }  // namespace SigProc
+
+}  // namespace WireCell
+#endif

sigproc/inc/WireCellSigProc/OmnibusSigProc.h

@@ -222,6 +222,8 @@ namespace WireCell {
 
             // average overall responses
             std::vector<Waveform::realseq_t> overall_resp[3];
+            // filters for overall responses
+            std::vector<std::string> m_filter_resps_tn{};
 
             // tag name for traces
             std::string m_wiener_tag{"wiener"};

sigproc/src/FilterResponse.cxx

@@ -0,0 +1,71 @@
+#include "WireCellSigProc/FilterResponse.h"
+
+#include "WireCellUtil/Persist.h"
+#include "WireCellUtil/Exceptions.h"
+#include "WireCellUtil/String.h"
+#include "WireCellUtil/Response.h"
+
+#include "WireCellUtil/NamedFactory.h"
+
+WIRECELL_FACTORY(FilterResponse, WireCell::SigProc::FilterResponse, WireCell::IChannelResponse,
+                 WireCell::IConfigurable)
+
+using namespace WireCell;
+
+SigProc::FilterResponse::FilterResponse(const char* filename, const int planeid)
+  : m_filename(filename), m_planeid(planeid)
+{
+}
+
+SigProc::FilterResponse::~FilterResponse() {}
+
+WireCell::Configuration SigProc::FilterResponse::default_configuration() const
+{
+    Configuration cfg;
+    cfg["filename"] = m_filename;
+    cfg["planeid"] = m_planeid;
+    return cfg;
+}
+
+void SigProc::FilterResponse::configure(const WireCell::Configuration& cfg)
+{
+
+    m_filename = get(cfg, "filename", m_filename);
+    if (m_filename.empty()) {
+        THROW(ValueError() << errmsg{"must supply a FilterResponse filename"});
+    }
+    m_planeid = get(cfg, "planeid", m_planeid);
+
+    auto top = Persist::load(m_filename);
+    const int nwires = top["nwires"].asInt();
+    const int nticks = top["nticks"].asInt();
+    if (!m_bins.nbins()) {  // first time
+        m_bins = Binning(nticks, 0, nticks); // tick range for filter
+    }
+    auto jflts = top["filters"];
+    if (jflts.isNull()) {
+        THROW(ValueError() << errmsg{"no filters given in file " + m_filename});
+    }
+
+    for (auto jflt : jflts) {
+        const int plane = jflt["plane"].asInt();
+        if (plane != m_planeid) continue;
+
+        const int wire = jflt["wire"].asInt();
+        Waveform::realseq_t resp(nticks, 0.0);
+        std::transform(jflt["values"].begin(), jflt["values"].end(), resp.begin(),
+                       [](const auto& v) { return v.asFloat(); });
+        m_cr[wire] = resp;
+    }
+}
+
+const Waveform::realseq_t& SigProc::FilterResponse::channel_response(int channel_ident) const
+{
+    const auto& it = m_cr.find(channel_ident);
+    if (it == m_cr.end()) {
+        THROW(KeyError() << errmsg{String::format("no response for channel %d", channel_ident)});
+    }
+    return it->second;
+}
+
+Binning SigProc::FilterResponse::channel_response_binning() const { return m_bins; }

sigproc/src/OmnibusSigProc.cxx

@@ -154,6 +154,13 @@ void OmnibusSigProc::configure(const WireCell::Configuration& config)
 
     m_charge_ch_offset = get(config, "charge_ch_offset", m_charge_ch_offset);
 
+    if (config.isMember("filter_responses_tn")) {
+        m_filter_resps_tn.clear();
+        for (auto tn: config["filter_responses_tn"]) {
+            m_filter_resps_tn.push_back(tn.asString());
+        }
+    }
+
     m_wiener_tag = get(config, "wiener_tag", m_wiener_tag);
     // m_wiener_threshold_tag = get(config, "wiener_threshold_tag", m_wiener_threshold_tag);
     if (! config["wiener_threshold_tag"].isNull()) {
@@ -1624,6 +1631,18 @@ bool OmnibusSigProc::operator()(const input_pointer& in, output_pointer& out)
         // load data into EIGEN matrices ...
         load_data(in, iplane);  // load into a large matrix
         // initial decon ...
+
+        // additional filters for overall resposne
+        if (!m_filter_resps_tn.empty()) {
+            for (size_t i = 0; i != overall_resp[iplane].size(); i++) {
+                auto fltresp = Factory::find_tn<IChannelResponse>(m_filter_resps_tn[iplane]);
+                const Waveform::realseq_t& flt = fltresp->channel_response(i); // filter at wire: i
+                for (int j = 0; j != std::min<int>(m_fft_nticks, flt.size()); j++) {
+                    overall_resp[iplane].at(i).at(j) *= flt.at(j); 
+                }
+            }
+        }
+
         decon_2D_init(iplane);  // decon in large matrix
         check_data(iplane, "after 2D init");
 
@@ -1782,6 +1801,13 @@ bool OmnibusSigProc::operator()(const input_pointer& in, output_pointer& out)
                 wiener_traces.insert(wiener_traces.end(), perframe.begin(), perframe.end());
             }
 
+            if (!m_filter_resps_tn.empty()) {
+                // reload data and field response
+                init_overall_response(in);
+                load_data(in, iplane); 
+                decon_2D_init(iplane);  // decon in large matrix
+            }
+
             decon_2D_charge(iplane);
             std::vector<double> dummy_thresholds;
             if (m_use_roi_debug_mode and !m_decon_charge_tag.empty()) {