Merge pull request #141 from wenqiang-gu/json_elecresp

add interface for tabulated electronics response

cfg/pgrapher/common/tools.jsonnet

@@ -80,6 +80,9 @@ function(params)
             shaping: params.elec.shaping,
             gain: params.elec.gain,
             postgain: params.elec.postgain,
+            filename: if std.objectHas(params.elec, "filename")
+                      then params.elec.filename
+                      else ""
         },
     },
 

cfg/pgrapher/experiment/pdsp/params.jsonnet

@@ -131,6 +131,11 @@ base {
                        // theoretical elec resp (14mV/fC): 36.6475 ADC*tick/1ke
       shaping: 2.2 * wc.us,
     },
+    // elec: super.elec {
+    //     type: "JsonElecResponse",
+    //     filename: "bnl-coldelec-response-gain14-shaping5.json.bz2",
+    //     postgain: 1.0,
+    // },
 
     sim: super.sim {
 

cfg/pgrapher/experiment/pdsp/wct-sim-check.jsonnet

@@ -49,8 +49,10 @@ local tracklist = [
 
   {
     time: 0 * wc.us, 
-    charge: -5000, // 5000 e/mm
-    ray: params.det.bounds,
+    // charge: -5000, // 5000 e/mm
+    // ray: params.det.bounds,
+    charge: -10000, // 5000 e/mm
+    ray: parallel,
   },
 
 ];
@@ -123,12 +125,12 @@ local parallel_pipes = [
   g.pipeline([
                sn_pipes[n],
                magnify_pipes[n],
-               nf_pipes[n],
-               magnify_pipes2[n],
-               sp_pipes[n],
-               magnify_pipes3[n],
-               magnify_pipes4[n],
-               magnify_pipes5[n],
+               // nf_pipes[n],
+               // magnify_pipes2[n],
+               // sp_pipes[n],
+               // magnify_pipes3[n],
+               // magnify_pipes4[n],
+               // magnify_pipes5[n],
              ],
              'parallel_pipe_%d' % n)
   for n in std.range(0, std.length(tools.anodes) - 1)

gen/inc/WireCellGen/JsonElecResponse.h

@@ -0,0 +1,44 @@
+/** An ElecResponse waveform is that due to sending a unit charge
+ * pulse into an amplifier with gain and shaping and a particular
+ * response function electronics. This JsonElecResponse defines a
+ * response function from a json input. */
+
+#ifndef WIRECELL_GEN_JSONELECRESPONSE
+#define WIRECELL_GEN_JSONELECRESPONSE
+
+#include "WireCellIface/IWaveform.h"
+#include "WireCellIface/IConfigurable.h"
+#include "WireCellUtil/Units.h"
+#include "WireCellUtil/Response.h"
+
+namespace WireCell {
+    namespace Gen {
+
+        class JsonElecResponse : public IWaveform, public IConfigurable {
+           public:
+            JsonElecResponse(int nticks = 10000, double t0 = 0,
+                             double postgain = 1.0, double tick = 0.5 * units::us,
+                             std::string filename="");
+            // IConfigurable interface
+            virtual void configure(const WireCell::Configuration& cfg);
+            virtual WireCell::Configuration default_configuration() const;
+
+            // The starting point of the sampling
+            virtual double waveform_start() const;
+            // The sampling period aka bin width
+            virtual double waveform_period() const;
+            // The collection of samples
+            virtual const sequence_type& waveform_samples() const;
+            // The collection of sample rebinned
+            virtual sequence_type waveform_samples(const WireCell::Binning& tbins) const;
+
+           private:
+            Configuration m_cfg;
+            mutable sequence_type m_wave;
+            Waveform::realseq_t m_times;
+            Waveform::realseq_t m_amps;
+        };
+    }  // namespace Gen
+}  // namespace WireCell
+
+#endif

gen/src/JsonElecResponse.cxx

@@ -0,0 +1,94 @@
+#include "WireCellGen/JsonElecResponse.h"
+
+#include "WireCellUtil/Persist.h"
+#include "WireCellUtil/NamedFactory.h"
+#include "WireCellUtil/Response.h"
+#include "WireCellUtil/Units.h"
+
+#include <iostream>
+#include <algorithm>
+
+WIRECELL_FACTORY(JsonElecResponse, WireCell::Gen::JsonElecResponse, WireCell::IWaveform, WireCell::IConfigurable)
+
+using namespace std;
+using namespace WireCell;
+
+// Generate waveform given a binning and a tablulated function
+WireCell::Waveform::realseq_t generate(const WireCell::Binning& tbins,
+    const Waveform::realseq_t& times, const Waveform::realseq_t& amps)
+{
+    const int nsamples = tbins.nbins();
+    WireCell::Waveform::realseq_t ret(nsamples, 0.0);
+    for (int ind = 0; ind < nsamples; ++ind) {
+        const double time = tbins.center(ind);
+        if (time < times.front() or time > times.back())
+            ret.at(ind) = 0;
+        else { // find the element that is greater or equal
+            auto it = std::lower_bound(times.begin(), times.end(), time);
+            auto k = it - times.begin();
+            if (*it == time) ret.at(ind) = amps.at(k);
+            else { // interpolation
+                ret.at(ind) = amps.at(k-1) + (amps.at(k)-amps.at(k-1))/(times.at(k)-times.at(k-1))*(time - times.at(k-1));
+            }
+        }
+    }
+    return ret;
+}
+
+Gen::JsonElecResponse::JsonElecResponse(int nticks,  double t0, double postgain,
+                                        double tick, std::string filename)
+{
+    m_cfg["filename"] = filename;
+    m_cfg["postgain"] = postgain;
+    m_cfg["start"] = t0;
+    m_cfg["tick"] = tick;
+    m_cfg["nticks"] = nticks;
+}
+WireCell::Configuration Gen::JsonElecResponse::default_configuration() const { return m_cfg; }
+void Gen::JsonElecResponse::configure(const WireCell::Configuration& cfg)
+{
+    m_cfg = cfg;
+    auto filename = m_cfg["filename"].asString();
+    if (filename.empty()) {
+        THROW(ValueError() << errmsg{"must supply a JsonElecResponse filename"});
+    }
+
+    auto top = Persist::load(filename);
+    auto jtimes = top["times"];
+    auto jamps = top["amplitudes"];
+    if (jtimes.size() != jamps.size()) {
+        THROW(ValueError() << errmsg{"inconsistent dimensions in file " + filename});
+    }
+    const int nsamp = jtimes.size();
+    m_times.resize(nsamp);
+    m_amps.resize(nsamp);
+
+    if (jtimes.isNull()) {
+        THROW(ValueError() << errmsg{"no values given in file " + filename});
+    }
+    const double tick = waveform_period();
+    for (int i=0; i<nsamp; i++) {
+        m_times.at(i) = jtimes[i].asFloat();
+        m_amps.at(i) = jamps[i].asFloat();
+    }
+
+    const int nbins = m_cfg["nticks"].asInt();
+    const double t0 = waveform_start();
+    Binning bins(nbins, t0, t0 + nbins * tick);
+    m_wave = generate(bins, m_times, m_amps);
+    Waveform::scale(m_wave, m_cfg["postgain"].asDouble());
+}
+
+double Gen::JsonElecResponse::waveform_start() const { return m_cfg["start"].asDouble(); }
+
+double Gen::JsonElecResponse::waveform_period() const { return m_cfg["tick"].asDouble(); }
+
+const IWaveform::sequence_type& Gen::JsonElecResponse::waveform_samples() const { return m_wave; }
+
+IWaveform::sequence_type Gen::JsonElecResponse::waveform_samples(const WireCell::Binning& tbins) const
+{
+    sequence_type rebinned_wave = generate(tbins, m_times, m_amps);
+    Waveform::scale(rebinned_wave, m_cfg["postgain"].asDouble());
+    return rebinned_wave;
+}
+