Merge pull request #180 from wenqiang-gu/rotate_deposet

Rotate deposet

cfg/pgrapher/experiment/pdsp/wct-sim-deposet-drifter.jsonnet

@@ -0,0 +1,157 @@
+local g = import 'pgraph.jsonnet';
+local f = import 'pgrapher/common/funcs.jsonnet';
+local wc = import 'wirecell.jsonnet';
+
+local io = import 'pgrapher/common/fileio.jsonnet';
+local tools_maker = import 'pgrapher/common/tools.jsonnet';
+local params = import 'pgrapher/experiment/pdsp/simparams.jsonnet';
+
+local tools = tools_maker(params);
+
+local sim_maker = import 'pgrapher/experiment/pdsp/sim.jsonnet';
+local sim = sim_maker(params, tools);
+
+local stubby = {
+  tail: wc.point(3000.0, 1000.0, 0.0, wc.mm),
+  head: wc.point(3000.0, 1000.0, 2300.0, wc.mm),
+};
+
+local tracklist = [
+
+  {
+    time: 0 * wc.us, 
+    charge: -5000, // 5000 e/mm
+    ray: stubby, // params.det.bounds,
+  },
+
+];
+local output = 'wct-sim-ideal-sig.npz';
+
+
+//local depos = g.join_sources(g.pnode({type:"DepoMerger", name:"BlipTrackJoiner"}, nin=2, nout=1),
+//                             [sim.ar39(), sim.tracks(tracklist)]);
+local depos = sim.tracks(tracklist, step=1.0 * wc.mm);
+local deposet = g.pnode({
+        type: "TrackDepoSet",
+        name: "trackdeposet",
+        data: {
+            step_size: 1.0*wc.mm,
+            tracks: tracklist,
+        }
+    }, nin=0, nout=1);
+
+
+//local deposio = io.numpy.depos(output);
+local drifter = sim.drifter;
+local setdrifter = g.pnode({
+            type: 'DepoSetDrifter',
+            data: {
+                drifter: "Drifter"
+            }
+        }, nin=1, nout=1, uses=[drifter]);
+local bagger = sim.make_bagger();
+// signal plus noise pipelines
+//local sn_pipes = sim.signal_pipelines;
+local sn_pipes = sim.splusn_pipelines;
+
+local multimagnify = import 'pgrapher/experiment/pdsp/multimagnify.jsonnet';
+local magoutput = 'case3.root';
+
+local multi_magnify = multimagnify('orig', tools, magoutput);
+local magnify_pipes = multi_magnify.magnify_pipelines;
+local multi_magnify2 = multimagnify('raw', tools, magoutput);
+local magnify_pipes2 = multi_magnify2.magnify_pipelines;
+local multi_magnify3 = multimagnify('gauss', tools, magoutput);
+local magnify_pipes3 = multi_magnify3.magnify_pipelines;
+local multi_magnify4 = multimagnify('wiener', tools, magoutput);
+local magnify_pipes4 = multi_magnify4.magnify_pipelines;
+local multi_magnify5 = multimagnify('threshold', tools, magoutput);
+local magnify_pipes5 = multi_magnify5.magnifysummaries_pipelines;
+
+local perfect = import 'pgrapher/experiment/pdsp/chndb-base.jsonnet';
+local chndb = [{
+  type: 'OmniChannelNoiseDB',
+  name: 'ocndbperfect%d' % n,
+  data: perfect(params, tools.anodes[n], tools.field, n),
+  uses: [tools.anodes[n], tools.field],  // pnode extension
+} for n in std.range(0, std.length(tools.anodes) - 1)];
+
+//local chndb_maker = import 'pgrapher/experiment/pdsp/chndb.jsonnet';
+//local noise_epoch = "perfect";
+//local noise_epoch = "after";
+//local chndb_pipes = [chndb_maker(params, tools.anodes[n], tools.fields[n]).wct(noise_epoch)
+//                for n in std.range(0, std.length(tools.anodes)-1)];
+local nf_maker = import 'pgrapher/experiment/pdsp/nf.jsonnet';
+// local nf_pipes = [nf_maker(params, tools.anodes[n], chndb_pipes[n]) for n in std.range(0, std.length(tools.anodes)-1)];
+local nf_pipes = [nf_maker(params, tools.anodes[n], chndb[n], n, name='nf%d' % n) for n in std.range(0, std.length(tools.anodes) - 1)];
+
+local sp_maker = import 'pgrapher/experiment/pdsp/sp.jsonnet';
+local sp = sp_maker(params, tools);
+local sp_pipes = [sp.make_sigproc(a) for a in tools.anodes];
+
+//local parallel_pipes = [g.pipeline([sn_pipes[n], magnify_pipes[n],
+//                                nf_pipes[n], magnify_pipes2[n] ], "parallel_pipe_%d"%n)
+//                    for n in std.range(0, std.length(tools.anodes)-1)];
+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],
+             ],
+             'parallel_pipe_%d' % n)
+  for n in std.range(0, std.length(tools.anodes) - 1)
+];
+local outtags = ['raw%d' % n for n in std.range(0, std.length(tools.anodes) - 1)];
+local parallel_graph = f.fanpipe('DepoSetFanout', parallel_pipes, 'FrameFanin', 'sn_mag_nf', outtags);
+
+
+//local frameio = io.numpy.frames(output);
+local sink = sim.frame_sink;
+
+// local graph = g.pipeline([depos, drifter, bagger, parallel_graph, sink]);
+local plainbagger = g.pnode({
+        type:'DepoBagger',
+        name:'plainbagger',
+        data: {
+            gate: [0, 0],
+        },
+    }, nin=1, nout=1);
+
+local deposet_rotate = g.pnode({
+        type:'DepoSetRotate',
+        name:'deposet_rotate',
+        data: {
+            rotate: true,
+            transpose: [1,0,2],
+            scale: [-1,1,1],
+        },
+    }, nin=1, nout=1);
+
+// FIXME: need a "bagger" after "setdrifter" for trimming out-of-window depos
+// local graph = g.pipeline([depos, plainbagger, setdrifter, parallel_graph, sink]);
+local graph = g.pipeline([depos, plainbagger, deposet_rotate, setdrifter, parallel_graph, sink]);
+
+local app = {
+  type: 'Pgrapher',
+  data: {
+    edges: g.edges(graph),
+  },
+};
+
+local cmdline = {
+    type: "wire-cell",
+    data: {
+        plugins: ["WireCellGen", "WireCellPgraph", "WireCellSio", "WireCellSigProc", "WireCellRoot"],
+        apps: ["Pgrapher"]
+    }
+};
+
+
+// Finally, the configuration sequence which is emitted.
+
+[cmdline] + g.uses(graph) + [app]

gen/inc/WireCellGen/DepoSetRotate.h

@@ -0,0 +1,55 @@
+//  This modifies coordinate of depos, in particular, when depo positions
+//  from upstream simulation (eg, LarSoft) has a non-X drift direction,
+//  one can rotate the coordinate and still has a X-drift internally in
+//  wire-cell-toolkit, and later rotate back to the original coordinate
+//  convention after the "drifter".
+//
+//  The current implementation allows a transpose of three axes as well as
+//  a scaling to mimic a rotation, eg, a transpose of [1,0,2] switches x-
+//  and y-axis, a scale of [1,-1, 1] maps y to minus-y. A combination of
+//  above two operations is equivalent to a 90-deg rotation around z-axis.
+//
+//  A more general rotation can be described in matrix, while the
+//  computation is expensive and can be implemented if needed.
+//
+
+#ifndef WIRECELLGEN_DEPOSETROTATE
+#define WIRECELLGEN_DEPOSETROTATE
+
+#include "WireCellIface/IDepoSetFilter.h"
+#include "WireCellIface/INamed.h"
+#include "WireCellIface/IConfigurable.h"
+#include "WireCellAux/Logger.h"
+
+namespace WireCell::Gen {
+
+    class DepoSetRotate : public Aux::Logger,
+                           public IDepoSetFilter, public IConfigurable {     
+      public:
+
+        DepoSetRotate();
+        virtual ~DepoSetRotate();
+
+        // IDepoSetFilter
+        virtual bool operator()(const input_pointer& in, output_pointer& out);
+
+        /// WireCell::IConfigurable interface.
+        virtual void configure(const WireCell::Configuration& config);
+        virtual WireCell::Configuration default_configuration() const;
+
+      private:
+
+        // Define how to transpose and scale coordinates
+        bool m_rotate{false};
+        std::tuple<int, int, int> m_transpose{0,1,2};
+        std::tuple<double, double, double> m_scale{1.,1.,1.};
+
+        size_t m_count{0};
+
+    };
+
+}
+
+
+
+#endif

gen/src/DepoSetRotate.cxx

@@ -0,0 +1,73 @@
+#include "WireCellGen/DepoSetRotate.h"
+#include "WireCellUtil/NamedFactory.h"
+#include "WireCellAux/SimpleDepoSet.h"
+#include "WireCellAux/SimpleDepo.h"
+
+WIRECELL_FACTORY(DepoSetRotate, WireCell::Gen::DepoSetRotate,
+                 WireCell::INamed,
+                 WireCell::IDepoSetFilter, WireCell::IConfigurable)
+
+
+using namespace WireCell;
+using namespace WireCell::Gen;
+
+DepoSetRotate::DepoSetRotate()
+    : Aux::Logger("DepoSetRotate", "gen")
+{
+}
+DepoSetRotate::~DepoSetRotate()
+{
+}
+
+WireCell::Configuration DepoSetRotate::default_configuration() const
+{
+    Configuration cfg;
+
+    // set the transpose and scale with a three element array
+    cfg["rotate"] = false;
+    cfg["transpose"] = Json::arrayValue;
+    cfg["scale"] = Json::arrayValue;
+    return cfg;
+}
+
+void DepoSetRotate::configure(const WireCell::Configuration& cfg)
+{
+    m_rotate = get(cfg, "rotate", false);
+
+    if (! cfg["transpose"].isNull()) {
+        m_transpose = std::tuple<int, int, int>(
+            cfg["transpose"][0].asInt(), cfg["transpose"][1].asInt(), cfg["transpose"][2].asInt());
+    }
+    if (! cfg["scale"].isNull()) {
+        m_scale = std::tuple<double, double, double>(
+            cfg["scale"][0].asDouble(), cfg["scale"][1].asDouble(), cfg["scale"][2].asDouble());
+    }
+}
+
+bool DepoSetRotate::operator()(const input_pointer& in, output_pointer& out)
+{
+    out = nullptr;
+    if (!in) {                  // EOS
+        log->debug("EOS at call={}", m_count);
+        return true;
+    }
+
+    IDepo::vector all_depos;
+    for (auto depo: *(in->depos())) {
+      if (m_rotate) {
+        auto pos = depo->pos();
+        WireCell::Point newpos;
+        newpos.x(pos[std::get<0>(m_transpose)] * std::get<0>(m_scale)); // apply transpose and scale
+        newpos.y(pos[std::get<1>(m_transpose)] * std::get<1>(m_scale));
+        newpos.z(pos[std::get<2>(m_transpose)] * std::get<2>(m_scale));
+        auto newdepo = std::make_shared<Aux::SimpleDepo>(depo->time(), newpos, depo->charge(), depo, depo->extent_long(), depo->extent_tran());
+        all_depos.push_back(newdepo);
+      }
+    }
+
+    log->debug("call={} rotated ndepos={}", m_count, all_depos.size());
+    out = std::make_shared<Aux::SimpleDepoSet>(m_count, all_depos);
+    ++m_count;
+    return true;
+}
+