diff --git a/aux/docs/frame-files.org b/aux/docs/frame-files.org index 797c9ca5c..2f3985103 100644 --- a/aux/docs/frame-files.org +++ b/aux/docs/frame-files.org @@ -53,6 +53,12 @@ summary__.npy frame__.npy channels__.npy tickinfo__.npy +frame__.npy +channels__.npy +tickinfo__.npy +frame_*_.npy +channels_*_.npy +tickinfo_*_.npy chanmask__.npy chanmask__.npy ... @@ -60,19 +66,23 @@ chanmask__.npy We will call "framelet" the trio of arrays of type ~frame~, ~channels~ and ~tickinfo~ and an optional fourth ~summary~ array that carry the same -~~ and frame ~~. To be valid, the size of the 1D ~channels~ -array, and the summary array if present, must be the same as the -number of (channel) rows in the ~frame~ array. The first two entries -(~time~ and ~tick~) of all ~tickinfo~ are expected to be identical but their -third (~tbin0~) may differ for each framelet. See [[file:tickinfo.org]]. - -The rows of the ~frame~ and ~channel~ arrays are added to the ~IFrame~ -traces collection in the order they are encountered in the file. -Otherwise, order does not matter except that all arrays for a given -frame ~~ must be contiguous. The ~~ portion of the framelet -array file names is used to add to the ~IFrame~ trace map map. In the -current version of the file, there is no support for representing -tagged trace indices. +~~ and frame ~~. + +To be valid, the size of the 1D ~channels~ array, and the summary array if +present, must be the same as the number of (channel) rows in the ~frame~ array. +The first two entries (~time~ and ~tick~) of all ~tickinfo~ must be be identical but +their third entry (~tbin0~) may differ for each framelet. See [[file:tickinfo.org]]. + +The rows of the ~frame~ array and the corresponding elements of the ~channel~ array +are added to the ~IFrame~ traces collection in the order they are encountered in +the file. Otherwise, order does not matter except that all arrays for a given +frame ~~ must be contiguous. + +The ~~ portion of the framelet provides a *trace tag* (and *not* a frame tag). +When the ~~ is unspecified (ie, there are a pair of subsequent underscores, +eg ~frame__0~~) or when the ~~ takes a special value of ~*~ (eg ~frame_*_0~) then +the arrays are interpreted to represent an portion of the frame consisting of +untagged traces. In such a case, any associated ~summary_~ is ignored. The ~chanmask~ arrays are optional and each provides one entry in the ~IFrame~ channel mask map. The key is take as the ~~ part of the diff --git a/cfg/layers/README.org b/cfg/layers/README.org index 9d1def939..54e6d178e 100644 --- a/cfg/layers/README.org +++ b/cfg/layers/README.org @@ -3,13 +3,13 @@ * TODO +- [ ] move the emerging common pattern of ~mids/*/variants.jsonnet~ files into ~mids.jsonnet~. - [ ] add ~mains~ section - [ ] add ~options~ - [ ] generate dots by cross product of detector type, variants and mains - [ ] add ~uboone~ - [ ] remove ~helpers~ - * Introduction The WCT configuration structure layers convention defines a Jsonnet diff --git a/cfg/layers/gen-mids.sh b/cfg/layers/gen-mids.sh index 4c93b5e53..70cb5062b 100755 --- a/cfg/layers/gen-mids.sh +++ b/cfg/layers/gen-mids.sh @@ -17,10 +17,16 @@ cat < "$out" EOF # this loop should keep path relative -for one in mids/*/mids.jsonnet +for api in mids/*/api.jsonnet do - det=$(basename $(dirname $one)) - echo " $det : import \"$one\"," >> "$out" + det=$(basename $(dirname $api)) + var="$(dirname $api)/variants.jsonnet" + cat <> "$out" + $det : { + api: import "$api", + variants: import "$var", + }, +EOF done echo "}" >> "$out" diff --git a/cfg/layers/gengeom.sh b/cfg/layers/gengeom.sh new file mode 100644 index 000000000..98bccc3fb --- /dev/null +++ b/cfg/layers/gengeom.sh @@ -0,0 +1,12 @@ +#!/usr/bin/bash +me="$(realpath "$BASH_SOURCE")" +mydir="$(dirname "$me")" + +cd "$mydir" + +for det in base pdsp uboone +do + wirecell-util wires-volumes $det \ + | jsonnetfmt -n 4 --max-blank-lines 1 - \ + > $mydir/mids/$det/variants/geometry.jsonnet +done diff --git a/cfg/layers/high.jsonnet b/cfg/layers/high.jsonnet index 8c2d4fa30..0bd076548 100644 --- a/cfg/layers/high.jsonnet +++ b/cfg/layers/high.jsonnet @@ -4,12 +4,11 @@ // A lookup by detector and variant. See gen-mids.sh. local mids = import "mids.jsonnet"; +local base = import "mids/base/api.jsonnet"; local svcs = import "high/svcs.jsonnet"; local midapi = import "midapi.jsonnet"; - - { // Forward the various "standard" wirecell utillities as a // convenience to end user. @@ -30,9 +29,50 @@ local midapi = import "midapi.jsonnet"; // the mid. services :: svcs, - // The mid-level API factory function. - mid :: function(detector, variant="nominal", services=svcs(), options={}) - midapi + mids[detector](services, variant, options=options), + // All known mid-layer detectors + mids :: mids, + + // The parameter object factory. + // + // The parameter object is selected from the detector's mid-level + // variants.jsonnet structure based on the variant name. + // + // All three arguments are saved into the selected parameter object with a + // key formed by appending "_name", eg "detector_name". + // + // The detector mid API may use "params.structure_name" to supply an API + // that results in variant DFP graph structure independent from the chosen + // parameter variant represented by "params.variant_name". + // + // The following variant and structure names SHALL be implemented according + // to these guiding definitions: + // + // - nominal :: In name only. This is for defining an idealized + // configuration lacking variance. It shall avoid introducing magic numbers + // such as arbitrary time offsets, bad or per channel variance and shall + // consider a uniform detector and resulting DFP graph structure. + // + // - actual :: As close to the needs of realism as can be. This + // configuration shall endeavor to produce a configuration that model the + // real world detector to the extent possible. Strive to construct + // "actual" parameter variant and structure that derive from nominal. + // + // Other variant and structure names are allowed for special purposes. + params :: function(detector, variant="nominal", structure=null) + mids[detector].variants[variant] { + detector_name: detector, + variant_name: variant, + structure_name: if std.type(structure) == "null" then variant else structure, + }, + + // Return a mid-level API implementation inside a midapi API. + api :: function(detector, params, services=svcs(), options={}) + local sv = svcs(); // fixme: need a way to choose eg GPU svcs + //local base = midapi(sv, params, options=options); + local def = base(sv, params, options=options); + local det = mids[detector].api(sv, params, options=options); + local imp = std.mergePatch(def, det); + midapi(imp), } diff --git a/cfg/layers/low.jsonnet b/cfg/layers/low.jsonnet index cb12a8830..17faf337a 100644 --- a/cfg/layers/low.jsonnet +++ b/cfg/layers/low.jsonnet @@ -6,9 +6,15 @@ util : import "low/util.jsonnet", gen : import "low/gen.jsonnet", - drifter : import "low/drifter.jsonnet", + anodes : import "low/anodes.jsonnet", + + drifter : import "low/drifter.jsonnet", + reframer : import "low/reframer.jsonnet", dnnroi : import "low/dnnroi.jsonnet", img : import "low/img.jsonnet", + + resps : import "low/resps.jsonnet", + ssss : import "low/ssss.jsonnet", } diff --git a/cfg/layers/low/drifter.jsonnet b/cfg/layers/low/drifter.jsonnet index cde0f792f..0b2ce06d4 100644 --- a/cfg/layers/low/drifter.jsonnet +++ b/cfg/layers/low/drifter.jsonnet @@ -27,15 +27,17 @@ function(random, xregions, lar, offset=0, fluctuate=true, name="") pg.pnode({ type: 'DepoSetDrifter', name: name, - data: { drifter: "Drifter" } - }, nin=1, nout=1, uses=[pg.pnode({ - type: 'Drifter', - name: name, - data: lar { - rng: wc.tn(random), - xregions: xregions, // for help, see: low.util.driftsToXregions(drifts) - time_offset: offset, - fluctuate: fluctuate, - }, - }, nin=1, nout=1, uses=[random])]) + data: { drifter: "Drifter:"+name } + }, nin=1, nout=1, uses = [ + pg.pnode({ + type: 'Drifter', + name: name, + data: lar { + rng: wc.tn(random), + xregions: xregions, // for help, see: low.util.driftsToXregions(drifts) + time_offset: offset, + fluctuate: fluctuate, + }, + }, nin=1, nout=1, uses=[random]) + ]) diff --git a/cfg/layers/low/img.jsonnet b/cfg/layers/low/img.jsonnet index 8a24fb217..996c00fda 100644 --- a/cfg/layers/low/img.jsonnet +++ b/cfg/layers/low/img.jsonnet @@ -5,10 +5,7 @@ local util = import "util.jsonnet"; local pg = import "pgraph.jsonnet"; local wc = import "wirecell.jsonnet"; -function (anode) { - - local ident = util.idents(anode), - +function (anode, name) { // The recommended slicer. There is almost no reasonable default // to each detector variant may as well create the MaskSlices @@ -17,11 +14,11 @@ function (anode) { span=4, active_planes=[0,1,2], masked_planes=[], dummy_planes=[]) pg.pnode({ type: "MaskSlices", - name: ident+ext, + name: name+ext, data: { - wiener_tag: "wiener"+ident, - charge_tag: "gauss"+ident, - error_tag: "gauss_error"+ident, + wiener_tag: "wiener"+ext, + charge_tag: "gauss"+ext, + error_tag: "gauss_error"+ext, tick_span: span, anode: wc.tn(anode), min_tbin: min_tbin, @@ -36,7 +33,7 @@ function (anode) { // This slicing does not handle charge uncertainty nor bad channels. simple_slicing :: function(span=4, tag="") pg.pnode({ type: "SumSlices", - name: ident, + name: name, data: { tag: tag, tick_span: span, @@ -49,7 +46,7 @@ function (anode) { // unique name extension "ext" must be given. single_tiling :: function(face, ext="") pg.pnode({ type: "GridTiling", - name: "%s-%d%s"%[ident, face, ext], + name: "%s-%d%s"%[name, face, ext], data: { anode: wc.tn(anode), face: face, @@ -63,7 +60,7 @@ function (anode) { local slice_fanout = pg.pnode({ type: "SliceFanout", - name: ident+ext, + name: name+ext, data: { multiplicity: 2 }, }, nin=1, nout=2), @@ -71,7 +68,7 @@ function (anode) { local blobsync = pg.pnode({ type: "BlobSetSync", - name: ident+ext, + name: name+ext, data: { multiplicity: 2 } }, nin=2, nout=1), @@ -82,7 +79,7 @@ function (anode) { edges= [pg.edge(slice_fanout, tilings[n], n, 0) for n in [0,1]] + [pg.edge(tilings[n], blobsync, 0, n) for n in [0,1]], - name=ident+ext), + name=name+ext), }.ret, // A tiling subgraph matching the anode faces. The "uname" must @@ -229,7 +226,7 @@ function (anode) { // of slice spans. Makes blob-blob edges clustering :: function(spans=1.0) pg.pnode({ type: "BlobClustering", - name: ident, + name: name, data: { spans : spans } }, nin=1, nout=1), @@ -237,7 +234,7 @@ function (anode) { // blob-measure edges. grouping :: function() pg.pnode({ type: "BlobGrouping", - name: ident, + name: name, data: { } }, nin=1, nout=1), @@ -248,7 +245,7 @@ function (anode) { whiten=true) pg.pnode({ type: "ChargeSolving", - name: ident, + name: name, data: { "meas_value_threshold": meas_val_thresh, "meas_error_threshold": meas_err_thresh, @@ -260,14 +257,14 @@ function (anode) { // This solver is simplistic, prefer charge solving. blob_solving :: function(threshold=0.0) pg.pnode({ type: "BlobSolving", - name: ident, + name: name, data: { threshold: threshold } }, nin=1, nout=1), // A function that projects blobs back to frames. reframing :: function(tag="") pg.pnode({ type: "BlobReframer", - name: ident, + name: name, data: { frame_tag: tag, } diff --git a/cfg/layers/low/reframer.jsonnet b/cfg/layers/low/reframer.jsonnet new file mode 100644 index 000000000..91e13cde5 --- /dev/null +++ b/cfg/layers/low/reframer.jsonnet @@ -0,0 +1,18 @@ +local wc = import "wirecell.jsonnet"; +local pg = import "pgraph.jsonnet"; +local util = import "util.jsonnet"; + +function(params, anode, tags=[], name=null) + pg.pnode({ + type: 'Reframer', + name: if std.type(name) == "null" then util.idents(anode) else name, + data: { + anode: wc.tn(anode), + tags: tags, + fill: 0.0, + toffset: 0, + tbin: params.ductor.tbin, + nticks: params.ductor.binning.nticks - self.tbin, + }, + }, nin=1, nout=1, uses=[anode]) + diff --git a/cfg/layers/mids/pdsp/api/resps.jsonnet b/cfg/layers/low/resps.jsonnet similarity index 100% rename from cfg/layers/mids/pdsp/api/resps.jsonnet rename to cfg/layers/low/resps.jsonnet diff --git a/cfg/layers/low/ssss.jsonnet b/cfg/layers/low/ssss.jsonnet new file mode 100644 index 000000000..d1381dcda --- /dev/null +++ b/cfg/layers/low/ssss.jsonnet @@ -0,0 +1,45 @@ +// The spdir and other tests make use of "ssss" (simple sim sig smeared splat something something) +// This helps generate the variant configs. +// +// Run wirecell-gen morse-* to find long/tran smearing numbers that match the +// extra spread the sigproc induces. CAVEAT: these results are indirectly tied +// to the field response used in that process. Defaults here are from PDSP. + +function(nominal, + long=[ + 2.691862363980221, + 2.6750200122535057, + 2.7137567141154055 + ], + tran=[ + 0.7377218875719689, + 0.7157764520393882, + 0.13980698710556544 + ]) +{ + local smeared = nominal + { + splat: super.splat + { + "smear_long": long, + "smear_tran": tran, + } + }, + + // Used for test-morse-pdsp + morse_nominal: nominal, + + // Used for test-ssss-pdsp + ssss_nominal: nominal, + ssss_smeared: smeared, + + // Used for test/scripts/spdir + spdir_lofr: smeared { + detector_data: super.detector_data { + fields: "%s-fields-lo.json.bz2" % nominal.detector_name + }, + }, + spdir_hifr: smeared { + detector: super.detector_data { + fields: "%s-fields-hi.json.bz2" % nominal.detector_name + }, + }, +} diff --git a/cfg/layers/low/util.jsonnet b/cfg/layers/low/util.jsonnet index 99e304ba2..d4ff3b554 100644 --- a/cfg/layers/low/util.jsonnet +++ b/cfg/layers/low/util.jsonnet @@ -8,4 +8,11 @@ driftsToXregions :: function(vols) std.set(std.prune(std.flattenArrays([v.faces for v in vols])), function(o) std.toString(o)), + + assure_name :: function(name, obj) + if std.type(name) == "null" + then $.idents(obj) + else name, + + } diff --git a/cfg/layers/midapi.jsonnet b/cfg/layers/midapi.jsonnet index 4ed0a7450..5b111a0f0 100644 --- a/cfg/layers/midapi.jsonnet +++ b/cfg/layers/midapi.jsonnet @@ -1,21 +1,41 @@ -// This file defines the "mid" layer API. +// This defines and enforces the API that every detector must provide. // -// It provides a "base class" with dummy methods that will assert if not -// overridden. +// It "works" merely by forwarding calls to the a detector implementation (imp) +// in order to enforce the arguments. For methods that return a DFP subgraph +// the comments def expectation for input/output edges. // -// Caution: due to laziness it is likely woefully incomplete. +// This file should not be modified except to extend support to novel subgraph +// types. // -// Each detector variant mid API object will override what is here. -// Thus, this object simply gives a place to define an "abstract base -// class". Any method not provided by a concrete implementation and -// which gets called will result in an error. +// See mids/base/api.jsonnet for a default implementation that also gives more +// guidance on what detector implementations must provide. +local util = import "low/util.jsonnet"; + +function (imp) + local call_it(meth)={ + [meth]: function(anode, name=null) + imp[meth](anode, name=util.assure_name(name, anode)) + }; + local call_it_tags(meth)={ + [meth]: function(anode, name=null, tags=[]) + imp[meth](anode, name=util.assure_name(name, anode), tags=tags) + }; { - // Return a drifter configured for a detector. For optimal - // runtime performance, the implementation should return a - // DepoSetDrifter and NOT a chain with singular Drifter followed - // by Bagger. Mid implementer should use low.drifter - drifter() :: std.assertEqual("drifter" == "implemented"), - + // Return list of anode configuration objects. + anodes() : imp.anodes(), + + // Input DepoSet (original depos) output DepoSet (drifted depos). + drifter(name="") : imp.drifter(name), } ++ call_it("signal") ++ call_it("noise") ++ call_it("digitizer") ++ call_it("splat") ++ call_it("nf") ++ call_it("sp") ++ call_it_tags("reframer") ++ call_it("img") + + diff --git a/cfg/layers/mids.jsonnet b/cfg/layers/mids.jsonnet index 90ca01a4e..5473e8e83 100644 --- a/cfg/layers/mids.jsonnet +++ b/cfg/layers/mids.jsonnet @@ -1,7 +1,17 @@ // This file is generated, do not edit. -// Generated on Fri Dec 8 11:53:11 AM EST 2023 by bv on haiku.home +// Generated on Fri Mar 1 03:29:06 PM EST 2024 by bv on haiku.home // with /home/bv/wrk/wct/spdir-metric/toolkit/cfg/layers/gen-mids.sh { - pdsp : import "mids/pdsp/mids.jsonnet", - uboone : import "mids/uboone/mids.jsonnet", + base : { + api: import "mids/base/api.jsonnet", + variants: import "mids/base/variants.jsonnet", + }, + pdsp : { + api: import "mids/pdsp/api.jsonnet", + variants: import "mids/pdsp/variants.jsonnet", + }, + uboone : { + api: import "mids/uboone/api.jsonnet", + variants: import "mids/uboone/variants.jsonnet", + }, } diff --git a/cfg/layers/mids/base/api.jsonnet b/cfg/layers/mids/base/api.jsonnet new file mode 100644 index 000000000..b08778fee --- /dev/null +++ b/cfg/layers/mids/base/api.jsonnet @@ -0,0 +1,81 @@ +// This file defines the base or default "mid" layer API implementation. +// +// People looking to configure specific detectors should read it but very likely +// should NOT modify it. The comments give expectation for each subgraph. In +// particular Frames passed between subgraphs HAVE NO TRACE TAGS. +// +// This is called through midapi.jsonnet to assure interface correctness. +// +// API methods with no implementation will ASSERT and thus must be provided by a +// detector implementation. + +local low = import "../../low.jsonnet"; +local tn = low.wc.tn; +local pg = low.pg; + +function(services, params, options={}) +{ + // Return a list of all anode config objects. + anodes() : + low.anodes(params.geometry.volumes, params.geometry.wires_file), + + + // A node that returns a depo set drifter. + drifter(name="") : + low.drifter(services.random, + low.util.driftsToXregions(params.geometry.volumes), + params.lar, name=name), + + // DepoSet to voltage Frame. NO TRACE TAGS. + signal(anode, name) : + std.assertEqual("signal" == "implemented"), + + // Input voltage Frame, output noisy voltage Frame. NO TRACE TAGS. + noise(anode, name) : + std.assertEqual("signal" == "implemented"), + + // Voltage Frame to ADC Frame. NO TRACE TAGS + digitizer(anode, name) : + pg.pnode({ + type: 'Digitizer', + name: name, + data: params.digi { + anode: tn(anode), + } + }, nin=1, nout=1, uses=[anode]), + + // Pull out the responses as a short hand variable. + local sim_res = low.resps(params).sim, + + // DepoSet to signal Frame (an ersatz sim+SP). NO TRACE TAGS. + splat(anode, name) : + pg.pnode({ + type: 'DepoFluxSplat', + name: name, + data: params.splat { + anode: tn(anode), + field_response: tn(sim_res.fr), // note, only uses (period, origin) + }, + }, nin=1, nout=1, uses=[anode, sim_res.fr]), + + // noisy ADC frame to quiet ADC frame. NO TRACE TAGS. + nf(anode, name) : + std.assertEqual("nf" == "implemented"), + + // ADC frame to signal frame. NO TRACE TAGS. + sp(anode, name) : + std.assertEqual("sp" == "implemented"), + + // Signal frame to better signal Frame. NO TRACE TAGS + dnnroi(anode, name) : + std.assertEqual("dnnroi" == "implemented"), + + // Frame to Frame + reframer(anode, name, tags=[]) : + low.reframer(params, anode, name=name, tags=tags), + + // frame to blob set + img(anode, name) : + std.assertEqual("img" == "implemented"), + +} diff --git a/cfg/layers/mids/base/variants.jsonnet b/cfg/layers/mids/base/variants.jsonnet new file mode 100644 index 000000000..163788942 --- /dev/null +++ b/cfg/layers/mids/base/variants.jsonnet @@ -0,0 +1,212 @@ +// The base/variants.jsonnet provides base for "nominal" and "actual" parameters. +// +// To define a detector variant param object through INHERITANCE do like: +// +// // In /variants.jsonnet +// local base = import "../../base/variants.jsonnet"; +// local nominal = import "variants/nominal.jsonnet"; +// local actual = import "variants/actual.jsonnet"; +// { +// nominal: base.nominal + nominal, +// actual: self.nominal + actual, +// } +// +// Or to use replacement semantics, change the body to: +// +// { +// nominal: std.mergePatch(base.nominal, nominal), +// actual: std.mergePatch(self.nominal, actual), +// } +// +// That is, /variants/nominal.jsonnet overrides some parts of +// base.nominal and then further /variants/actual.jsonnet overrides +// some parts of the detector's nominal. + +local low = import "../../low.jsonnet"; +local wc = low.wc; +local detectors = import "detectors.jsonnet"; + +local nominal = { + + // Note, when constructed via high.params() these may be overridden before + // the resulting params object is passed to the detector mid api. + detector_name: "base", + variant_name: "nominal", + structure_name: "nominal", + + // The high level object summarizing the contents of wire-cell-data + // corresponding to a canonical detector name. + detector_data: detectors[self.detector_name], + + // The nominal base uses the same sample time domain binning everywhere. + // It provides the default binning in more local contexts (such as sim/sp). + binning : { + tick: 500*wc.ns, + nticks: 4096, // 2ms is enough for anyone! + }, + + // Nominal LAr properties. Almost certainly a detector should override + // some or all of these. + lar : { + // Longitudinal diffusion constant + DL : 4.0 * wc.cm2/wc.s, + // Transverse diffusion constant + DT : 8.8 * wc.cm2/wc.s, + // Electron lifetime + lifetime : 10*wc.ms, + // Electron drift speed, assumes a certain applied E-field + drift_speed : 1.6*wc.mm/wc.us, // approx at 500 V/cm + // LAr density + density: 1.389*wc.g/wc.centimeter3, + // Decay rate per mass for natural Ar39. + ar39activity: 1*wc.Bq/wc.kg, + }, + + + // This substructure describes the detector geometry. Almost certainly the + // geometry.volumes must be supplied by a specific detector. + // + // It's main chunk is the "volumes" substructure which must be consistent + // with the location of the wires in the "wires file". This can be + // challenging to get correct by hand. It can partly be derived from the + // wires file but that is too slow to do every load into WCT so it is + // recommended to generate it once like: + // + // $ wirecell-util wires-volumes | jsonnetfmt - > geometry.jsonnet + // + // The "anode", "response" and "cathode" planes can be given on that command + // line but here we override it to match other values from the config. + geometry_data: import "variants/geometry.jsonnet", + geometry : std.mergePatch($.geometry_data, { + xplanes: { + // Distances between collection plane and a logical plane. + danode: 0.0, + dresponse: 10*wc.cm, + dcathode: nominal.binning.nticks * nominal.binning.tick / nominal.lar.drift_speed, + } + }), + + + // The electronics response parameters. Almost certainly, the inheriting + // detector should override this section. + elec : { + // The type of electronics + type: "cold", + + // The FE amplifier gain in units of Voltage/Charge. + gain : 14.0*wc.mV/wc.fC, + + // The shaping (aka peaking) time of the amplifier shaper. + shaping : 2.0*wc.us, + + // A realtive gain applied after shaping and before ADC. + // + // Don't use this to fix wrong sign depos. If you neeed to + // fix sign of eg larsoft depos its input component has a + // "scale" parameter which can be given a -1. Or better, fix + // the C++ that sets the wrong sign to begin with. + postgain: 1.0, + }, + + // The "RC" response. The detector should override and some may require + // DFP structural changes to implement RCRC + rc : { + width: 1.0*wc.ms, + }, + + // The ductor transforms drifted depos to currents + ductor: { + + // The (single) field file. + field_file: $.detector_data.field, + + // Note, an extension to support "actual" detector may require multiple + // fields, as in the case of uboone. The above variable should remain + // singular. Access multiple files through another variable such as: + // field_files: $.detector_data.fields, // (plural) + + // Nominally, we assume simulation time starts at t=0. However, any + // ionization deposited between "response" and "anode" planes are + // "backed up" in space and time. In order to not trim them out we must + // start the simulation earlier so that a depo at t=0 at the "anode" + // plane will have its response in tick >= 0. + start_time : -($.geometry.xplanes.dresponse - $.geometry.xplanes.danode)/$.lar.drift_speed, + + // Most simply, no custom binning. + tbin: 0, + binning: $.binning, + readout_time : self.binning.nticks * self.binning.tick, + }, + + // Simulating noise + noise : { + model: { + spectra_file: $.detector_data.noise, + + // These are frequency space binning which are not necessarily + // same as some time binning - but here they are. + period: $.binning.tick, // 1/frequency + nsamples: $.binning.nticks, // + // Optimize binning + wire_length_scale: 1.0*wc.cm, + }, + // Optimize use of randoms + replacement_percentage: 0.02, + }, + + // Digitization model parameters. + digi : { + // A relative gain applied just prior to digitization. This + // is not FE gain, see elec for that. + gain: 1.0, + + // Voltage baselines added to any input voltage signal listed + // in a per plan (U,V,W) array. + baselines: [1*wc.volt, 1*wc.volt, 0.5*wc.volt], + + // The resolution (bits) of the ADC + resolution: 12, + + // The voltage range as [min,max] of the ADC, eg min voltage + // counts 0 ADC, max counts 2^resolution-1. + fullscale: [0*wc.volt, 2*wc.volt], + }, + + // Configure the "splat" (DepoFluxSplat) is an approximation to the + // combination of ductor+sigproc. + splat : { + sparse: true, + tick: $.binning.tick, + window_start: 0, + window_duration: self.tick * $.binning.nticks, + reference_time: 0.0, + }, + + // The noise filter parameter pack + nf : { + // In principle, may use a different field response. + field_file : $.ductor.field_file, + + binning: $.binning, + + // The name of the ChannelNoiseDb + chndb: "perfect", + filters: { + channel: ["single"], // sticky, single, gaincalib + grouped: [], // grouped, + status: [], + }, + }, + + // Signal processing parameter pack + sp: { + // In principle, may use a different field response. + field_file : $.ductor.field_file, + binning : $.binning, + }, +}; + +{ + nominal: nominal, + actual: nominal, +} + low.ssss(nominal) diff --git a/cfg/layers/mids/base/variants/geometry.jsonnet b/cfg/layers/mids/base/variants/geometry.jsonnet new file mode 100644 index 000000000..90413e1e0 --- /dev/null +++ b/cfg/layers/mids/base/variants/geometry.jsonnet @@ -0,0 +1,155 @@ +local wc = import 'wirecell.jsonnet'; +{ + + wires_file: 'protodune-wires-larsoft-v4.json.bz2', + + // distance between collection wire plane and a plane. + xplanes: { + danode: 10.0 * wc.mm, + dresponse: 100.0 * wc.mm, + dcathode: 1000.0 * wc.mm, + }, + local xplanes = self.xplanes, // to make available below + + volumes: [ + + { + wires: 0, + xcenter: -3636.9 * wc.mm, // absolute center of APA + local xcenter = self.xcenter, // to make available below. + faces: [ + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter + (xplanes.danode + dcollection), + response: xcenter + (xplanes.dresponse + dcollection), + cathode: xcenter + (xplanes.dcathode + dcollection), + }, + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter - (xplanes.danode + dcollection), + response: xcenter - (xplanes.dresponse + dcollection), + cathode: xcenter - (xplanes.dcathode + dcollection), + }, + + ], + }, + + { + wires: 1, + xcenter: 3636.9 * wc.mm, // absolute center of APA + local xcenter = self.xcenter, // to make available below. + faces: [ + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter + (xplanes.danode + dcollection), + response: xcenter + (xplanes.dresponse + dcollection), + cathode: xcenter + (xplanes.dcathode + dcollection), + }, + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter - (xplanes.danode + dcollection), + response: xcenter - (xplanes.dresponse + dcollection), + cathode: xcenter - (xplanes.dcathode + dcollection), + }, + + ], + }, + + { + wires: 2, + xcenter: -3636.9 * wc.mm, // absolute center of APA + local xcenter = self.xcenter, // to make available below. + faces: [ + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter + (xplanes.danode + dcollection), + response: xcenter + (xplanes.dresponse + dcollection), + cathode: xcenter + (xplanes.dcathode + dcollection), + }, + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter - (xplanes.danode + dcollection), + response: xcenter - (xplanes.dresponse + dcollection), + cathode: xcenter - (xplanes.dcathode + dcollection), + }, + + ], + }, + + { + wires: 3, + xcenter: 3636.9 * wc.mm, // absolute center of APA + local xcenter = self.xcenter, // to make available below. + faces: [ + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter + (xplanes.danode + dcollection), + response: xcenter + (xplanes.dresponse + dcollection), + cathode: xcenter + (xplanes.dcathode + dcollection), + }, + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter - (xplanes.danode + dcollection), + response: xcenter - (xplanes.dresponse + dcollection), + cathode: xcenter - (xplanes.dcathode + dcollection), + }, + + ], + }, + + { + wires: 4, + xcenter: -3636.9 * wc.mm, // absolute center of APA + local xcenter = self.xcenter, // to make available below. + faces: [ + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter + (xplanes.danode + dcollection), + response: xcenter + (xplanes.dresponse + dcollection), + cathode: xcenter + (xplanes.dcathode + dcollection), + }, + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter - (xplanes.danode + dcollection), + response: xcenter - (xplanes.dresponse + dcollection), + cathode: xcenter - (xplanes.dcathode + dcollection), + }, + + ], + }, + + { + wires: 5, + xcenter: 3636.9 * wc.mm, // absolute center of APA + local xcenter = self.xcenter, // to make available below. + faces: [ + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter + (xplanes.danode + dcollection), + response: xcenter + (xplanes.dresponse + dcollection), + cathode: xcenter + (xplanes.dcathode + dcollection), + }, + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter - (xplanes.danode + dcollection), + response: xcenter - (xplanes.dresponse + dcollection), + cathode: xcenter - (xplanes.dcathode + dcollection), + }, + + ], + }, + + ], +} diff --git a/cfg/layers/mids/pdsp/api.jsonnet b/cfg/layers/mids/pdsp/api.jsonnet index 23cfbabb3..f5e9dddee 100644 --- a/cfg/layers/mids/pdsp/api.jsonnet +++ b/cfg/layers/mids/pdsp/api.jsonnet @@ -1,28 +1,13 @@ // return a mid-level API object -local low = import "../../low.jsonnet"; +local nf = import "api/nf.jsonnet"; +local sp = import "api/sp.jsonnet"; local sim = import "api/sim.jsonnet"; -local sigproc = import "api/sigproc.jsonnet"; local img = import "api/img.jsonnet"; // Create a mid API object. No options supported. function(services, params, options={}) { - - anodes :: function() - low.anodes(params.geometry.drifts, params.geometry.wires_file), - - drifter :: function(name="") - low.drifter(services.random, - low.util.driftsToXregions(params.geometry.drifts), - params.lar, name=name), - - // track_depos, signal, noise, digitizer - sim :: sim(services, params), - - // nf, sp, dnnroi - sigproc :: sigproc(services, params, options), - - img :: img(services, params), -} - - + nf : nf(services, params, options), + sp : sp(services, params, options), + img : img(services, params) +} + sim(services, params) diff --git a/cfg/layers/mids/pdsp/api/img.jsonnet b/cfg/layers/mids/pdsp/api/img.jsonnet index 54e1802ae..609529d88 100644 --- a/cfg/layers/mids/pdsp/api/img.jsonnet +++ b/cfg/layers/mids/pdsp/api/img.jsonnet @@ -4,13 +4,12 @@ local low = import "../../../low.jsonnet"; local pg = low.pg; local wc = low.wc; -function(services, params) function(anode) - local ident = low.util.idents(anode); - local img = low.img(anode); +function(services, params) function(anode, name) + local img = low.img(anode, name); local wfm = { type: 'WaveformMap', - name: ident, + name: name, data: { filename: params.img.charge_error_file, }, @@ -20,11 +19,11 @@ function(services, params) function(anode) active_planes=[0,1,2], masked_planes=[], dummy_planes=[]) = pg.pnode({ type: "MaskSlices", - name: ident+ext, + name: name+ext, data: { - wiener_tag: "wiener"+ident, - charge_tag: "gauss"+ident, - error_tag: "gauss_error"+ident, + wiener_tag: "wiener", + charge_tag: "gauss", + error_tag: "gauss_error", tick_span: span, anode: wc.tn(anode), min_tbin: min_tbin, @@ -40,10 +39,10 @@ function(services, params) function(anode) local sigunc = pg.pnode({ type: 'ChargeErrorFrameEstimator', - name: ident, + name: name, data: { - intag: "gauss"+ident, - outtag: 'gauss_error'+ident, + intag: "gauss", + outtag: 'gauss_error', anode: wc.tn(anode), rebin: 4, // this number should be consistent with the waveform_map choice fudge_factors: [2.31, 2.31, 1.1], // fudge factors for each plane [0,1,2] @@ -59,6 +58,6 @@ function(services, params) function(anode) img.tiling(), img.clustering(), img.grouping(), - img.charge_solving(), // fixme: a few optins we may want to allow to specify in variant params - ], ident) + img.charge_solving(), + ], name=name) diff --git a/cfg/layers/mids/pdsp/api/nf.jsonnet b/cfg/layers/mids/pdsp/api/nf.jsonnet index f239404dc..6d6fd853b 100644 --- a/cfg/layers/mids/pdsp/api/nf.jsonnet +++ b/cfg/layers/mids/pdsp/api/nf.jsonnet @@ -2,6 +2,7 @@ local low = import "../../../low.jsonnet"; local tn = low.wc.tn; +local get = low.wc.get; local chndbs = { "actual": import "chndb-actual.jsonnet", @@ -10,16 +11,11 @@ local chndbs = { local filters = import "noise-filters.jsonnet"; local frs = import "frs.jsonnet"; -function(services, params, options) function(anode) - local pars = std.mergePatch(params, std.get(options, "params", {})); - - local ident = low.util.idents(anode); - - // nf uses same fr as sp - local fr = frs(pars).nf; +function(services, params, options) function(anode, name) + local pars = std.mergePatch(params, get(options, "params", {})); local chndb = chndbs[pars.nf.chndb](anode, pars.nf.binning, - fr, services.dft); + frs(pars).nf, services.dft); // filter look up local flu = filters(services, pars, anode, chndb); @@ -34,7 +30,7 @@ function(services, params, options) function(anode) // the NF low.pg.pnode({ type: 'OmnibusNoiseFilter', - name: ident, + name: name, data: { // Nonzero forces the number of ticks in the waveform @@ -47,8 +43,10 @@ function(services, params, options) function(anode) grouped_filters: [tn(obj) for obj in fobj.grouped], channel_status_filters: [tn(obj) for obj in fobj.status], noisedb: tn(chndb), - intraces: 'orig' + ident, - outtraces: 'raw' + ident, + intraces: '', + outtraces: '', + // intraces: 'orig' + ident, + // outtraces: 'raw' + ident, }, }, nin=1, nout=1, uses = fobj.channel + fobj.grouped + fobj.status) diff --git a/cfg/layers/mids/pdsp/api/sigproc.jsonnet b/cfg/layers/mids/pdsp/api/sigproc.jsonnet deleted file mode 100644 index d1d433f8f..000000000 --- a/cfg/layers/mids/pdsp/api/sigproc.jsonnet +++ /dev/null @@ -1,21 +0,0 @@ -local low = import "../../../low.jsonnet"; - -// These are big and unique to PDSP so we have to write them out -// exhaustively. -local nf = import "nf.jsonnet"; -local sp = import "sp.jsonnet"; - -function(services, params, options={}) { - - // API method - nf :: nf(services, params, options), - - // API method - sp :: sp(services, params, options), - - // API method - dnnroi :: function(anode) - low.dnnroi(anode, params.dnnroi.model_file, - services.platform, params.dnnroi.output_scale), - -} diff --git a/cfg/layers/mids/pdsp/api/sim.jsonnet b/cfg/layers/mids/pdsp/api/sim.jsonnet index 2c1ff892f..952eee263 100644 --- a/cfg/layers/mids/pdsp/api/sim.jsonnet +++ b/cfg/layers/mids/pdsp/api/sim.jsonnet @@ -11,14 +11,9 @@ local wc = low.wc; local pg = low.pg; local idents = low.util.idents; -local resps = import "resps.jsonnet"; - function(services, params, options={}) { - // Signal binning may be extended from nominal. - local sig_binning = params.ductor.binning, - - local res = resps(params).sim, + local res = low.resps(params).sim, // some have more than one local short_responses = [ res.er, ], @@ -28,7 +23,7 @@ function(services, params, options={}) { type: 'PlaneImpactResponse', name: std.toString(plane), uses: [res.fr, services.dft] + short_responses + long_responses, - data: sig_binning { + data: params.ductor.binning { plane: plane, dft: wc.tn(services.dft), field_response: wc.tn(res.fr), @@ -41,7 +36,7 @@ function(services, params, options={}) { // API method sim.track_depos: subgraph making some depos from // ideal tracks in the detector. - track_depos :: function(tracklist = [{ + track_depos : function(tracklist = [{ time: 0, charge: -5000, ray: { @@ -54,28 +49,14 @@ function(services, params, options={}) { params.ductor.start_time) ]), - // API method sim.reframer - reframer :: function(anode, tags=[], name=null) - pg.pnode({ - type: 'Reframer', - name: if std.type(name) == "null" then idents(anode) else name, - data: { - anode: wc.tn(anode), - tags: tags, - fill: 0.0, - toffset: 0, - tbin: params.ductor.tbin, - nticks: sig_binning.nticks - self.tbin, - }, - }, nin=1, nout=1, uses=[anode]), // API method sim.signal: subgraph making pure signal voltage from // depos. - signal :: function(anode, tags=[]) + signal : function(anode, name=null) pg.pipeline([ pg.pnode({ type: 'DepoTransform', - name: idents(anode), + name: name, data: { rng: wc.tn(services.random), dft: wc.tn(services.dft), @@ -86,17 +67,17 @@ function(services, params, options={}) { first_frame_number: 0, readout_time: params.ductor.readout_time, start_time: params.ductor.start_time, - tick: sig_binning.tick, + tick: params.ductor.binning.tick, nsigma: 3, }, }, nin=1, nout=1, uses = pirs + [anode, services.random, services.dft]), - $.reframer(anode, tags=tags, name=idents(anode))]), + low.reframer(params, anode, name=name)]), // API method sim.noise: subgraph adding noise to voltage - noise :: function(anode) + noise : function(anode, name) local model = { type: 'EmpiricalNoiseModel', - name: idents(anode), + name: name, data: params.noise.model { anode: wc.tn(anode), dft: wc.tn(services.dft), @@ -106,7 +87,7 @@ function(services, params, options={}) { }; pg.pnode({ type: 'AddNoise', - name: idents(anode), + name: name, data: { rng: wc.tn(services.random), dft: wc.tn(services.dft), @@ -115,57 +96,5 @@ function(services, params, options={}) { replacement_percentage: params.noise.replacement_percentage, }}, nin=1, nout=1, uses=[services.random, services.dft, model]), - // API method sim.digitizer: return subgraph adding digitization - // of voltage to produce ADC - digitizer :: function(anode) - pg.pnode({ - type: 'Digitizer', - name: idents(anode), - data: params.digi { - anode: wc.tn(anode), - frame_tag: "orig" + idents(anode), - } - }, nin=1, nout=1, uses=[anode]), - - // The approximated sim+sigproc - splat :: function(anode, name=null) - pg.pnode({ - type: 'DepoFluxSplat', - name: if std.type(name) == "null" then idents(anode) else name, - data: params.splat { - anode: wc.tn(anode), - field_response: wc.tn(res.fr), - }, - }, nin=1, nout=1, uses=[anode, res.fr]), - - // Construct obsolete single-depo, not generic "splat" - singledeposplat :: function(anode) - pg.pnode({ - type: 'DuctorFramer', - name: idents(anode), - data: { - ductor: 'DepoSplat:' + idents(anode), - fanin: 'FrameFanin:' + idents(anode), - }, - }, nin=1, nout=1, uses=[ { - type: 'FrameFanin', - name: idents(anode), - data: { - multiplicity: 0, // "dynamic" - } - }, { - type: 'DepoSplat', - name: idents(anode), - data: { - anode: wc.tn(anode), - nsigma: 3.0, - start_time: params.ductor.start_time, - readout_time: params.ductor.readout_time, - tick: params.ductor.binning.tick, - continuous: true, // see DepoSplat comments - fixed: false, - drift_speed: params.lar.drift_speed, - } }]), - } diff --git a/cfg/layers/mids/pdsp/api/sp.jsonnet b/cfg/layers/mids/pdsp/api/sp.jsonnet index 71f9773bf..dd4f089f2 100644 --- a/cfg/layers/mids/pdsp/api/sp.jsonnet +++ b/cfg/layers/mids/pdsp/api/sp.jsonnet @@ -7,23 +7,20 @@ local spfilt = import "sp-filters.jsonnet"; local low = import "../../../low.jsonnet"; local wc = low.wc; -local resps = import "resps.jsonnet"; - // Allow an optional argument "sparse" as this is really an end-user // decision. Higher layers may expose this option to the TLA. -function(services, params, options={}) function(anode) +function(services, params, options={}) function(anode, name) local pars = std.mergePatch(params, std.get(options, "params", {})); local opts = {sparse:true} + options; - local ident = low.util.idents(anode); local resolution = pars.digi.resolution; local fullscale = pars.digi.fullscale[1] - pars.digi.fullscale[0]; local ADC_mV_ratio = ((1 << resolution) - 1 ) / fullscale; - local res = resps(params).sp; + local res = low.resps(params).sp; low.pg.pnode({ type: 'OmnibusSigProc', - name: ident, + name: name, data: { /** @@ -56,23 +53,25 @@ function(services, params, options={}) function(anode) r_sep_peak: 6.0, // default 6.0 r_low_peak_sep_threshold_pre: 1200, // default 1200 - // frame tags - wiener_tag: 'wiener' + ident, - decon_charge_tag: 'decon_charge' + ident, - gauss_tag: 'gauss' + ident, + /// These defaults are okay, even with out anode ident qualifier + // wiener_tag: 'wiener', + // gauss_tag: 'gauss', + // FIXME: we should not need to "unset" these. They should empty by + // default and use_roi_debug_mode should be removed. use_roi_debug_mode: false, - tight_lf_tag: 'tight_lf' + ident, - loose_lf_tag: 'loose_lf' + ident, - cleanup_roi_tag: 'cleanup_roi' + ident, - break_roi_loop1_tag: 'break_roi_1st' + ident, - break_roi_loop2_tag: 'break_roi_2nd' + ident, - shrink_roi_tag: 'shrink_roi' + ident, - extend_roi_tag: 'extend_roi' + ident, + decon_charge_tag: '', + tight_lf_tag: '', + loose_lf_tag: '', + cleanup_roi_tag: '', + break_roi_loop1_tag: '', + break_roi_loop2_tag: '', + shrink_roi_tag: '', + extend_roi_tag: '', use_multi_plane_protection: false, - mp3_roi_tag: 'mp3_roi' + ident, - mp2_roi_tag: 'mp2_roi' + ident, + mp3_roi_tag: '', + mp2_roi_tag: '', isWrapped: false, sparse : opts.sparse, diff --git a/cfg/layers/mids/pdsp/mids.jsonnet b/cfg/layers/mids/pdsp/mids.jsonnet index ce660be0d..790343ddd 100644 --- a/cfg/layers/mids/pdsp/mids.jsonnet +++ b/cfg/layers/mids/pdsp/mids.jsonnet @@ -3,11 +3,7 @@ // The "real" code is elsewhere local api = import "api.jsonnet"; - -// Import all variant parameter packs to allow for dict like lookup. -local variants = { - nominal: import "variants/nominal.jsonnet", -} + import "variants/ssss.jsonnet"; +local variants = import "variants.jsonnet"; function(services, variant="", options={}) api(services, variants[variant], options) diff --git a/cfg/layers/mids/pdsp/variants.jsonnet b/cfg/layers/mids/pdsp/variants.jsonnet new file mode 100644 index 000000000..31392c490 --- /dev/null +++ b/cfg/layers/mids/pdsp/variants.jsonnet @@ -0,0 +1,11 @@ +// Import all variant parameter packs to allow for dict like lookup. + +local ssss = import "variants/ssss.jsonnet"; +local resample = import "variants/resample.jsonnet" ; +local nominal = import "variants/nominal.jsonnet"; +local actual = import "variants/actual.jsonnet"; + +{ + nominal: nominal, + actual: actual +} + ssss + resample diff --git a/cfg/layers/mids/pdsp/variants/actual.jsonnet b/cfg/layers/mids/pdsp/variants/actual.jsonnet new file mode 100644 index 000000000..434a764ef --- /dev/null +++ b/cfg/layers/mids/pdsp/variants/actual.jsonnet @@ -0,0 +1,44 @@ +local nominal = import "nominal.jsonnet"; +nominal { + + // The ductor transforms drifted depos to currents + ductor: { + + # field_file: "dune-garfield-1d565.json.bz2", + field_file: $.detector.fields, + + // The distance from the anode centerline to where the field + // response calculation begins drifting. Take care that field + // response calculations are not always done with the + // thickness of the anode in mind and may report their + // "response plane" measured relative to a number of other + // monuments. This value is almost certainly required in the + // "geometry" section. + response_plane : 10*wc.cm, + + // The "absolute" time (ie, in G4 time) that the lower edge of + // of final readout tick #0 should correspond to. This is a + // "fixed" notion. + local tzero = -250*wc.us, + + // Simulate extra time so we can see the tail of response + // from activity from depos that fall outside the readout + // window. This amount of extra is somewhat arbitrary but + // the duration of the field response calculation is a + // good measure. + local pre = self.response_plane / $.lar.drift_speed, + + // The time bin where the readout should be considered to + // actually start given the pre-signal simulated. + tbin : wc.roundToInt(pre / self.binning.tick), + binning: { + // Ductor ticks at same speed as ADC + tick : $.binning.tick, + // Over the somewhat enlarged domain + nticks : $.ductor.tbin + $.binning.nticks, + }, + start_time : tzero - pre, + readout_time : self.binning.nticks * self.binning.tick, + }, + +} diff --git a/cfg/layers/mids/pdsp/variants/geometry.jsonnet b/cfg/layers/mids/pdsp/variants/geometry.jsonnet new file mode 100644 index 000000000..90413e1e0 --- /dev/null +++ b/cfg/layers/mids/pdsp/variants/geometry.jsonnet @@ -0,0 +1,155 @@ +local wc = import 'wirecell.jsonnet'; +{ + + wires_file: 'protodune-wires-larsoft-v4.json.bz2', + + // distance between collection wire plane and a plane. + xplanes: { + danode: 10.0 * wc.mm, + dresponse: 100.0 * wc.mm, + dcathode: 1000.0 * wc.mm, + }, + local xplanes = self.xplanes, // to make available below + + volumes: [ + + { + wires: 0, + xcenter: -3636.9 * wc.mm, // absolute center of APA + local xcenter = self.xcenter, // to make available below. + faces: [ + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter + (xplanes.danode + dcollection), + response: xcenter + (xplanes.dresponse + dcollection), + cathode: xcenter + (xplanes.dcathode + dcollection), + }, + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter - (xplanes.danode + dcollection), + response: xcenter - (xplanes.dresponse + dcollection), + cathode: xcenter - (xplanes.dcathode + dcollection), + }, + + ], + }, + + { + wires: 1, + xcenter: 3636.9 * wc.mm, // absolute center of APA + local xcenter = self.xcenter, // to make available below. + faces: [ + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter + (xplanes.danode + dcollection), + response: xcenter + (xplanes.dresponse + dcollection), + cathode: xcenter + (xplanes.dcathode + dcollection), + }, + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter - (xplanes.danode + dcollection), + response: xcenter - (xplanes.dresponse + dcollection), + cathode: xcenter - (xplanes.dcathode + dcollection), + }, + + ], + }, + + { + wires: 2, + xcenter: -3636.9 * wc.mm, // absolute center of APA + local xcenter = self.xcenter, // to make available below. + faces: [ + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter + (xplanes.danode + dcollection), + response: xcenter + (xplanes.dresponse + dcollection), + cathode: xcenter + (xplanes.dcathode + dcollection), + }, + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter - (xplanes.danode + dcollection), + response: xcenter - (xplanes.dresponse + dcollection), + cathode: xcenter - (xplanes.dcathode + dcollection), + }, + + ], + }, + + { + wires: 3, + xcenter: 3636.9 * wc.mm, // absolute center of APA + local xcenter = self.xcenter, // to make available below. + faces: [ + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter + (xplanes.danode + dcollection), + response: xcenter + (xplanes.dresponse + dcollection), + cathode: xcenter + (xplanes.dcathode + dcollection), + }, + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter - (xplanes.danode + dcollection), + response: xcenter - (xplanes.dresponse + dcollection), + cathode: xcenter - (xplanes.dcathode + dcollection), + }, + + ], + }, + + { + wires: 4, + xcenter: -3636.9 * wc.mm, // absolute center of APA + local xcenter = self.xcenter, // to make available below. + faces: [ + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter + (xplanes.danode + dcollection), + response: xcenter + (xplanes.dresponse + dcollection), + cathode: xcenter + (xplanes.dcathode + dcollection), + }, + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter - (xplanes.danode + dcollection), + response: xcenter - (xplanes.dresponse + dcollection), + cathode: xcenter - (xplanes.dcathode + dcollection), + }, + + ], + }, + + { + wires: 5, + xcenter: 3636.9 * wc.mm, // absolute center of APA + local xcenter = self.xcenter, // to make available below. + faces: [ + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter + (xplanes.danode + dcollection), + response: xcenter + (xplanes.dresponse + dcollection), + cathode: xcenter + (xplanes.dcathode + dcollection), + }, + + { + local dcollection = 52.3 * wc.mm, + anode: xcenter - (xplanes.danode + dcollection), + response: xcenter - (xplanes.dresponse + dcollection), + cathode: xcenter - (xplanes.dcathode + dcollection), + }, + + ], + }, + + ], +} diff --git a/cfg/layers/mids/pdsp/variants/nominal.jsonnet b/cfg/layers/mids/pdsp/variants/nominal.jsonnet index 67cddac98..e3da9ae7d 100644 --- a/cfg/layers/mids/pdsp/variants/nominal.jsonnet +++ b/cfg/layers/mids/pdsp/variants/nominal.jsonnet @@ -19,132 +19,19 @@ // 5. Add this new object to the variants object in mids.jsonnet. local wc = import "wirecell.jsonnet"; +local base_variants = import "../../base/variants.jsonnet"; -local detectors = import "detectors.jsonnet"; -local mydet = detectors.pdsp; +base_variants.nominal { -{ - // Define the LAr properties. In principle, this is NOT subject - // to variance though in principle it could be (eg, top/bottom of - // DUNE-VD may have different LAr temps?). We include this here - // to keep data out of the the mid code and because this info is - // used in a few places. - lar : { - // Longitudinal diffusion constant - DL : 4.0 * wc.cm2/wc.s, - // Transverse diffusion constant - DT : 8.8 * wc.cm2/wc.s, + lar : super.lar { // Electron lifetime lifetime : 10.4*wc.ms, // Electron drift speed, assumes a certain applied E-field drift_speed : 1.565*wc.mm/wc.us, // at 500 V/cm - // LAr density - density: 1.389*wc.g/wc.centimeter3, - // Decay rate per mass for natural Ar39. - ar39activity: 1*wc.Bq/wc.kg, }, - // Extract all the location for planes used in geometry - planes : { - apa_cpa: 3.63075*wc.m, // LArSoft - // apa_cpa = 3.637*wc.m, // DocDB 203 - cpa_thick: 3.175*wc.mm, // 1/8", from Bo Yu (BNL) and confirmed with LArSoft - // cpa_thick = 50.8*wc.mm, // DocDB 203 - apa_w2w: 85.725*wc.mm, // between W planes, DocDB 203 - apa_g2g: 114.3*wc.mm, - // measured from grid to "anode" plane along drift direction - offset: 4.76*wc.mm - }, - - // This function returns the minimal geometry needed by WCT. In - // the transverse plane the geometry is defined by the "wires - // file" and in the drift direction by the "drifts" structure - // which defines a trio of planes: "anode", "response" and - // "cathode" which are perpendicular to drift. Each are assumed - // approximately equipotential along their plane. If your - // detector fails this, then use a "variant params" pack for each - // detector region which does meet this requirement. - // - geometry : { - - // The exhaustive list of the location of every single wire - // (or strip) segment. - // wires_file: "protodune-wires-larsoft-v4.json.bz2", - wires_file: mydet.wires, - - // Comments on how to chose the "anode" plane location: The - // "anode" cut off plane, here measured from APA centerline, - // determines how close to the wires do we consider any depo. - // Anything closer will simply be discarded, else it will - // either be drifted or "backed up" to the response plane. - // This is somewhat arbitrary choice. The "backing up" assume - // parallel drift points and so some small error is made when - // applied to ionization produce very close to the wires. To - // reduce this small error the "anode" plane can be moved - // closer to the "response" plane but at the cost of not - // simulate the volume left beyond the "anode" plane. - - // local apa_plane = 0.5*apa_g2g, // put anode plane at grid plane - local apa_plane = 0.5*$.planes.apa_g2g - $.planes.offset, // put anode plane at first induction plane - - // Response plane MUST be located so that the wire plane - // locations used in the FR calculation matches the wire plane - // locations given to WCT. - local res_plane = 0.5*$.planes.apa_w2w + $.ductor.response_plane, - - // The cathode plane is like the anode cut off plane. Any - // depo not between the two is dropped prior to drifting. It - // is demarks the face of the CPA. - local cpa_plane = $.planes.apa_cpa - 0.5*$.planes.cpa_thick, - - // The drifts are then defined in terms of these above - // numbers. You can use "wirecell-util wires-info" or - // "wirecell-util wires-volumes" or others to understand the - // mapping of anode number to the 6 locations in global X and - // Z coordinates. For Larsoft wires the numbering is column - // major starting at small X and Z so the centerline is - // -/+/-/+/-/+. Also important is that the faces are listed - // "front" first. Front is the one with the more positive X - // coordinates and if we want to ignore a face it is made - // null. - drifts : [{ - local sign = 2*(n%2)-1, - local centerline = sign*$.planes.apa_cpa, - wires: n, // anode number - name: "apa%d"%n, - faces: - // top, front face is against cryo wall - if sign > 0 - then [ - { - anode: centerline + apa_plane, - response: centerline + res_plane, - cathode: centerline + cpa_plane, - }, - { - anode: centerline - apa_plane, - response: centerline - res_plane, - cathode: centerline - cpa_plane, - } - ] - // bottom, back face is against cryo wall - else [ - { - anode: centerline + apa_plane, - response: centerline + res_plane, - cathode: centerline + cpa_plane, - }, - { - anode: centerline - apa_plane, - response: centerline - res_plane, - cathode: centerline - cpa_plane, - } - - ], - } for n in std.range(0,5)], - }, + geometry_data : import "geometry.jsonnet", - // The nominal time binning for data produced by the detector. binning : { tick: 0.5*wc.us, nticks: 6000, @@ -175,73 +62,6 @@ local mydet = detectors.pdsp; width: 1.1*wc.ms, }, - // The ductor transforms drifted depos to currents - ductor: { - - # field_file: "dune-garfield-1d565.json.bz2", - field_file: mydet.fields, - - // The distance from the anode centerline to where the field - // response calculation begins drifting. Take care that field - // response calculations are not always done with the - // thickness of the anode in mind and may report their - // "response plane" measured relative to a number of other - // monuments. This value is almost certainly required in the - // "geometry" section. - response_plane : 10*wc.cm, - - // The "absolute" time (ie, in G4 time) that the lower edge of - // of final readout tick #0 should correspond to. This is a - // "fixed" notion. - local tzero = -250*wc.us, - - // Simulate extra time so we can see the tail of response - // from activity from depos that fall outside the readout - // window. This amount of extra is somewhat arbitrary but - // the duration of the field response calculation is a - // good measure. - local pre = self.response_plane / $.lar.drift_speed, - - // The time bin where the readout should be considered to - // actually start given the pre-signal simulated. - tbin : wc.roundToInt(pre / self.binning.tick), - binning: { - // Ductor ticks at same speed as ADC - tick : $.binning.tick, - // Over the somewhat enlarged domain - nticks : $.ductor.tbin + $.binning.nticks, - }, - start_time : tzero - pre, - readout_time : self.binning.nticks * self.binning.tick, - }, - - // A "splat" (DepoFluxSplat) is an approximation to the combination of - // ductor+sigproc - splat : { - sparse: true, - tick: $.ductor.binning.tick, - window_start: $.ductor.start_time, - window_duration: $.ductor.readout_time, - reference_time: 0.0, - - }, - - // Simulating noise - noise : { - model: { - #spectra_file: "protodune-noise-spectra-v1.json.bz2", - spectra_file: mydet.noise, - - // These are frequency space binning which are not necessarily - // same as some time binning - but here they are. - period: $.binning.tick, // 1/frequency - nsamples: $.binning.nticks, // - // Optimize binning - wire_length_scale: 1.0*wc.cm, - }, - // Optimize use of randoms - replacement_percentage: 0.02, - }, // Digitization model parameters. digi : { @@ -285,41 +105,6 @@ local mydet = detectors.pdsp; }, }, - // Signal processing parameter pack - sp: { - // In principle, may use a different field response. - field_file : $.ductor.field_file, - binning : $.binning, - }, - - dnnroi: { - // this is a total fudge factor - output_scale : 1.0, - - // unet-l23-cosmic500-e50.ts was made for uboone.... - model_file : "unet-l23-cosmic500-e50.ts", - }, - - - // Imaging paramter pack - img : { - // For now we use MicroBooNE's - #"charge_error_file": "microboone-charge-error.json.bz2", - "charge_error_file": mydet.qerr, - - // Number of ticks to collect into one time slice span - span: 4, - - binning : $.binning, - - // The tiling strategy refers to how to deal with substantial - // dead channels. The "perfect" strategy simply performs so - // called "3-plane" tiling. The "mixed" strategy also - // includes so called "2-plane" tiling, producing more and - // more ambiguous blobs but will not have "gaps" due to dead - // channels. - tiling_strategy: "perfect", - }, } diff --git a/cfg/layers/mids/pdsp/variants/resample.jsonnet b/cfg/layers/mids/pdsp/variants/resample.jsonnet new file mode 100644 index 000000000..d54be964f --- /dev/null +++ b/cfg/layers/mids/pdsp/variants/resample.jsonnet @@ -0,0 +1,26 @@ +local wc = import "wirecell.jsonnet"; +local nominal = import "nominal.jsonnet"; + +{ + resample_500ns : nominal + { + binning : { + tick: 500 * wc.ns, + nticks: 6000, + }, + + ductor: super.ductor { + field_file: "dune-garfield-1d565-100ns.json.bz2", + }, + }, + resample_512ns : nominal + { + binning : { + tick: 512 * wc.ns, + nticks: 6000 * 512 / 500, + }, + + ductor: super.ductor { + field_file: "dune-garfield-1d565-64ns.json.bz2", + }, + + }, +} diff --git a/cfg/layers/mids/pdsp/variants/simple.jsonnet b/cfg/layers/mids/pdsp/variants/simple.jsonnet deleted file mode 100644 index 8ca144938..000000000 --- a/cfg/layers/mids/pdsp/variants/simple.jsonnet +++ /dev/null @@ -1,321 +0,0 @@ -// This is the file cfg/layers/mids/pdsp/variant/simple.jsonnet. -// If this file is found at any other path, I will delete it without notice. - -// This files defines a simple VARIANT for pdsp which is NOT MEANT TO BE -// CORRECT. Do not "improve" it. It is meant to provide simple parameters to -// make debugging easier. - -// If you require a new variant, read and follow the comments at the top of -// "nominal.jsonnet". - -local wc = import "wirecell.jsonnet"; - -local detectors = import "detectors.jsonnet"; -local mydet = detectors.pdsp; -local nominal = import "nominal.jsonnet"; - -local simple = { - // Define the LAr properties. In principle, this is NOT subject - // to variance though in principle it could be (eg, top/bottom of - // DUNE-VD may have different LAr temps?). We include this here - // to keep data out of the the mid code and because this info is - // used in a few places. - lar : { - // Longitudinal diffusion constant - DL : 4.0 * wc.cm2/wc.s, - // Transverse diffusion constant - DT : 8.8 * wc.cm2/wc.s, - // Electron lifetime - lifetime : 10.4*wc.ms, - // Electron drift speed, assumes a certain applied E-field - // See https://github.com/WireCell/wire-cell-toolkit/issues/266 - // We do not pick a "simple" value (1.6) in order that we match - // what is *likely* in the FR file. - drift_speed : 1.565*wc.mm/wc.us, // at 500 V/cm - // LAr density - density: 1.389*wc.g/wc.centimeter3, - // Decay rate per mass for natural Ar39. - ar39activity: 1*wc.Bq/wc.kg, - }, - - // Extract all the location for planes used in geometry - planes : { - apa_cpa: 3.63075*wc.m, // LArSoft - // apa_cpa = 3.637*wc.m, // DocDB 203 - cpa_thick: 3.175*wc.mm, // 1/8", from Bo Yu (BNL) and confirmed with LArSoft - // cpa_thick = 50.8*wc.mm, // DocDB 203 - apa_w2w: 85.725*wc.mm, // between W planes, DocDB 203 - apa_g2g: 114.3*wc.mm, - // measured from grid to "anode" plane along drift direction - offset: 4.76*wc.mm - }, - - // This function returns the minimal geometry needed by WCT. In - // the transverse plane the geometry is defined by the "wires - // file" and in the drift direction by the "drifts" structure - // which defines a trio of planes: "anode", "response" and - // "cathode" which are perpendicular to drift. Each are assumed - // approximately equipotential along their plane. If your - // detector fails this, then use a "variant params" pack for each - // detector region which does meet this requirement. - // - geometry : { - - // The exhaustive list of the location of every single wire - // (or strip) segment. - // wires_file: "protodune-wires-larsoft-v4.json.bz2", - wires_file: mydet.wires, - - // Comments on how to chose the "anode" plane location: The - // "anode" cut off plane, here measured from APA centerline, - // determines how close to the wires do we consider any depo. - // Anything closer will simply be discarded, else it will - // either be drifted or "backed up" to the response plane. - // This is somewhat arbitrary choice. The "backing up" assume - // parallel drift points and so some small error is made when - // applied to ionization produce very close to the wires. To - // reduce this small error the "anode" plane can be moved - // closer to the "response" plane but at the cost of not - // simulate the volume left beyond the "anode" plane. - - // local apa_plane = 0.5*apa_g2g, // put anode plane at grid plane - local apa_plane = 0.5*$.planes.apa_g2g - $.planes.offset, // put anode plane at first induction plane - - // Response plane MUST be located so that the wire plane - // locations used in the FR calculation matches the wire plane - // locations given to WCT. - local res_plane = 0.5*$.planes.apa_w2w + $.ductor.response_plane, - - // The cathode plane is like the anode cut off plane. Any - // depo not between the two is dropped prior to drifting. It - // is demarks the face of the CPA. - local cpa_plane = $.planes.apa_cpa - 0.5*$.planes.cpa_thick, - - // The drifts are then defined in terms of these above - // numbers. You can use "wirecell-util wires-info" or - // "wirecell-util wires-volumes" or others to understand the - // mapping of anode number to the 6 locations in global X and - // Z coordinates. For Larsoft wires the numbering is column - // major starting at small X and Z so the centerline is - // -/+/-/+/-/+. Also important is that the faces are listed - // "front" first. Front is the one with the more positive X - // coordinates and if we want to ignore a face it is made - // null. - drifts : [{ - local sign = 2*(n%2)-1, - local centerline = sign*$.planes.apa_cpa, - wires: n, // anode number - name: "apa%d"%n, - faces: - // top, front face is against cryo wall - if sign > 0 - then [ - { - anode: centerline + apa_plane, - response: centerline + res_plane, - cathode: centerline + cpa_plane, - }, - { - anode: centerline - apa_plane, - response: centerline - res_plane, - cathode: centerline - cpa_plane, - } - ] - // bottom, back face is against cryo wall - else [ - { - anode: centerline + apa_plane, - response: centerline + res_plane, - cathode: centerline + cpa_plane, - }, - { - anode: centerline - apa_plane, - response: centerline - res_plane, - cathode: centerline - cpa_plane, - } - - ], - } for n in std.range(0,5)], - }, - - // The nominal time binning for data produced by the detector. - binning : { - tick: 0.5*wc.us, - nticks: 6000, - }, - - // The electronics response parameters. - elec : { - // The FE amplifier gain in units of Voltage/Charge. - gain : 14.0*wc.mV/wc.fC, - - // The shaping (aka peaking) time of the amplifier shaper. - shaping : 2.2*wc.us, - - // A realtive gain applied after shaping and before ADC. - // - // Don't use this to fix wrong sign depos. If you neeed to - // fix sign of eg larsoft depos its input component has a - // "scale" parameter which can be given a -1. Or better, fix - // the C++ that sets the wrong sign to begin with. - // - // Pulser calibration: 41.649 ADC*tick/1ke - // theoretical elec resp (14mV/fC): 36.6475 ADC*tick/1ke - postgain: 1.1365, - }, - - // The "RC" response - rc : { - width: 1.1*wc.ms, - }, - - // The ductor transforms drifted depos to currents - ductor: { - - # field_file: "dune-garfield-1d565.json.bz2", - field_file: mydet.fields, - - // The distance from the anode centerline to where the field - // response calculation begins drifting. Take care that field - // response calculations are not always done with the - // thickness of the anode in mind and may report their - // "response plane" measured relative to a number of other - // monuments. This value is almost certainly required in the - // "geometry" section. - response_plane : 10*wc.cm, - - // The time bin where the readout should be considered to - // actually start given the pre-signal simulated. - tbin : 0, - binning: { - // Ductor ticks at same speed as ADC - tick : $.binning.tick, - // Over the somewhat enlarged domain - nticks : $.ductor.tbin + $.binning.nticks, - }, - start_time : -self.response_plane / $.lar.drift_speed, - readout_time : self.binning.nticks * self.binning.tick, - }, - - // A "splat" (DepoFluxSplat) is an approximation to the combination of - // ductor+sigproc - splat : { - sparse: true, - tick: $.ductor.binning.tick, - window_start: $.ductor.start_time, - window_duration: $.ductor.readout_time, - reference_time: 0.0, - }, - - // Simulating noise - noise : { - model: { - #spectra_file: "protodune-noise-spectra-v1.json.bz2", - spectra_file: mydet.noise, - - // These are frequency space binning which are not necessarily - // same as some time binning - but here they are. - period: $.binning.tick, // 1/frequency - nsamples: $.binning.nticks, // - // Optimize binning - wire_length_scale: 1.0*wc.cm, - }, - // Optimize use of randoms - replacement_percentage: 0.02, - }, - - // Digitization model parameters. - digi : { - // A relative gain applied just prior to digitization. This - // is not FE gain, see elec for that. - gain: 1.0, - - // Voltage baselines added to any input voltage signal listed - // in a per plan (U,V,W) array. - // - // Per tdr, chapter 2 - // induction plane: 2350 ADC, collection plane: 900 ADC - baselines: [1003.4*wc.millivolt,1003.4*wc.millivolt,507.7*wc.millivolt], - - // The resolution (bits) of the ADC - resolution: 12, - - // The voltage range as [min,max] of the ADC, eg min voltage - // counts 0 ADC, max counts 2^resolution-1. - // - // The tdr says, "The ADC ASIC has an input - // buffer with offset compensation to match the output of the - // FE ASIC. The input buffer first samples the input signal - // (with a range of 0.2 V to 1.6 V)..." - fullscale: [0.2*wc.volt, 1.6*wc.volt], - }, - - // The noise filter parameter pack - nf : { - // In principle, may use a different field response. - field_file : $.ductor.field_file, - - binning: $.binning, - - // The name of the ChannelNoiseDb - chndb: "actual", - filters: { - channel: ["single"], // sticky, single, gaincalib - grouped: [], // grouped, - status: [], - }, - }, - - // Signal processing parameter pack - sp: { - // In principle, may use a different field response. - field_file : $.ductor.field_file, - binning : $.binning, - }, - - dnnroi: { - // this is a total fudge factor - output_scale : 1.0, - - // unet-l23-cosmic500-e50.ts was made for uboone.... - model_file : "unet-l23-cosmic500-e50.ts", - }, - - - // Imaging paramter pack - img : { - // For now we use MicroBooNE's - #"charge_error_file": "microboone-charge-error.json.bz2", - "charge_error_file": mydet.qerr, - - // Number of ticks to collect into one time slice span - span: 4, - - binning : $.binning, - - // The tiling strategy refers to how to deal with substantial - // dead channels. The "perfect" strategy simply performs so - // called "3-plane" tiling. The "mixed" strategy also - // includes so called "2-plane" tiling, producing more and - // more ambiguous blobs but will not have "gaps" due to dead - // channels. - tiling_strategy: "perfect", - }, - -}; - -local override = nominal + { - ductor: super.ductor + { - tbin : 0, - binning: { - tick : $.binning.tick, - nticks : $.ductor.tbin + $.binning.nticks, - }, - start_time : -self.response_plane / $.lar.drift_speed, - readout_time : self.binning.nticks * self.binning.tick, - } -}; - -// simple -override - diff --git a/cfg/layers/mids/pdsp/variants/ssss.jsonnet b/cfg/layers/mids/pdsp/variants/ssss.jsonnet index 8d6cef948..e3a136812 100644 --- a/cfg/layers/mids/pdsp/variants/ssss.jsonnet +++ b/cfg/layers/mids/pdsp/variants/ssss.jsonnet @@ -1,26 +1,13 @@ local wc = import "wirecell.jsonnet"; -local detectors = import "detectors.jsonnet"; -local mydet = detectors.pdsp; local nominal = import "nominal.jsonnet"; -local override = nominal + { - ductor: super.ductor + { - tbin : 0, - binning: { - tick : $.binning.tick, - nticks : $.ductor.tbin + $.binning.nticks, - }, - start_time : -self.response_plane / $.lar.drift_speed, - readout_time : self.binning.nticks * self.binning.tick, - } -}; - -local smeared = override + { +local smeared = nominal + { splat: super.splat + { // Run wirecell-gen morse-* to find these numbers that match the extra - // spread the sigproc induces. + // spread the sigproc induces. CAVEAT: these results are inderectly + // tied to the field response used. "smear_long": [ 2.691862363980221, 2.6750200122535057, @@ -37,12 +24,21 @@ local smeared = override + { { // Used for test-morse-pdsp - morse_nominal: override, + morse_nominal: nominal, // Used for test-ssss-pdsp - ssss_nominal: override, + ssss_nominal: nominal, ssss_smeared: smeared, - // Used for test-spdir - spdir: smeared + // Used for test/scripts/spdir. Note, these files are generated. + spdir_lofr: smeared { + detector_data: super.detector_data { + field: "pdsp-fields-lo.json.bz2" + }, + }, + spdir_hifr: smeared { + detector_data: super.detector_data { + field: "pdsp-fields-hi.json.bz2" + }, + }, } diff --git a/cfg/layers/mids/uboone/api.jsonnet b/cfg/layers/mids/uboone/api.jsonnet index a809ae16e..4b70115d9 100644 --- a/cfg/layers/mids/uboone/api.jsonnet +++ b/cfg/layers/mids/uboone/api.jsonnet @@ -1,25 +1,13 @@ // return a mid-level API object -local low = import "../../low.jsonnet"; +local nf = import "api/nf.jsonnet"; +local sp = import "api/sp.jsonnet"; local sim = import "api/sim.jsonnet"; -local sigproc = import "api/sigproc.jsonnet"; local img = import "api/img.jsonnet"; function(services, params, options) { + nf : nf(services, params), + sp : sp(services, params, options), + img : img(services, params), +} + sim(services, params) - anodes :: function() - low.anodes(params.geometry.drifts, params.geometry.wires_file), - - drifter :: function(name="") - low.drifter(services.random, - low.util.driftsToXregions(params.geometry.drifts), - params.lar, name=name), - - // track_depos, signal, noise, digitizer - sim :: sim(services, params), - - // nf, sp, dnnroi - sigproc :: sigproc(services, params, options), - - img :: img(services, params), -} diff --git a/cfg/layers/mids/uboone/api/frs.jsonnet b/cfg/layers/mids/uboone/api/frs.jsonnet index d519806bb..8901e4518 100644 --- a/cfg/layers/mids/uboone/api/frs.jsonnet +++ b/cfg/layers/mids/uboone/api/frs.jsonnet @@ -2,17 +2,22 @@ function(params) { sp: { type: "FieldResponse", name: "sp", - data: { filename: params.sp.field_file } + data: { filename: params.sp.field_file } // singular }, nf: { type: "FieldResponse", name: "nf", - data: { filename: params.nf.field_file } + data: { filename: params.nf.field_file } // singular }, - sims: [{ + sim : { + type: "FieldResponse", + name: params.ductor.field_file, + data: { filename: params.ductor.field_file } + }, + sims : [{ type: "FieldResponse", name: fn, data: { filename: fn } - } for fn in params.ductor.field_files], + } for fn in params.ductor.field_files], // plural } diff --git a/cfg/layers/mids/uboone/api/img.jsonnet b/cfg/layers/mids/uboone/api/img.jsonnet index 65d6f9b0c..531e4738e 100644 --- a/cfg/layers/mids/uboone/api/img.jsonnet +++ b/cfg/layers/mids/uboone/api/img.jsonnet @@ -1,193 +1,64 @@ -// Produce config for WC 3D imaging subgraph for microboone - -// This mostly transcribes: -// https://github.com/HaiwangYu/wcp-porting-img/blob/main/wct-celltree-img.jsonnet - -// There are lots of "magic" numbers here that likely should be in a -// variant parameter pack. - +// Produce config for WC 3D imaging subgraph local low = import "../../../low.jsonnet"; local pg = low.pg; local wc = low.wc; -local high = import "layers/high.jsonnet"; -// Return a subgraph that spans from frame to cluster. -// Shape of subgraph is controlled by params.img.slice_strategy -function(services, params) function(anode) +function(services, params) function(anode) local ident = low.util.idents(anode); local img = low.img(anode); - - //-------------- haiwangs - -local img = import "img.jsonnet"; - -local celltreesource = pg.pnode({ - type: "CelltreeSource", - name: "celltreesource", - data: { - filename: "celltreeOVERLAY.root", - EventNo: 6501, - // in_branch_base_names: raw [default], calibGaussian, calibWiener - in_branch_base_names: ["calibWiener", "calibGaussian"], - out_trace_tags: ["wiener", "gauss"], // orig, gauss, wiener - in_branch_thresholds: ["channelThreshold", ""], - }, - }, nin=0, nout=1); - -local dumpframes = pg.pnode({ - type: "DumpFrames", - name: 'dumpframe', - }, nin=1, nout=0); - -local magdecon = pg.pnode({ - type: 'MagnifySink', - name: 'magdecon', - data: { - output_filename: "mag.root", - root_file_mode: 'UPDATE', - frames: ['gauss', 'wiener', 'gauss_error'], - cmmtree: [['bad','bad']], - summaries: ['gauss', 'wiener', 'gauss_error'], - trace_has_tag: true, - anode: wc.tn(anodes[0]), - }, - }, nin=1, nout=1, uses=[anodes[0]]); - -local waveform_map = { - type: 'WaveformMap', - name: 'wfm', - data: { - filename: "microboone-charge-error.json.bz2", - }, - uses: [], + local wfm = { + type: 'WaveformMap', + name: ident, + data: { + filename: params.img.charge_error_file, + }, }; -local charge_err = pg.pnode({ - type: 'ChargeErrorFrameEstimator', - name: 'cefe', - data: { - intag: "gauss", - outtag: 'gauss_error', - anode: wc.tn(anodes[0]), - rebin: 4, // this number should be consistent with the waveform_map choice - fudge_factors: [2.31, 2.31, 1.1], // fudge factors for each plane [0,1,2] - time_limits: [12, 800], // the unit of this is in ticks - errors: wc.tn(waveform_map), - }, - }, nin=1, nout=1, uses=[waveform_map, anodes[0]]); - -local cmm_mod = pg.pnode({ - type: 'CMMModifier', - name: '', - data: { - cm_tag: "bad", - trace_tag: "gauss", - anode: wc.tn(anodes[0]), - start: 0, // start veto ... - end: 9592, // end of veto - ncount_cont_ch: 2, - cont_ch_llimit: [296, 2336+4800 ], // veto if continues bad channels - cont_ch_hlimit: [671, 2463+4800 ], - ncount_veto_ch: 1, - veto_ch_llimit: [3684], // direct veto these channels - veto_ch_hlimit: [3699], - dead_ch_ncount: 10, - dead_ch_charge: 1000, - ncount_dead_ch: 2, - dead_ch_llimit: [2160, 2080], // veto according to the charge size for dead channels - dead_ch_hlimit: [2176, 2096], - ncount_org: 5, // organize the dead channel ranges according to these boundaries - org_llimit: [0 , 1920, 3840, 5760, 7680], // must be ordered ... - org_hlimit: [1919, 3839, 5759, 7679, 9592], // must be ordered ... - }, - }, nin=1, nout=1, uses=[anodes[0]]); - -local frame_quality_tagging = pg.pnode({ - type: 'FrameQualityTagging', - name: '', - data: { - trace_tag: 'gauss', - anode: wc.tn(anodes[0]), - nrebin: 4, // rebin count ... - length_cut: 3, - time_cut: 3, - ch_threshold: 100, - n_cover_cut1: 12, - n_fire_cut1: 14, - n_cover_cut2: 6, - n_fire_cut2: 6, - fire_threshold: 0.22, - n_cover_cut3: [1200, 1200, 1800 ], - percent_threshold: [0.25, 0.25, 0.2 ], - threshold1: [300, 300, 360 ], - threshold2: [150, 150, 180 ], - min_time: 3180, - max_time: 7870, - flag_corr: 1, - }, - }, nin=1, nout=1, uses=[anodes[0]]); - -local frame_masking = pg.pnode({ - type: 'FrameMasking', - name: '', - data: { - cm_tag: "bad", - trace_tags: ['gauss','wiener'], - anode: wc.tn(anodes[0]), - }, - }, nin=1, nout=1, uses=[anodes[0]]); - -local anode = anodes[0]; -// multi slicing includes 2-view tiling and dead tiling -local active_planes = [[0,1,2],[0,1],[1,2],[0,2],]; -local masked_planes = [[],[2],[0],[1]]; -// single, multi, active, masked -function(slicing_strategy = "single") - - local imgpipe = { - single: pg.pipeline([ - img.slicing(anode, anode.name, "gauss", 109, active_planes=[0,1,2], masked_planes=[],dummy_planes=[]), // 109*22*4 - img.tiling(anode, anode.name), - img.solving(anode, anode.name), - // img.clustering(anode, anode.name), - img.dump_new(anode, anode.name, params.lar.drift_speed), - ], "img-" + anode.name), - active: pg.pipeline([ - img.multi_active_slicing_tiling(anode, anode.name+"-ms-active", "gauss", 4), - img.solving(anode, anode.name+"-ms-active"), - img.dump(anode, anode.name+"-ms-active", params.lar.drift_speed)]), - masked: pg.pipeline([ - // img.multi_masked_slicing_tiling(anode, anode.name+"-ms-masked", "gauss", 109), - img.multi_masked_2view_slicing_tiling(anode, anode.name+"-ms-masked", "gauss", 109), - img.clustering(anode, anode.name+"-ms-masked"), - img.dump(anode, anode.name+"-ms-masked", params.lar.drift_speed)]), - multi: { - local active_fork = pg.pipeline([ - img.multi_active_slicing_tiling(anode, anode.name+"-ms-active", "gauss", 4), - img.solving(anode, anode.name+"-ms-active"), - img.dump(anode, anode.name+"-ms-active", params.lar.drift_speed), - ]), - local masked_fork = pg.pipeline([ - // img.multi_masked_slicing_tiling(anode, anode.name+"-ms-masked", "gauss", 109), - img.multi_masked_2view_slicing_tiling(anode, anode.name+"-ms-masked", "gauss", 109), - img.clustering(anode, anode.name+"-ms-masked"), - img.dump(anode, anode.name+"-ms-masked", params.lar.drift_speed), - ]), - ret: pg.fan.fanout("FrameFanout",[active_fork,masked_fork], "fan_active_masked"), - }.ret, - }; - local graph = pg.pipeline([ - celltreesource, - // frame_quality_tagging, // event level tagging - // cmm_mod, // CMM modification - // frame_masking, // apply CMM - charge_err, // calculate charge error - // magdecon, // magnify out - // dumpframes, - imgpipe[slicing_strategy], - ], "main"); + local slicing(min_tbin=0, max_tbin=0, ext="", span=4, + active_planes=[0,1,2], masked_planes=[], dummy_planes=[]) = + pg.pnode({ + type: "MaskSlices", + name: ident+ext, + data: { + wiener_tag: "wiener"+ident, + charge_tag: "gauss"+ident, + error_tag: "gauss_error"+ident, + tick_span: span, + anode: wc.tn(anode), + min_tbin: min_tbin, + max_tbin: max_tbin, + active_planes: active_planes, + masked_planes: masked_planes, + dummy_planes: dummy_planes, + }, + }, nin=1, nout=1, uses=[anode]); + + + // Inject signal uncertainty. + local sigunc = + pg.pnode({ + type: 'ChargeErrorFrameEstimator', + name: ident, + data: { + intag: "gauss"+ident, + outtag: 'gauss_error'+ident, + anode: wc.tn(anode), + rebin: 4, // this number should be consistent with the waveform_map choice + fudge_factors: [2.31, 2.31, 1.1], // fudge factors for each plane [0,1,2] + time_limits: [12, 800], // the unit of this is in ticks + errors: wc.tn(wfm), + }, + }, nin=1, nout=1, uses=[wfm, anode]); + + + low.pg.pipeline([ + sigunc, + slicing(span=params.img.span), + img.tiling(), + img.clustering(), + img.grouping(), + img.charge_solving(), // fixme: a few optins we may want to allow to specify in variant params + ], ident) - high.main(graph, "TbbFlow") - diff --git a/cfg/layers/mids/uboone/api/nf.jsonnet b/cfg/layers/mids/uboone/api/nf.jsonnet index 08dfd064c..f8433864e 100644 --- a/cfg/layers/mids/uboone/api/nf.jsonnet +++ b/cfg/layers/mids/uboone/api/nf.jsonnet @@ -6,8 +6,7 @@ local wc = low.wc; local pg = low.pg; local chndb = import "details/chndb.jsonnet"; -function(services, params) function(anode) - local name = "%d"%anode.data.ident; +function(services, params) function(anode, name) local chndbobj = chndb(services, params, anode); @@ -78,8 +77,8 @@ function(services, params) function(anode) wc.tn(status), ], noisedb: wc.tn(chndbobj), - intraces: "orig", - outtraces: "quiet", + intraces: "", + outtraces: "", } }, uses=[chndbobj, anode, single, grouped, bitshift, status], nin=1, nout=1); @@ -87,8 +86,8 @@ function(services, params) function(anode) type: "OmnibusPMTNoiseFilter", name:name, data: { - intraces: "quiet", - outtraces: "raw", + intraces: "", + outtraces: "", anode: wc.tn(anode), } }, nin=1, nout=1, uses=[anode]); diff --git a/cfg/layers/mids/uboone/api/sigproc.jsonnet b/cfg/layers/mids/uboone/api/sigproc.jsonnet deleted file mode 100644 index b2ca5b784..000000000 --- a/cfg/layers/mids/uboone/api/sigproc.jsonnet +++ /dev/null @@ -1,21 +0,0 @@ -local low = import "../../../low.jsonnet"; - -// These are big and unique to uboone so we have to write them out -// exhaustively. -local nf = import "nf.jsonnet"; -local sp = import "sp.jsonnet"; - -function(services, params, options={}) { - - // API method - nf :: nf(services, params), - - // API method - sp :: sp(services, params, options), - - // API method - dnnroi :: function(anode) - low.dnnroi(anode, params.dnnroi.model_file, - services.platform, params.dnnroi.output_scale), - -} diff --git a/cfg/layers/mids/uboone/api/sim.jsonnet b/cfg/layers/mids/uboone/api/sim.jsonnet index 11db507e0..1f4826566 100644 --- a/cfg/layers/mids/uboone/api/sim.jsonnet +++ b/cfg/layers/mids/uboone/api/sim.jsonnet @@ -13,12 +13,11 @@ local idents = low.util.idents; function(services, params, options={}) { + // local res = low.resps(params).sim, + // Signal binning may be extended from nominal. local sig_binning = params.ductor.binning, - // Must be *3* field responses - local fr = frs(params).sims, - local broken_detector = [{ regions: sr.uv+sr.vy, mode:"reject", @@ -43,19 +42,20 @@ function(services, params, options={}) { data: sig_binning { width: params.rc.width } }, - // "sys status" is false for nominal. No sys_resp. local short_responses = [ cer, ], local long_responses = [ rc, rc ], local make_pirs = function(index) [{ + // A set of multiple FRs. 0-element is "nominal" + local uboone_frs = frs(params).sims, type: 'PlaneImpactResponse', name: "p"+std.toString(plane)+"i"+std.toString(index), - uses: [fr[index], services.dft] + short_responses + long_responses, + uses: [uboone_frs[index], services.dft] + short_responses + long_responses, data: sig_binning { plane: plane, dft: wc.tn(services.dft), - field_response: wc.tn(fr[index]), + field_response: wc.tn(uboone_frs[index]), short_responses: [wc.tn(sr) for sr in short_responses], long_responses: [wc.tn(lr) for lr in long_responses], overall_short_padding: 100*wc.us, @@ -63,7 +63,8 @@ function(services, params, options={}) { }, } for plane in [0,1,2]], - local transforms = function(anode, index, pirs) { + // Make a DepoTransform parameterized by its FR. + local transforms = function(anode, pirs, index=0) { local common = { rng: wc.tn(services.random), dft: wc.tn(services.dft), @@ -78,22 +79,57 @@ function(services, params, options={}) { nsigma: 3, }, local uses = pirs + [anode, services.random, services.dft], - standard: - pg.pnode({ - type: 'DepoTransform', - name: "a"+idents(anode)+"i"+std.toString(index), - data: common, - }, nin=1, nout=1, uses = uses), + standard: pg.pnode({ + type: 'DepoTransform', + name: "a"+idents(anode)+"i"+std.toString(index), + data: common, + }, nin=1, nout=1, uses = uses), // BIG FAT WARNING: this will not work in standard WCT and // requires std.extVar's to be set..... - overlay: - pg.pnode({ - type:'wclsReweightedDepoTransform', - name: "a"+idents(anode)+"i"+std.toString(index), - data: common + params.ductor.overlay, - }, nin=1, nout=1, uses = uses), + overlay: pg.pnode({ + type:'wclsReweightedDepoTransform', + name: "a"+idents(anode)+"i"+std.toString(index), + data: common + params.ductor.overlay, + }, nin=1, nout=1, uses = uses), }, + + local sys = { + type: 'ResponseSys', + data: sig_binning { + // overall_short_padding should take into account this offset + // "start". currently all "start" should be the same cause we only + // have an overall time offset compensated in reframer. These values + // correspond to files.fields[0, 1, 2] e.g. normal, shorted U, and + // shorted Y + start: [-10*wc.us, -10*wc.us, -10*wc.us], + magnitude: [1.0, 1.0, 1.0], + time_smear: [0.0*wc.us, 0.0*wc.us, 0.0*wc.us], + } + }, + + // Make some depos. Note, this does NOT have a bagger. + track_depos : function(tracklist = [{ + time: 0, + charge: -5000, + ray: { + tail: wc.point(-4.0, 0.0, 0.0, wc.m), + head: wc.point(+4.0, 6.1, 7.0, wc.m), + }}]) + low.gen.track_depos(tracklist), + + // Simple signal simulation with a single FR. + signal_single(anode, name) : pg.pipeline([ + local pirs = make_pirs(index=0); + pg.pipeline([ + low.drifter(services.random, + low.util.driftsToXregions(params.geometry.volumes), + params.lar, name=name), + transforms(anode, pirs)[params.ductor.transform], + low.reframer(params, anode, name=name), + ])]), + + // For signal_multiple(). local make_branch = function(anode, index) local pirs = make_pirs(index); pg.pipeline([ @@ -113,52 +149,36 @@ function(services, params, options={}) { mode: broken_detector[index].mode, }, }, nin=1, nout=1, uses = [anode])]), - transforms(anode, index, pirs)[params.ductor.transform], + transforms(anode, pirs, index)[params.ductor.transform], ]), - - local sys = { - type: 'ResponseSys', - data: sig_binning { - // overall_short_padding should take into account this offset "start". - // currently all "start" should be the same cause we only have an overall time offset - // compensated in reframer - // These values correspond to files.fields[0, 1, 2] - // e.g. normal, shorted U, and shorted Y - start: [-10*wc.us, -10*wc.us, -10*wc.us], - magnitude: [1.0, 1.0, 1.0], - time_smear: [0.0*wc.us, 0.0*wc.us, 0.0*wc.us], - } - }, - - // Make some depos. Note, this does NOT have a bagger. - track_depos :: function(tracklist = [{ - time: 0, - charge: -5000, - ray: { - tail: wc.point(-4.0, 0.0, 0.0, wc.m), - head: wc.point(+4.0, 6.1, 7.0, wc.m), - }}]) - low.gen.track_depos(tracklist), - - local ubsigtags = ['ubsig%d'%n for n in [0,1,2]], - - // API method sim.signal: subgraph making pure signal voltage from - // depos. - signal :: function(anode) pg.pipeline([ - pg.fan.pipe('DepoSetFanout', [make_branch(anode, index) for index in [0,1,2]], + // More complex signal subgraph needed for "actual". + signal_multiple(anode, name) : + local ubsigtags = ['ubsig%d'%n for n in [0,1,2]]; + pg.pipeline([ + // Super position of multiple FRs + pg.fan.pipe('DepoSetFanout', + [make_branch(anode, index) for index in [0,1,2]], 'FrameFanin', 'ubsigraph', ubsigtags), - pg.pnode({ + // A special reframer with a magic number time offset. + pg.pnode({ type: 'Reframer', - name: idents(anode), - data: params.reframer { + name: name, + data: { anode: wc.tn(anode), + local wtf_is_this = 81*wc.us, + tbin: wtf_is_this/(params.ductor.binning.tick), + nticks: params.ductor.binning.nticks, + toffset: params.ductor.drift_dt - wtf_is_this, tags: ubsigtags, fill: 0.0, }, }, nin=1, nout=1, uses=[anode])]), - + + signal : self["signal_" + params.ductor.universe], + + local csdb = { type: "StaticChannelStatus", name: "", @@ -173,10 +193,11 @@ function(services, params, options={}) { } }, - noise :: function(anode) + + noise(anode, name) : local model = { type: 'EmpiricalNoiseModel', - name: idents(anode), + name: name, data: params.noise.model { anode: wc.tn(anode), dft: wc.tn(services.dft), @@ -186,7 +207,7 @@ function(services, params, options={}) { }; pg.pnode({ type: 'AddNoise', - name: idents(anode), + name: name, data: { rng: wc.tn(services.random), dft: wc.tn(services.dft), @@ -195,15 +216,6 @@ function(services, params, options={}) { replacement_percentage: params.noise.replacement_percentage, }}, nin=1, nout=1, uses=[services.random, services.dft, model]), - digitizer :: function(anode) - pg.pnode({ - type: 'Digitizer', - name: idents(anode), - data: params.digi { - anode: wc.tn(anode), - frame_tag: "orig" + idents(anode), - } - }, nin=1, nout=1, uses=[anode]), } diff --git a/cfg/layers/mids/uboone/api/sp.jsonnet b/cfg/layers/mids/uboone/api/sp.jsonnet index 25506d4ce..307db27fb 100644 --- a/cfg/layers/mids/uboone/api/sp.jsonnet +++ b/cfg/layers/mids/uboone/api/sp.jsonnet @@ -12,16 +12,18 @@ local pg = low.pg; local frs = import "frs.jsonnet"; -function(services, params, options={}) function(anode) +function(services, params, options={}) function(anode, name) local fr = frs(params).nf; local cer = params.ductor.binning { type: "ColdElecResponse", + name: name, data: params.elec, }; local perchanresp = { type: "PerChannelResponse", + name: name, data: { filename: params.nf.chresp_file, } @@ -29,6 +31,7 @@ function(services, params, options={}) function(anode) local sigproc_perchan = pg.pnode({ type: "OmnibusSigProc", + name: name, data: { // This class has a HUGE set of parameters. See // OmnibusSigProc.h for the list. For here, for now, we @@ -49,6 +52,7 @@ function(services, params, options={}) function(anode) local sigproc_uniform = pg.pnode({ type: "OmnibusSigProc", + name: name, data: { anode: wc.tn(anode), dft: wc.tn(services.dft), @@ -76,16 +80,17 @@ function(services, params, options={}) function(anode) local rawsplit = pg.pnode({ type: "FrameSplitter", - name: "rawsplitter" + name: "rawsplitter"+name, }, nin=1, nout=2); local sigsplit = pg.pnode({ type: "FrameSplitter", - name: "sigsplitter" + name: "sigsplitter"+name }, nin=1, nout=2); local chsel = pg.pnode({ type: "ChannelSelector", + name: name, data: { // channels that will get L1SP applied channels: std.range(3566,4305), @@ -135,14 +140,14 @@ function(services, params, options={}) function(anode) // from normal SP for input to L1SP local rawsigmerge = pg.pnode({ type: "FrameMerger", - name: "rawsigmerge", + name: "rawsigmerge"+name, data: { rule: "replace", // note: the first two need to match the order of what data is // fed to ports 0 and 1 of this component in the pgraph below! mergemap: [ - ["raw","raw","raw"], + ["","","raw"], ["gauss","gauss","gauss"], ], } @@ -153,14 +158,14 @@ function(services, params, options={}) function(anode) // to be found. local l1merge = pg.pnode({ type: "FrameMerger", - name: "l1merge", + name: "l1merge"+name, data: { rule: "replace", // note: the first two need to match the order of what data is // fed to ports 0 and 1 of this component in the pgraph below! mergemap: [ - ["raw","raw","raw"], + // ["raw","raw","raw"], ["l1sp","gauss","gauss"], ["l1sp","wiener","wiener"], ], @@ -179,4 +184,4 @@ function(services, params, options={}) function(anode) pg.edge(chsel, l1spfilter), pg.edge(l1spfilter, l1merge), ], - name="L1SP") + name="L1SP"+name) diff --git a/cfg/layers/mids/uboone/mids.jsonnet b/cfg/layers/mids/uboone/mids.jsonnet index 1fd60396b..fe581bf72 100644 --- a/cfg/layers/mids/uboone/mids.jsonnet +++ b/cfg/layers/mids/uboone/mids.jsonnet @@ -4,9 +4,8 @@ local api = import "api.jsonnet"; // Import all variant parameter packs to allow for dict like lookup. -local variants = { - nominal: import "variants/nominal.jsonnet", -}; +local variants = import "variants.jsonnet"; function(services, variant="", options={}) api(services, variants[variant], options) + diff --git a/cfg/layers/mids/uboone/variants.jsonnet b/cfg/layers/mids/uboone/variants.jsonnet new file mode 100644 index 000000000..a0ae02f86 --- /dev/null +++ b/cfg/layers/mids/uboone/variants.jsonnet @@ -0,0 +1,12 @@ +// Import all variant parameter packs to allow for dict like lookup. + +local low = import "../../low.jsonnet"; +//local ssss = import "variants/ssss.jsonnet"; + +local nominal = import "variants/nominal.jsonnet"; +local actual = import "variants/actual.jsonnet"; + +{ + nominal: nominal, + actual: actual, +} + low.ssss(nominal) diff --git a/cfg/layers/mids/uboone/variants/actual.jsonnet b/cfg/layers/mids/uboone/variants/actual.jsonnet new file mode 100644 index 000000000..8e4869370 --- /dev/null +++ b/cfg/layers/mids/uboone/variants/actual.jsonnet @@ -0,0 +1,142 @@ +// This sort of follows +// uboonedata/WireCellData/pgrapher/experiment/uboone/*.jsonnet +// v08_00_00. However, it assumes a "pure" WCT implementation. See +// the "overlay" variant for what is needed to configure the official +// MicroBooNE "overlay simulation". + +local wc = import "wirecell.jsonnet"; +local nominal = import "nominal.jsonnet"; + +nominal { + + binning : { + tick: 0.5*wc.us, + nticks: 9595, + }, + + // These parameters only make sense for running WCT simulation on + // microboone in larsoft. The "trigger" section is not + // "standard". This section is just a set up for use below in + // "sim". There is no trigger, per se, in the simulation but + // rather a contract between the generators of energy depositions + // (ie, LarG4) and the drift and induction simulation (WCT). For + // details of this contract see: + // https://microboone-docdb.fnal.gov/cgi-bin/private/ShowDocument?docid=12290 + trigger : { + + // A hardware trigger occurs at some "absolute" time but near + // 0.0 for every "event". It is measured in "simulation" time + // which is the same clock used for expressing energy + // deposition times. The following values are from table 3 of + // DocDB 12290. + hardware_times: { + none: 0.0, + + // BNB hardware trigger time. Note interactions + // associated with BNB neutrinos should all be produced + // starting in the beam gate which begins at 3125ns and is + // 1600ns in duration. + bnb : -31.25*wc.ns, + + // Like above but for NUMI. It's gate is 9600ns long starting + // at 4687.5ns. + numi : -54.8675*wc.ns, + ext : -414.0625*wc.ns, + mucs: -405.25*wc.ns, + }, + hw_time: self.hardware_times["bnb"], + + // Measured relative to the hardware trigger time above is a + // time offset to the time that the first tick of the readout + // should sample. This is apparently fixed for all hardware + // trigger types (?). + time_offset: -1.6*wc.ms, + + time: self.hw_time + self.time_offset, + }, + + geometry : super.geometry + { + xplanes: super.xplanes { + // Comments from old config: "plane, arbitrary choice. Microboone + // wires put collection plane at absolute x=-6mm, response.plane_dx + // is measured relative to collection plane wires.". I (bv) do not + // know why this means subtracting 6mm.... + dresponse: 10*wc.cm - 6*wc.mm, + } + }, + + ductor: super.ductor { + + // Note, in another variant, perhaps one named 'larsoft', the transform + // value may be set to "overlay". But, not here. + transform: "standard", + + // Tickle API to set up a super-position of multiple FRs. + universe: "multiple", + + response_plane : $.geometry.xplanes.dresponse, + + local tzero = 0*wc.us, + drift_dt : self.response_plane / $.lar.drift_speed, + tbin : wc.roundToInt(self.drift_dt / self.binning.tick), + local nticks = super.binning.nticks, + binning: super.binning { + nticks : $.ductor.tbin + nticks, + }, + start_time : tzero - self.drift_dt + $.trigger.time, + readout_time : self.binning.nticks * self.binning.tick, + }, + + // Simulating noise + noise : { + model: { + spectra_file: "microboone-noise-spectra-v2.json.bz2", + // These are frequency space binning which are not necessarily + // same as some time binning - but here they are. + period: $.binning.tick, // 1/frequency The N number of samples + // in MicroBooNE noise filtering spectra is defined independently + // from how many ticks are actually read out in ata. + nsamples: 9592, + // Optimize binning + wire_length_scale: 1.0*wc.cm, + }, + // Optimize use of randoms + replacement_percentage: 0.02, + }, + + digi : { + // There are post-FE amplifiers. Should this be + // elec.postgain? (could be, but it's fine here. There are 2 + // relative gains available). + gain: 1.0, + + // The resolution (bits) of the ADC + resolution: 12, + + // These values are pre-amp (ASIC) pedestals. Two additional AC coupling + // would change the baseline and it also depends on the bias voltage from regulator. + // +/-10 ADC? variation is expected. + //baselines: [900*wc.millivolt,900*wc.millivolt,200*wc.millivolt], + // Actuall MicroBooNE baselines are 2046 ADC for induction planes and 473 for collection. + baselines: [999.3*wc.millivolt,999.3*wc.millivolt,231.02*wc.millivolt], + + fullscale: [0*wc.volt, 2.0*wc.volt], + }, + + // Describe misconfigured channels. Used by NF's chndb and sim's + // static channel status. Nominally, everything is perfect. See + // other variants for non-perfect content. + misconfigured: { + // The list of misconfigured channel IDs + channels: [], + // Misconfigured parameters, only relevant if there are + // channels. Nominal case is there is no misconfig so set + // same as nominal config. + gain: $.elec.gain, + shaping: $.elec.shaping, + }, + + +} + + diff --git a/cfg/layers/mids/uboone/variants/data.jsonnet b/cfg/layers/mids/uboone/variants/data.jsonnet new file mode 100644 index 000000000..e69de29bb diff --git a/cfg/layers/mids/uboone/variants/geometry.jsonnet b/cfg/layers/mids/uboone/variants/geometry.jsonnet new file mode 100644 index 000000000..2cd7d96f0 --- /dev/null +++ b/cfg/layers/mids/uboone/variants/geometry.jsonnet @@ -0,0 +1,33 @@ +local wc = import 'wirecell.jsonnet'; +{ + + wires_file: 'microboone-celltree-wires-v2.1.json.bz2', + + // distance between collection wire plane and a plane. + xplanes: { + danode: 10.0 * wc.mm, + dresponse: 100.0 * wc.mm, + dcathode: 1000.0 * wc.mm, + }, + local xplanes = self.xplanes, // to make available below + + volumes: [ + + { + wires: 0, + xcenter: -6.0 * wc.mm, // absolute center of APA + local xcenter = self.xcenter, // to make available below. + faces: [ + + { + local dcollection = 6.0 * wc.mm, + anode: xcenter + (xplanes.danode + dcollection), + response: xcenter + (xplanes.dresponse + dcollection), + cathode: xcenter + (xplanes.dcathode + dcollection), + }, + null, + ], + }, + + ], +} diff --git a/cfg/layers/mids/uboone/variants/nominal.jsonnet b/cfg/layers/mids/uboone/variants/nominal.jsonnet index e35662b6a..870286039 100644 --- a/cfg/layers/mids/uboone/variants/nominal.jsonnet +++ b/cfg/layers/mids/uboone/variants/nominal.jsonnet @@ -1,12 +1,10 @@ -// This sort of follows -// uboonedata/WireCellData/pgrapher/experiment/uboone/*.jsonnet -// v08_00_00. However, it assumes a "pure" WCT implementation. See -// the "overlay" variant for what is needed to configure the official -// MicroBooNE "overlay simulation". +// Nominal is in name only. See actual for more realistic. local wc = import "wirecell.jsonnet"; +local base_variants = import "../../base/variants.jsonnet"; + +base_variants.nominal { -{ lar : { // Longitudinal diffusion constant DL : 6.4 * wc.cm2/wc.s, @@ -19,38 +17,25 @@ local wc = import "wirecell.jsonnet"; drift_speed : 1.098*wc.mm/wc.us, }, - geometry: { - wires_file: "microboone-celltree-wires-v2.1.json.bz2", - drifts: [{ - wires: 0, - name: "uboone", - faces: [ - { - // drop any depos w/in this plane. The exact - // choice represents some trade off in - // approximations. - anode: 0.0*wc.mm, - // plane, arbitrary choice. Microboone wires - // put collection plane at absolute x=-6mm, - // response.plane_dx is measured relative to - // collection plane wires. - response: $.ductor.response_plane - 6*wc.mm, - // Location of cathode measured from - // collection is based on a dump of - // ubcore/v06_83_00/gdml/microboonev11.gdml by - // Matt Toups - cathode: 2.5480*wc.m, - }, - null - ], - }] + // $ wirecell-util wires-volumes uboone \ + // | jsonnetfmt -n 4 --max-blank-lines 1 - \ + // > cfg/layers/mids/uboone/geometry.jsonnet + geometry_data : import "geometry.jsonnet", + geometry : $.geometry_data + { + xplanes: { + // Distances between collection plane and a logical plane. + danode: 0.0, + dresponse: 10*wc.cm, + // Location of cathode measured from collection is based on a dump + // of ubcore/v06_83_00/gdml/microboonev11.gdml by Matt Toups + dcathode: 2.5480*wc.m, + } }, // The nominal time binning for data produced by the detector. binning : { tick: 0.5*wc.us, - // Real detector makes 9595. - nticks: 9595, + nticks: 8192, }, elec : { @@ -63,78 +48,24 @@ local wc = import "wirecell.jsonnet"; width: 1.0*wc.ms, }, - // These parameters only make sense for running WCT simulation on - // microboone in larsoft. The "trigger" section is not - // "standard". This section is just a set up for use below in - // "sim". There is no trigger, per se, in the simulation but - // rather a contract between the generators of energy depositions - // (ie, LarG4) and the drift and induction simulation (WCT). For - // details of this contract see: - // https://microboone-docdb.fnal.gov/cgi-bin/private/ShowDocument?docid=12290 - trigger : { - - // A hardware trigger occurs at some "absolute" time but near - // 0.0 for every "event". It is measured in "simulation" time - // which is the same clock used for expressing energy - // deposition times. The following values are from table 3 of - // DocDB 12290. - hardware_times: { - none: 0.0, - - // BNB hardware trigger time. Note interactions - // associated with BNB neutrinos should all be produced - // starting in the beam gate which begins at 3125ns and is - // 1600ns in duration. - bnb : -31.25*wc.ns, - - // Like above but for NUMI. It's gate is 9600ns long starting - // at 4687.5ns. - numi : -54.8675*wc.ns, - ext : -414.0625*wc.ns, - mucs: -405.25*wc.ns, - }, - hw_time: self.hardware_times["bnb"], - - // Measured relative to the hardware trigger time above is a - // time offset to the time that the first tick of the readout - // should sample. This is apparently fixed for all hardware - // trigger types (?). - time_offset: -1.6*wc.ms, - - time: self.hw_time + self.time_offset, - }, - - ductor: { - + ductor: super.ductor { + // These tickle the API transform: "standard", // see overlay + universe: "single", // no super-position of different FRs - // We trifurcate the detector to deal with different - // problematic wire regions. - field_files: ["ub-10-half.json.bz2", - "ub-10-uv-ground-tuned-half.json.bz2", - "ub-10-vy-ground-tuned-half.json.bz2"], - - response_plane : 10*wc.cm, - - - local tzero = 0*wc.us, - drift_dt : self.response_plane / $.lar.drift_speed, - tbin : wc.roundToInt(self.drift_dt / self.binning.tick), - binning: { - tick : $.binning.tick, + // MB uses one or three field files. The standard variable is + // "field_file" (singular) but to simplify structure we define the full + // set even for nominal. For nominal ("single" universe) the first in + // the list is used. + field_files: $.detector_data.fields, - nticks : $.ductor.tbin + $.binning.nticks, - }, - start_time : tzero - self.drift_dt + $.trigger.time, - readout_time : self.binning.nticks * self.binning.tick, }, - // nominal: "sys_status" is FALSE and there is no "sys_resp". - - reframer : { - tbin: (81*wc.us)/($.binning.tick), - nticks: $.binning.nticks, - toffset: $.ductor.drift_dt - 81*wc.us, + splat : { + sparse: true, + tick: $.ductor.binning.tick, + window_start: $.ductor.start_time, + window_duration: $.ductor.readout_time, }, // Simulating noise @@ -147,7 +78,7 @@ local wc = import "wirecell.jsonnet"; // The MicroBooNE noise filtering spectra runs with a // different number of bins in frequency space than the // nominal number of ticks in time. - nsamples: 9592, + nsamples: 8192, // Optimize binning wire_length_scale: 1.0*wc.cm, }, @@ -175,9 +106,9 @@ local wc = import "wirecell.jsonnet"; }, nf: { - field_file: "ub-10-half.json.bz2", + field_file: $.detector_data.field, // "ub-10-half.json.bz2", binning: $.binning, - chresp_file: "microboone-channel-responses-v1.json.bz2", + chresp_file: $.detector_data.chresp, // The MicroBooNE noise filtering spectra runs with a // different number of bins in frequency space than the @@ -191,6 +122,11 @@ local wc = import "wirecell.jsonnet"; chndb_epoch: "perfect", }, + sp: { + field_file: $.detector_data.field, // "ub-10-half.json.bz2", + }, + + // FIXME: excise this from nominal. // Describe misconfigured channels. Used by NF's chndb and sim's // static channel status. Nominally, everything is perfect. See // other variants for non-perfect content. @@ -204,20 +140,11 @@ local wc = import "wirecell.jsonnet"; shaping: $.elec.shaping, }, - sp: { - field_file: "ub-10-half.json.bz2", - }, - // Imaging paramter pack img : { - charge_error_file: "microboone-charge-error.json.bz2", - - // Number of ticks to collect into one time slice span - span: 4, - + charge_error_file: $.detector_data.qerr, + span: 4, // Number of ticks to collect into one time slice span binning : $.binning, - - // fixme: remove old tiling strategy "perfect" and add slicing strategies }, } diff --git a/cfg/layers/mids/uboone/variants/resample.jsonnet b/cfg/layers/mids/uboone/variants/resample.jsonnet new file mode 100644 index 000000000..006ee9725 --- /dev/null +++ b/cfg/layers/mids/uboone/variants/resample.jsonnet @@ -0,0 +1,26 @@ +local wc = import "wirecell.jsonnet"; +local nominal = import "nominal.jsonnet"; + +{ + resample_500ns : nominal + { + binning : { + tick: 500 * wc.ns, + nticks: 6000, + }, + + ductor: super.ductor { + field_file: "uboone-100ns.json.bz2", + }, + }, + resample_512ns : nominal + { + binning : { + tick: 512 * wc.ns, + nticks: 6000 * 512 / 500, + }, + + ductor: super.ductor { + field_file: "uboone-64ns.json.bz2", + }, + + }, +} diff --git a/cfg/layers/mids/uboone/variants/ssss.jsonnet b/cfg/layers/mids/uboone/variants/ssss.jsonnet new file mode 100644 index 000000000..8bf045722 --- /dev/null +++ b/cfg/layers/mids/uboone/variants/ssss.jsonnet @@ -0,0 +1,45 @@ +local wc = import "wirecell.jsonnet"; + +local nominal = import "nominal.jsonnet"; + +local smeared = nominal + { + splat: super.splat + { + + // Run wirecell-gen morse-* to find these numbers that match the extra + // spread the sigproc induces. CAVEAT: these results are inderectly + // tied to the field response used. + // CAVEAT2 and FIXME: these from PDSP + "smear_long": [ + 2.691862363980221, + 2.6750200122535057, + 2.7137567141154055 + ], + "smear_tran": [ + 0.7377218875719689, + 0.7157764520393882, + 0.13980698710556544 + ] + } +}; + + +{ + // Used for test-morse-pdsp + morse_nominal: nominal, + + // Used for test-ssss-pdsp + ssss_nominal: nominal, + ssss_smeared: smeared, + + // Used for test/scripts/spdir + spdir_lofr: smeared { + detector: super.detector { + fields: "uboone-fields-lo.json.bz2" + }, + }, + spdir_hifr: smeared { + detector: super.detector { + fields: "uboone-fields-hi.json.bz2" + }, + }, +} diff --git a/cfg/layers/omnijob.jsonnet b/cfg/layers/omnijob.jsonnet new file mode 100644 index 000000000..fff7ba30d --- /dev/null +++ b/cfg/layers/omnijob.jsonnet @@ -0,0 +1,137 @@ +// This is a top level config file given to wire-cell to run segments of +// processing to calcualte the "standard" SP metrics (eff, bias, res) vs +// direction of ideal tracks. +// +// The caller defines a pipeline of tasks that default to sim,nf,sp. A subset +// of this task list may be given. Eg, UVCGAN may be inserted to perform domain +// translations and that may reasonably be done before or after the nf stage. +// +// Example usage: +// +// wire-cell -c spdir-metric.jsonnet \ +// -A detector=pdsp \ +// -A variant=nominal \ +// -A input=depos.npz \ +// -A output=sp.npz \ +// -A tasks=sim,nf,sp + +local high = import "layers/high.jsonnet"; +local wc = high.wc; +local pg = high.pg; + +// Main pipeline elements. +local builder(mid, anode, stages, outputs, dense=true) = { + local last_stage = stages[std.length(stages)-1], + main : { + drift: [ + // explicitly filter depos in anode. + pg.pnode({ + type:"DepoSetFilter", + name: "predrift", + data:{ + anode: wc.tn(anode), + }, + }, nin=1, nout=1, uses=[anode]), + mid.drifter(), + ], + splat: [ + // nothing inline, the pipline with the splat itself is part of the + // sink branch of the tap + ], + sim: [ + mid.signal(anode), + mid.noise(anode), + mid.digitizer(anode), + ], + sig: [ + mid.signal(anode), + ], + noi: [ + mid.noise(anode), + ], + dig: [ + mid.digitizer(anode), + ], + nf: [ + mid.nf(anode), + ], + sp: [ + mid.sp(anode), + ] + }, + pre_sink(stage) :: + if stage == "splat" + then [ mid.splat(anode) ] + else [], + + reframer(stage) :: + local reframers = { + splat: [mid.reframer(anode, name=outputs[stage])], + sp: [mid.reframer(anode, name=outputs[stage], tags=["gauss"])], + }; + if dense + then std.get(reframers, stage, []) + else [], + + file_sink(stage) :: [ + if stage == "drift" + then high.fio.depo_file_sink(outputs.drift) + else high.fio.frame_file_sink(outputs[stage]) + ], + + sink(stage) :: + pg.pipeline(self.pre_sink(stage) + self.reframer(stage) + self.file_sink(stage)), + + tap_or_sink(stage): + if stage == last_stage then + self.sink(stage) + else + high.fio.tap(if stage == "drift" || stage == "splat" + then "DepoSetFanout" + else "FrameFanout", + self.sink(stage), name=outputs[stage]), + + get_stage(stage): + if std.objectHas(outputs, stage) then + self.main[stage] + [self.tap_or_sink(stage)] + else + self.main[stage], +}; + + +local source(stage, input) = + if std.member(["drift","splat","sim"], stage) + then high.fio.depo_file_source(input) + else high.fio.frame_file_source(input); + + +local output_objectify(stages, output) = + local last_stage = stages[std.length(stages)-1]; + if std.type(output) == "string" then + {[last_stage]: output} + else output; + +// Return configuration for single-APA job. +// - detector :: canonical name of supported detector (pdsp, uboone, etc) +// - input :: name of file provding data to input to first task +// - ouput :: name of file to receive output of last task or object mapping task to output file +// - tasks :: array or comma separated list of task names +// - dense :: if false, save frames sparsely, else add reframers to make dense +// - variant :: the layers mids detector variant name +function (detector, input, output, tasks="drift,splat,sim,nf,sp", dense=true, variant="nominal") + local params = high.params(detector, variant); + local mid = high.api(detector, params); + local stages = wc.listify(tasks); + local outfiles = output_objectify(stages, output); // stage->filename + + local anodes = mid.anodes(); + local anode = anodes[0]; + + local b = builder(mid, anode, stages, outfiles, dense); + + local head = [source(stages[0], input)]; + local guts = [b.get_stage(stage) for stage in stages]; + local body = std.flattenArrays(guts); + + local graph = pg.pipeline(head + body); + high.main(graph, "Pgrapher") diff --git a/cfg/layers/test/test-api.jsonnet b/cfg/layers/test/test-api.jsonnet new file mode 100644 index 000000000..090014c78 --- /dev/null +++ b/cfg/layers/test/test-api.jsonnet @@ -0,0 +1,13 @@ +local high = import "layers/high.jsonnet"; +local detectors = ["pdsp", "uboone"]; +local variants = ["nominal", "spdir_hifr"]; + +{ + [det] : { + params :: high.params(det, "nominal"), + api :: high.api(det, self.params), + methods : std.objectFieldsAll(self.api), + nmethods: std.length(self.methods), + nvariants: std.length([high.params(det, var) for var in variants]), + } for det in detectors +} diff --git a/cfg/layers/test/test-base.jsonnet b/cfg/layers/test/test-base.jsonnet new file mode 100644 index 000000000..931e3a9b9 --- /dev/null +++ b/cfg/layers/test/test-base.jsonnet @@ -0,0 +1,9 @@ +local basebase = import "../mids/base/variants.jsonnet"; +local base = basebase.nominal { + detector_name: "base", + variant_name: "nominal", + structure_name: "nominal" +}; + +std.assertEqual("base", base.detector_data.detname) + diff --git a/cfg/layers/test/test-drifter.jsonnet b/cfg/layers/test/test-drifter.jsonnet new file mode 100644 index 000000000..9ed7d7ff5 --- /dev/null +++ b/cfg/layers/test/test-drifter.jsonnet @@ -0,0 +1,19 @@ + +local high = import "layers/high.jsonnet"; +local low = import "layers/low.jsonnet"; +local dr = import "layers/low/drifter.jsonnet"; + +{ + // params :: high.params("pdsp", "nominal"), + + // api :: high.api("pdsp", self.params), + + ran : high.services.random, + // xreg : low.util.driftsToXregions(self.params.geometry.drifts), + // lar: self.params.lar, + + // drifter: self.api.drifter() + // drifter : dr(high.services.random, + // low.util.driftsToXregions(self.params.geometry.drifts), + // self.params.lar), +} diff --git a/cfg/layers/test/test-params.jsonnet b/cfg/layers/test/test-params.jsonnet new file mode 100644 index 000000000..9bb50443f --- /dev/null +++ b/cfg/layers/test/test-params.jsonnet @@ -0,0 +1,16 @@ +local high = import "../high.jsonnet"; + +local un = high.params("uboone", "nominal"); +local ua = high.params("uboone", "actual"); + +[ + std.assertEqual(un.detector_name, "uboone"), + std.assertEqual(un.variant_name, "nominal"), + std.assertEqual(un.detector_data.wires, "microboone-celltree-wires-v2.1.json.bz2"), + + std.assertEqual(ua.detector_name, "uboone"), + std.assertEqual(ua.variant_name, "actual"), + std.assertEqual(ua.detector_data.wires, "microboone-celltree-wires-v2.1.json.bz2"), + +] + diff --git a/cfg/test/test-layer-actual.bats b/cfg/test/test-layer-actual.bats new file mode 100644 index 000000000..54defe1ca --- /dev/null +++ b/cfg/test/test-layer-actual.bats @@ -0,0 +1,19 @@ +#!/usr/bin/env bats + +# This tests the "cfg/layers/" config against the original "cfg/pgrapher/" + +bats_load_library wct-bats.sh + +@test "test layer actual wcls sim nf sp" { + cd_tmp file + + wcsonnet -o pgrapher-wct-sim-ideal-sn-nf-sp.json \ + pgrapher/experiment/uboone/wct-sim-ideal-sn-nf-sp.jsonnet + wcsonnet -o layers-wcls-sim-nf-sp.json layers/omnijob.jsonnet \ + -A tasks=sim,nf,sp \ + -A input=input.npz \ + -A output=output.npz \ + -A detector=uboone \ + -A variant=actual + +} diff --git a/cfg/test/test-mergepatch.jsonnet b/cfg/test/test-mergepatch.jsonnet new file mode 100644 index 000000000..7e1cd6dfe --- /dev/null +++ b/cfg/test/test-mergepatch.jsonnet @@ -0,0 +1,51 @@ +local wc = import "wirecell.jsonnet"; + +local g1 = { + wires_file: 'microboone-celltree-wires-v2.1.json.bz2', + + // distance between collection wire plane and a plane. + xplanes: { + danode: 10.0 * wc.mm, + dresponse: 100.0 * wc.mm, + dcathode: 1000.0 * wc.mm, + }, + local xplanes = self.xplanes, // to make available below + + volumes: [ + + { + wires: 0, + xcenter: -6.0 * wc.mm, // absolute center of APA + local xcenter = self.xcenter, // to make available below. + faces: [ + + { + local dcollection = 6.0 * wc.mm, + anode: xcenter + (xplanes.danode + dcollection), + response: xcenter + (xplanes.dresponse + dcollection), + cathode: xcenter + (xplanes.dcathode + dcollection), + }, + null, + ], + }, + + ], +}; + +local g2 = { + xplanes: { + // Distances between collection plane and a logical plane. + danode: 0.0, + dresponse: 10*wc.cm, + // Location of cathode measured from collection is based on a dump + // of ubcore/v06_83_00/gdml/microboonev11.gdml by Matt Toups + dcathode: 2.5480*wc.m, + } +}; + +{ + "g1+g2": g1+g2, + "g2+g1": g2+g1, + "merge(g1,g2)": std.mergePatch(g1,g2), + "merge(g2,g1)": std.mergePatch(g2,g1), +} diff --git a/cfg/wirecell.jsonnet b/cfg/wirecell.jsonnet index bc0ed59c5..22a43d77f 100644 --- a/cfg/wirecell.jsonnet +++ b/cfg/wirecell.jsonnet @@ -388,6 +388,11 @@ }, + + // This is std.get from 0.18.0 + get(o, f, default=null, inc_hidden=true):: + if std.objectHasEx(o, f, inc_hidden) then o[f] else default, + } diff --git a/gen/docs/talks/lmn-fr/lmn-fr.org b/gen/docs/talks/lmn-fr/lmn-fr.org new file mode 100644 index 000000000..2b70400a3 --- /dev/null +++ b/gen/docs/talks/lmn-fr/lmn-fr.org @@ -0,0 +1,302 @@ +#+title: LMN on FR +#+setupfile: ~/sync/talks/setup.org +#+LATEX_HEADER: \usepackage{siunitx} + +#+setupfile: "./lmn-fr-pdsp-0.org" + +* Meta :noexport: + +Requires: +- ~wire-cell-python~ installed and in environment. +- ~wire-cell-data~ repo in ~WIRECELL_PATH~. + +Run ~lmn-fr-plots.sh~ to generate ~.pdf~ and ~.org~ files in same directory as this +file and export to beamer PDF (~C-c C-e l P~). + +* Overview + +We want to test LMN resampling on realistic waveforms + +- Need "full simulation" but want to avoid new field calculations. +- Apply LMN resampling to existing field responses (FRs): + +| | | | | +| | Field response | simulation | ADC | +| | sampling | downsampling | sampling | +|-------+----------------+--------------+----------| +| have: | \qty{100}{ns} | \times 5 | 500 ns | +| want: | \qty{64}{ns} | \times 8 | 512 ns | + + +- Understand FR related units and normalization. + +- PDSP FR issues + +- Apply to resampled FR in WCT simulation + + +* Sampling interpretation and resampling normalization $N_s \to N_r$ + +** Instantaneous + +Samples are continuous function values, resampling is *interpolation*, signal strength preserved. + +- Norm: $A_0 = \frac{N_r}{N_s}$ (normalizing on forward $DFT$ is correct) + +** Integrating + +Samples integrate continuous function over sample period, signal sum preserved. + +- Norm: $A_1 = \frac{N_s}{N_r}A_0 = 1$ (normalizing on inverse $DFT^{-1}$ is correct) + +** Energy + +Samples provide "probability amplitude" or "field" values, signal Parseval energy (sum of squares) is preserved. + +- Norm: $A_2 = \sqrt{\frac{N_r}{N_s}}$ (symmetric $DFT$ normalization is correct) + +** Initial discussion details are in the paper draft. + + +* Resampling current field response ($I$) + +** :B_columns: +:PROPERTIES: +:BEAMER_env: columns +:END: + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.6 +:END: + +{{{fig-fr-pdsp-0}}} + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.4 +:END: + + +We interpret $I$ as an *instantaneous electric current*. + +- \small Resampling is an *interpolation*. +- \small Normalize the DFT round trip by the *initial size*. + + +** Note, all frequency spectra here are $|DFT(x)|$ with no normalization applied. + +* Zoom in - time-domain ringing + +** :B_columns: +:PROPERTIES: +:BEAMER_env: columns +:END: + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.6 +:END: + +{{{fig-fr-zoom-pdsp-0}}} + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.4 +:END: + +*Small ringing* in the time-domain. +- \small Due to non-zero power at the original Nyquist frequency (\qty{5}{MHz}). +- Indicates the $I$ is *undersampled* at the original \qty{100}{ns}. +- The eventual convolution with the slower *electronics response* attenuates the ringing. +- Note, frequency-domain wiggles are actual features, not ringing. + +** *We must be cognizant of ringing when applying LMN to ADC waveforms.* + +* Cold electronics response ($E$) + +** :B_columns: +:PROPERTIES: +:BEAMER_env: columns +:END: + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.6 +:END: + +{{{fig-cer-pdsp-0}}} + +#+begin_center +\small $gain=\qty{14}{mV/fC},\ shaping=\qty{2}{us}$ +#+end_center + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.4 +:END: + +$E$ is an *instantaneous* signal. +- \small Resampling is an *interpolation*. + +- \small Gain sets peak value and must be retained by the resampling. + +* Detector response ($I \otimes E$) + +** :B_columns: +:PROPERTIES: +:BEAMER_env: columns +:END: + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.6 +:END: + +{{{fig-fr-er-pdsp-0}}} + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.4 +:END: + + +We say *detector response* is $I \otimes E$. + +- \small Naively, we get "weird units" which are not $[voltage]$. +- \small The sim multiplies by the sample period to give the missing $[time]$. + +What normalization is correct? +- \small Do we resample before we convolve or vice versa? + +* Charge field response resampling + +** :B_columns: +:PROPERTIES: +:BEAMER_env: columns +:END: + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.6 +:END: + +{{{fig-q-pdsp-0}}} + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.4 +:END: + + + +We should have defined detector response as: \[I \to Q = I \cdot T\] +Resampling $Q$ requires an *integral interpretation*. +- \small Time-domain samples are scaled to preserve total integral. +- Trivially, + +\small $[charge]\cdot[voltage/charge]=[voltage]$. + +* Detector voltage response ($Q \otimes E$) + + +** :B_columns: +:PROPERTIES: +:BEAMER_env: columns +:END: + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.6 +:END: + +{{{fig-v-pdsp-0}}} + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.4 +:END: + +Normalization is now clear. +- \small $Q$ is integral-resampled. +- \small $E$ is interpolated-resampled. +- \small $V$ is (effectively) interpolated-resampled. +WCT ~Digitizer~ component will then produce interpolated ADC. + + +* PDSP FR issue: sample period not exactly \qty{100}{ns} + +The file ~dune-garfield-1d565.json.bz2~ seems to be current ProtoDUNE-SP +- \small Is it still best? + +It has $T=\qty{99.998998998999}{ns}$ which fails LMN "rationality condition". +- \small This is not FP round-off or JSON munging (100.0 is exact in JSON). +- \small For here/now, *force* $T = \qty{100}{ns}$. + +Already mentioned that our FRs are in general *undersampled*. + +* PDSP FR issue: time domain is too short + + +** :B_columns: +:PROPERTIES: +:BEAMER_env: columns +:END: + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.6 +:END: + +{{{fig-fr-front-pdsp-0}}} + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.4 +:END: + + +Early response is ignored. + +- \small FR starts low with a down step. +- Non-zero DC component. +- (ringing in the resampled). + +Note, FR sign opposite from ADC. +- \small Induction has positive net signal. +- $\sum Q / \sum |Q| \approx \qty{0.3}{\%}$ + +\scriptsize This problem was recently "discovered" by UVCGAN. + +* WCT ~DepoTransform~ simulation + +Folds multiple responses: +- field, electronics and RC (or RCRC). +- *downsamples* in time domain from ~fr.period~ to ~tick~. +Should naturally handle \qty{64}{ns} FR. + +* TODO + +- [ ] write 64ns FR file +- [ ] "diagonal" track sim with 512ns and 500ns ADC +- [ ] resample 512ns to 500ns and compare to native 500ns +- [ ] sigproc on resampled-500ns and native-500ns + + + +* + +#+begin_center +\Huge $\mathcal{FIN}$ +#+end_center + + +* ICARUS + +I see ~garfield-icarus-fnal-rev2.json.bz2~ mentioned in ~cfg/~ but only +~garfield-icarus-fnal-rev1.json.bz2~ exists in ~wire-cell-data~. + + +# Local Variables: +# eval: (fix-latex-previews) +# End: diff --git a/gen/inc/WireCellGen/GaussianDiffusion.h b/gen/inc/WireCellGen/GaussianDiffusion.h index 4e6b35cf1..6003d90d3 100644 --- a/gen/inc/WireCellGen/GaussianDiffusion.h +++ b/gen/inc/WireCellGen/GaussianDiffusion.h @@ -1,5 +1,5 @@ -#ifndef WIRECELLGEN_GaussianDiffusion -#define WIRECELLGEN_GaussianDiffusion +#ifndef WIRECELLGEN_GAUSSIANDIFFUSION +#define WIRECELLGEN_GAUSSIANDIFFUSION #include "WireCellUtil/Array.h" #include "WireCellUtil/Binning.h" @@ -7,7 +7,6 @@ #include "WireCellIface/IRandom.h" #include -#include namespace WireCell { namespace Gen { diff --git a/gen/inc/WireCellGen/Reframer.h b/gen/inc/WireCellGen/Reframer.h index 965287aa9..ad444e985 100644 --- a/gen/inc/WireCellGen/Reframer.h +++ b/gen/inc/WireCellGen/Reframer.h @@ -1,9 +1,9 @@ /** - A reframer takes makes a "rectangular" frame filled with samples - from the tagged traces of its input. This new frame has exactly - one trace for each channel and each trace is padded to span a - uniform duration. These configuration paramters control how the - new frame is shaped: + A reframer makes a "rectangular" frame filled with samples from the tagged + traces of its input. This new frame has exactly one trace for each pair of + (tag, channel) and each trace is padded to span the channel and time range. + + These configuration paramters control how the new frame is shaped: - tags :: list of trace tags to consider. an empty list means to use all traces. Default: []. diff --git a/gen/inc/WireCellGen/Retagger.h b/gen/inc/WireCellGen/Retagger.h index f466322c0..f787b1246 100644 --- a/gen/inc/WireCellGen/Retagger.h +++ b/gen/inc/WireCellGen/Retagger.h @@ -1,6 +1,6 @@ /* The retagger support frame and trace tag rule sets like FrameFanin * and FrameFanout but also support a "merge" ruleset which can - * combine trace sets. For exaple an element of the tag_rules in + * combine trace sets. For example an element of the tag_rules in * configuration may look like: { type: "Retagger", diff --git a/gen/inc/WireCellGen/WireBoundedDepos.h b/gen/inc/WireCellGen/WireBoundedDepos.h index 7b0e65ec4..e04d3953e 100644 --- a/gen/inc/WireCellGen/WireBoundedDepos.h +++ b/gen/inc/WireCellGen/WireBoundedDepos.h @@ -5,64 +5,63 @@ #include "WireCellIface/IConfigurable.h" #include "WireCellIface/IAnodePlane.h" #include "WireCellUtil/Units.h" +#include "WireCellAux/Logger.h" #include #include -namespace WireCell { +namespace WireCell::Gen { - namespace Gen { + /** + WireBoundedDepos outputs depos based on which wires they + "land". To "land" means to drift antiparallel to the + X-axis. The set of wire regions on which depos may or may + not "land" is configured as an array of objects with three + integer attributes: - /** - WireBoundedDepos outputs depos based on which wires they - "land". To "land" means to drift antiparallel to the - X-axis. The set of wire regions on which depos may or may - not "land" is configured as an array of objects with three - integer attributes: + [{ + plane: , + min: , + max: , - min: , - max: m_pimpos; - private: - IAnodePlane::pointer m_anode; - bool m_accept; - std::vector m_pimpos; + typedef std::tuple wire_bounds_t; + typedef std::vector wire_region_t; + std::vector m_regions; + }; - typedef std::tuple wire_bounds_t; - typedef std::vector wire_region_t; - std::vector m_regions; - }; - } // namespace Gen -} // namespace WireCell +} // namespace WireCell::Gen #endif diff --git a/gen/src/DepoSetFilter.cxx b/gen/src/DepoSetFilter.cxx index 5a8146887..a7487718e 100644 --- a/gen/src/DepoSetFilter.cxx +++ b/gen/src/DepoSetFilter.cxx @@ -33,9 +33,11 @@ void DepoSetFilter::configure(const WireCell::Configuration& cfg) if (anode == nullptr) { THROW(ValueError() << errmsg{"Input anode is a nullptr"}); } - IAnodeFace::vector abode_faces = anode->faces(); - for (auto face : abode_faces) { - m_boxes.push_back(face->sensitive()); + IAnodeFace::vector anode_faces = anode->faces(); + for (auto face : anode_faces) { + auto bb = face->sensitive(); + m_boxes.push_back(bb); + log->debug("face={} bounds={}", face->ident(), bb.bounds()); } } @@ -43,25 +45,25 @@ bool DepoSetFilter::operator()(const input_pointer& in, output_pointer& out) { out = nullptr; if (!in) { - log->debug("DepoSetFilter fail with no input on call = {}", m_count); + log->debug("EOS at call={}", m_count); return true; } IDepo::vector output_depos; - for (auto idepo : *(in->depos())) { - bool pass = false; + auto indepos = in->depos(); + for (auto idepo : *indepos) { for (auto box : m_boxes) { WireCell::Ray r = box.bounds(); - - if (box.inside(idepo->pos())) { - pass = true; - break; + if (! box.inside(idepo->pos())) { + continue; } - } - if (pass) { output_depos.push_back(idepo); + break; } } - log->debug("call={} Number of Depos for a give APA={}", m_count, output_depos.size()); + const int nin = indepos->size(); + const int nout = output_depos.size(); + log->debug("call={} depos: input={} output={} dropped={}", + m_count, nin, nout, nout-nin); out = std::make_shared(m_count, output_depos); ++m_count; return true; diff --git a/gen/src/EmpiricalNoiseModel.cxx b/gen/src/EmpiricalNoiseModel.cxx index 8dcb1851a..0b85d4af8 100644 --- a/gen/src/EmpiricalNoiseModel.cxx +++ b/gen/src/EmpiricalNoiseModel.cxx @@ -93,8 +93,8 @@ WireCell::Configuration Gen::EmpiricalNoiseModel::default_configuration() const void Gen::EmpiricalNoiseModel::resample(NoiseSpectrum& spectrum) const { - log->debug("m_nsamples={} fft_length={} m_period={} spec: nsamples={} period={}", - m_nsamples, m_fft_length, m_period, spectrum.nsamples, spectrum.period); + // log->debug("m_nsamples={} fft_length={} m_period={} spec: nsamples={} period={}", + // m_nsamples, m_fft_length, m_period, spectrum.nsamples, spectrum.period); if (spectrum.nsamples == m_fft_length && spectrum.period == m_period) { return; // natural sampling is what is requested. @@ -221,8 +221,8 @@ void Gen::EmpiricalNoiseModel::configure(const WireCell::Configuration& cfg) resample(*nsptr); m_spectral_data[nsptr->plane].push_back(nsptr); // assumes ordered by wire length! - log->debug("nwanted={} plane={} ntold={} ngot={} ninput={}", - m_nsamples, nsptr->plane, nsptr->nsamples, nsptr->amps.size(), nfreqs); + // log->debug("nwanted={} plane={} ntold={} ngot={} ninput={}", + // m_nsamples, nsptr->plane, nsptr->nsamples, nsptr->amps.size(), nfreqs); } if (errors) { diff --git a/gen/src/WireBoundedDepos.cxx b/gen/src/WireBoundedDepos.cxx index 74e922ae3..b72928f0c 100644 --- a/gen/src/WireBoundedDepos.cxx +++ b/gen/src/WireBoundedDepos.cxx @@ -9,20 +9,28 @@ #include -#include -#include -WIRECELL_FACTORY(WireBoundedDepos, WireCell::Gen::WireBoundedDepos, WireCell::IDrifter, WireCell::IConfigurable) +WIRECELL_FACTORY(WireBoundedDepos, + WireCell::Gen::WireBoundedDepos, + WireCell::INamed, WireCell::IDrifter, WireCell::IConfigurable) using namespace std; using namespace WireCell; +Gen::WireBoundedDepos::WireBoundedDepos() + : Aux::Logger("WireBoundedDepos", "gen") +{ +} +Gen::WireBoundedDepos::~WireBoundedDepos() +{ +} + WireCell::Configuration Gen::WireBoundedDepos::default_configuration() const { Configuration cfg; cfg["anode"] = ""; - // A the list of wire regions + // List of wire regions // [{ // plane: , // min: , @@ -43,7 +51,6 @@ void Gen::WireBoundedDepos::configure(const WireCell::Configuration& cfg) for (auto face : m_anode->faces()) { if (face->planes().empty()) { - std::cerr << "Gen::WireBoundedDepos: not given multi-plane AnodeFace for face " << face->ident() << "\n"; continue; } for (auto plane : face->planes()) { @@ -55,7 +62,7 @@ void Gen::WireBoundedDepos::configure(const WireCell::Configuration& cfg) } break; // fixme: } - m_accept = cfg["mode"].asString() == "accept"; + m_accept = get(cfg, "mode", "accept") == "accept"; auto jregions = cfg["regions"]; for (auto jregion : jregions) { @@ -66,8 +73,7 @@ void Gen::WireBoundedDepos::configure(const WireCell::Configuration& cfg) m_regions.push_back(wr); } - std::cerr << "WireBoundedDepos: " << cfg["mode"] << " with " << m_regions.size() << " wires in " << m_pimpos.size() - << " planes\n"; + log->debug("accept={} nregions={} nplanes={}", m_accept, m_regions.size(), m_pimpos.size()); } bool Gen::WireBoundedDepos::operator()(const input_pointer& depo, output_queue& outq) @@ -117,8 +123,3 @@ bool Gen::WireBoundedDepos::operator()(const input_pointer& depo, output_queue& return true; } -Gen::WireBoundedDepos::WireBoundedDepos() - : m_accept(true) -{ -} -Gen::WireBoundedDepos::~WireBoundedDepos() {} diff --git a/gen/test/test-lmn-fr-pdsp.bats b/gen/test/test-lmn-fr-pdsp.bats new file mode 100644 index 000000000..495fbdb9a --- /dev/null +++ b/gen/test/test-lmn-fr-pdsp.bats @@ -0,0 +1,211 @@ +#!/usr/bin/env bats + +# Run things for LMN FR resampling test for detector: +DETECTOR=pdsp + +# The field sampling period in ns which we want. +FIELDS_TDE_TICK=100 + +# ADC sampling period in ns that we want. +ADC_TDE_TICK=500 + +# The field sampling period in ns that we will resample. +FIELDS_BDE_TICK=64 + +# Original ADC sampling period in ns that we will resample. +ADC_BDE_TICK=512 + + + +bats_load_library wct-bats.sh + + +setup_file () { + cd_tmp + + # prime idempotent running + cat < depo-line-args-tmp +depo-line \ +--electron-density -5000 \ +--first -3578.36*mm,76.1*mm,3.3497*mm \ +--last -1.5875*mm,6060*mm,2303.37*mm \ +--sigma 1*mm,1*mm \ +--output depos.npz +EOF + mv_if_diff depo-line-args-tmp depo-line-args + + cat < preplots-args-tmp +lmn-fr-plots \ +--period 100*ns \ +--zoom-start 50*us \ +--plane 0 \ +--detector-name pdsp \ +-o lmn-fr-pdsp-0.pdf \ +-O lmn-fr-pdsp-0.org +EOF + mv_if_diff preplots-args-tmp preplots-args +} + + +@test "lmn fr preplots" { + cd_tmp file + + read -a args < preplots-args + + run_idempotently -s preplots-args -t lmn-fr-pdsp-0.pdf -t lmn-fr-pdsp-0.org -- \ + wirecell-resp "${args[@]}" + +} + + +fields_original="$(wirecell-util resolve -k fields "$DETECTOR" | head -1)" +fields_name="$(basename "$fields_original" .json.bz2)" + +# These file names MUST match those set by the "resample_*" PDSP variants. +fields_tde="${fields_name}-${FIELDS_TDE_TICK}ns.json.bz2" +fields_bde="${fields_name}-${FIELDS_BDE_TICK}ns.json.bz2" + + +@test "lmn fr resample fr" { + cd_tmp file + + # Produce field file names that match the resample_* PDSP variants. + + # default PDSP fr.period is not exactly 100ns so we rewrite it with 100ns + # forced and also put the result locally. + run_idempotently -s "$files_original" -t "$fields_sampled" -- \ + wirecell-resp condition \ + -P "$FIELDS_TDE_TICK"'*ns' \ + -o "$fields_tde" \ + "$fields_original" + + # Resample 100ns TDE fields to construct a 64ns BDE fields + run_idempotently -s "$fields_tde" -t "$fields_bde" -- \ + wirecell-resp resample \ + -t "$FIELDS_BDE_TICK"'*ns' \ + -o "$fields_bde" \ + "$fields_tde" +} + +@test "lmn fr gen depos" { + cd_tmp file + + read -a args < depo-line-args + + run_idempotently -s depo-line-args -t depos.npz -- \ + wirecell-gen "${args[@]}" +} + +tier_name () { + local tier="$1" ; shift + local tick="$1" ; shift + local ext="${1:-npz}" + echo "${DETECTOR}-${tier}-${tick}ns.${ext}" +} + +@test "lmn fr sim" { + cd_tmp file + + # We simulate with at 100ns->500ns and 64ns->512ns. + for tick in $ADC_TDE_TICK $ADC_BDE_TICK + do + log="$(tier_name wct-sim $tick log)" + out="$(tier_name digits $tick)" + rm -f "$log" + + run_idempotently -s depos.npz -t "$log" -t "$out" -- \ + wire-cell -l "$log" -L debug \ + -A input="depos.npz" \ + --tla-code output="{sim:\"$out\"}" \ + -A detector="${DETECTOR}" \ + -A variant="resample_${tick}ns" \ + -A tasks="drift,sim" \ + layers/omnijob.jsonnet + done +} + +@test "lmn fr resample digits" { + cd_tmp file + + # We symlink the TDE digits to a new name to match later pattern + ln -sf \ + $(tier_name digits $ADC_TDE_TICK) \ + $(tier_name digits-$ADC_TDE_TICK $ADC_TDE_TICK) + + + # Resample BDE->TDE. Note, eventually this should be done inside NF in + # order to save one FFT round trip per channel. + digin=$(tier_name digits $ADC_BDE_TICK) + digout=$(tier_name digits-$ADC_BDE_TICK $ADC_TDE_TICK) + run_idempotently -s $digin -t $digout -- \ + wirecell-util resample \ + -t "$ADC_TDE_TICK"'*ns' \ + -o $digout \ + $digin +} + +select_channels_near='335:350,1500:1548,2080:2100' +select_channels_far='0:20,1200:1220,2540:2560' + +signal_trange_near='175*us,300*us' +signal_trange_far='2.3*ms,2.5*ms' +digit_trange_near='175*us,300*us' +digit_trange_far='2.3*ms,2.5*ms' + +@test "lmn fr adc plots" { + cd_tmp file + + dig1=$(tier_name digits $ADC_TDE_TICK) + dig2=$(tier_name digits $ADC_BDE_TICK) + dig3=$(tier_name digits-$ADC_BDE_TICK $ADC_TDE_TICK) + + out=channels-adc.pdf + run_idempotently -s "$dig1" -s "$dig2" -s "$dig3" -t "$out" -- \ + wirecell-plot channels \ + -o "$out" \ + --trange '300*us,400*us' \ + --channel 335,1500,2100 \ + "$dig1" "$dig2" "$dig3" + +} + + +@test "lmn fr sigproc" { + cd_tmp file + + local name=pdsp + + # Do normal signal processing using normal tick on normally sampled and + # resampled digits + for tick in $ADC_TDE_TICK $ADC_BDE_TICK + do + log="$(tier_name wct-sp-$tick $ADC_TDE_TICK log)" + digin=$(tier_name digits-$tick $ADC_TDE_TICK) + sigout=$(tier_name signals-$tick $ADC_TDE_TICK) + run_idempotently -s "$digin" -t "$log" -t "$sigout" -- \ + wire-cell -l "$log" -L debug \ + -A input="$digin" \ + --tla-code output='{sp:"'$sigout'"}' \ + -A detector="${DETECTOR}" \ + -A variant="resample_${ADC_TDE_TICK}ns" \ + -A tasks="sp" \ + layers/omnijob.jsonnet + done + +} + +@test "lmn fr sig plots" { + cd_tmp file + + sig1=$(tier_name signals-$ADC_TDE_TICK $ADC_TDE_TICK) + sig2=$(tier_name signals-$ADC_BDE_TICK $ADC_TDE_TICK) + + out=channels-sig.pdf + run_idempotently -s "$sig1" -s "$sig2" -t "$out" -- \ + wirecell-plot channels \ + -o "$out" \ + --trange '250*us,300*us' \ + --channel 335,1500,2100 \ + "$sig1" "$sig2" + +} diff --git a/sigproc/inc/WireCellSigProc/FrameMerger.h b/sigproc/inc/WireCellSigProc/FrameMerger.h index 8c1e2f17c..130ec5d58 100644 --- a/sigproc/inc/WireCellSigProc/FrameMerger.h +++ b/sigproc/inc/WireCellSigProc/FrameMerger.h @@ -8,25 +8,26 @@ #include "WireCellIface/IConfigurable.h" #include "WireCellIface/IFrameJoiner.h" #include "WireCellUtil/Configuration.h" +#include "WireCellAux/Logger.h" -namespace WireCell { - namespace SigProc { - class FrameMerger : public IFrameJoiner, public IConfigurable { - public: - FrameMerger(); - virtual ~FrameMerger(); - - // IJoinNode - virtual bool operator()(const input_tuple_type& intup, output_pointer& out); - - // IConfigurable - virtual void configure(const WireCell::Configuration& config); - virtual WireCell::Configuration default_configuration() const; - - private: - Configuration m_cfg; - }; - } // namespace SigProc -} // namespace WireCell +namespace WireCell::SigProc { + class FrameMerger : public Aux::Logger, public IFrameJoiner, public IConfigurable { + public: + FrameMerger(); + virtual ~FrameMerger(); + + // IJoinNode + virtual bool operator()(const input_tuple_type& intup, output_pointer& out); + + // IConfigurable + virtual void configure(const WireCell::Configuration& config); + virtual WireCell::Configuration default_configuration() const; + + private: + Configuration m_cfg; + size_t m_count{0}; + }; + +} // namespace WireCell::SigProc #endif diff --git a/sigproc/inc/WireCellSigProc/FrameSplitter.h b/sigproc/inc/WireCellSigProc/FrameSplitter.h index 1a201e636..637c2fa78 100644 --- a/sigproc/inc/WireCellSigProc/FrameSplitter.h +++ b/sigproc/inc/WireCellSigProc/FrameSplitter.h @@ -2,21 +2,23 @@ #define WIRECELLSIGPROC_FRAMESPLITTER #include "WireCellIface/IFrameSplitter.h" +#include "WireCellAux/Logger.h" -namespace WireCell { - namespace SigProc { +namespace WireCell::SigProc { - // fixme/note: this class is generic and nothing to do with - // sigproc per se. There are a few such frame related nodes - // also in gen. They should be moved into some neutral package. - class FrameSplitter : public WireCell::IFrameSplitter { - public: - FrameSplitter(); - virtual ~FrameSplitter(); + // fixme/note: this class is generic and nothing to do with + // sigproc per se. There are a few such frame related nodes + // also in gen. They should be moved into some neutral package. + class FrameSplitter : public Aux::Logger, public WireCell::IFrameSplitter { + public: + FrameSplitter(); + virtual ~FrameSplitter(); - virtual bool operator()(const input_pointer& in, output_tuple_type& out); - }; - } // namespace SigProc -} // namespace WireCell + virtual bool operator()(const input_pointer& in, output_tuple_type& out); + private: + size_t m_count{0}; + }; + +} // namespace WireCell::SigProc #endif diff --git a/sigproc/inc/WireCellSigProc/L1SPFilter.h b/sigproc/inc/WireCellSigProc/L1SPFilter.h index c7e7f4e3c..e71a3763f 100644 --- a/sigproc/inc/WireCellSigProc/L1SPFilter.h +++ b/sigproc/inc/WireCellSigProc/L1SPFilter.h @@ -11,13 +11,14 @@ #include "WireCellIface/IDFT.h" #include "WireCellAux/SimpleTrace.h" +#include "WireCellAux/Logger.h" #include "WireCellUtil/Interpolate.h" namespace WireCell { namespace SigProc { - class L1SPFilter : public IFrameFilter, public IConfigurable { + class L1SPFilter : public Aux::Logger, public IFrameFilter, public IConfigurable { public: L1SPFilter(double gain = 14.0 * units::mV / units::fC, double shaping = 2.2 * units::microsecond, double postgain = 1.2, double ADC_mV = 4096 / (2000. * units::mV), @@ -52,6 +53,8 @@ namespace WireCell { linterp* lin_V; linterp* lin_W; + + size_t m_count{0}; }; } // namespace SigProc } // namespace WireCell diff --git a/sigproc/src/FrameMerger.cxx b/sigproc/src/FrameMerger.cxx index 1647f6456..225cd4d75 100644 --- a/sigproc/src/FrameMerger.cxx +++ b/sigproc/src/FrameMerger.cxx @@ -7,7 +7,9 @@ #include -WIRECELL_FACTORY(FrameMerger, WireCell::SigProc::FrameMerger, WireCell::IFrameJoiner, WireCell::IConfigurable) +WIRECELL_FACTORY(FrameMerger, + WireCell::SigProc::FrameMerger, + WireCell::INamed, WireCell::IFrameJoiner, WireCell::IConfigurable) using namespace WireCell; @@ -15,16 +17,16 @@ Configuration SigProc::FrameMerger::default_configuration() const { Configuration cfg; - // The merge map specifies a list of triples of tags. The first - // tag in the pair corresponds to a trace tag on the first frame, - // and likewise the second. The third is a tag that will be - // placed on the merged set of traces in the output frame. The - // merge is performed considering only traces in the two input - // frames that respectively match their tags. As tags may have - // overlapping sets of traces, it is important to recognize that - // the merge is progressive and in order of the merge map. If the - // merge map is empty then all traces in frame 1 are merged with - // all traces in frame 2 irrespective of any tags. + // The merge map specifies a list of triples of tags. The first (second) + // tag is compared to trace tags in the frame from port 0 (1). Either or + // both of the first two tags may be empty strings which match untagged + // traces. The third tag placed on the merged set of traces from these two + // input in the output frame. The merge is performed considering only + // traces in the two input frames that respectively match their tags. As + // tags may have overlapping sets of traces, it is important to recognize + // that the merge is progressive and in order of the merge map. If the + // merge map is empty then all traces in frame 1 are merged with all traces + // in frame 2 irrespective of any tags. cfg["mergemap"] = Json::arrayValue; // The rule determins the algorithm employed in the merge. @@ -49,23 +51,24 @@ bool SigProc::FrameMerger::operator()(const input_tuple_type& intup, output_poin auto one = std::get<0>(intup); auto two = std::get<1>(intup); if (!one or !two) { - std::cerr << "FrameMerger: EOS\n"; + log->debug("EOS at call={}", m_count++); return true; } auto jmergemap = m_cfg["mergemap"]; const int nsets = jmergemap.size(); + log->debug("call={} frame1: {}", m_count, WireCell::Aux::taginfo(one)); + log->debug("call={} frame2: {}", m_count, WireCell::Aux::taginfo(two)); + // collect traces into a vector of vector whether we are dealling // with all traces or honoring tags. std::vector tracesv1, tracesv2; if (!nsets) { tracesv1.push_back(Aux::untagged_traces(one)); tracesv2.push_back(Aux::untagged_traces(two)); - std::cerr << "FrameMerger: see frame: " << one->ident() << " no tags, whole frame\n"; } else { - std::cerr << "FrameMerger: see frame: " << one->ident() << " with tags:\n"; for (int ind = 0; ind < nsets; ++ind) { auto jtags = jmergemap[ind]; std::string tag1 = jtags[0].asString(); @@ -73,9 +76,10 @@ bool SigProc::FrameMerger::operator()(const input_tuple_type& intup, output_poin std::string tag3 = jtags[2].asString(); tracesv1.push_back(Aux::tagged_traces(one, tag1)); tracesv2.push_back(Aux::tagged_traces(two, tag2)); - std::cerr << "\ttags: " << tag1 << "[" << tracesv1.back().size() << "]" - << " + " << tag2 << "[" << tracesv2.back().size() << "]" - << " -> " << tag3 << "\n"; + log->debug("call={} tags: {}[{}] + {}[{}] -> {}", + m_count, + tag1, tracesv1.back().size(), + tag2, tracesv2.back().size(), tag3); } } @@ -136,6 +140,8 @@ bool SigProc::FrameMerger::operator()(const input_tuple_type& intup, output_poin return true; } -SigProc::FrameMerger::FrameMerger() {} +SigProc::FrameMerger::FrameMerger() + : Aux::Logger("FrameMerger", "sigproc") +{} SigProc::FrameMerger::~FrameMerger() {} diff --git a/sigproc/src/FrameSplitter.cxx b/sigproc/src/FrameSplitter.cxx index 6c459af09..db31ef93a 100644 --- a/sigproc/src/FrameSplitter.cxx +++ b/sigproc/src/FrameSplitter.cxx @@ -3,31 +3,30 @@ #include "WireCellUtil/NamedFactory.h" #include "WireCellAux/FrameTools.h" -#include -WIRECELL_FACTORY(FrameSplitter, WireCell::SigProc::FrameSplitter, WireCell::IFrameSplitter) +WIRECELL_FACTORY(FrameSplitter, + WireCell::SigProc::FrameSplitter, + WireCell::INamed, WireCell::IFrameSplitter) using namespace WireCell::SigProc; -FrameSplitter::FrameSplitter() {} +FrameSplitter::FrameSplitter() + : Aux::Logger("FrameSplitter", "sigproc") +{ +} + FrameSplitter::~FrameSplitter() {} bool FrameSplitter::operator()(const input_pointer& in, output_tuple_type& out) { - if (!in) { - std::cerr << "FrameSplitter: passing on EOS\n"; - } - else { - std::cerr << "FrameSplitter: passing on frame: " << in->ident() << ":"; - for (auto tag : in->trace_tags()) { - auto tt = Aux::tagged_traces(in, tag); - std::cerr << " " << tag << "[" << tt.size() << "]"; - } - std::cerr << std::endl; - } - get<0>(out) = in; get<1>(out) = in; + if (!in) { + log->debug("EOS at call={}", m_count++); + return true; + } + + log->debug("call={} {}", m_count++, WireCell::Aux::taginfo(in)); return true; } diff --git a/sigproc/src/L1SPFilter.cxx b/sigproc/src/L1SPFilter.cxx index ac6d58978..b80f5b12c 100644 --- a/sigproc/src/L1SPFilter.cxx +++ b/sigproc/src/L1SPFilter.cxx @@ -14,7 +14,9 @@ #include #include -WIRECELL_FACTORY(L1SPFilter, WireCell::SigProc::L1SPFilter, WireCell::IFrameFilter, WireCell::IConfigurable) +WIRECELL_FACTORY(L1SPFilter, + WireCell::SigProc::L1SPFilter, + WireCell::INamed, WireCell::IFrameFilter, WireCell::IConfigurable) using namespace Eigen; using namespace WireCell; @@ -26,7 +28,8 @@ using WireCell::Aux::DftTools::inv_c2r; L1SPFilter::L1SPFilter(double gain, double shaping, double postgain, double ADC_mV, double fine_time_offset, double coarse_time_offset) - : m_gain(gain) + : Aux::Logger("L1SPFilter", "sigproc") + , m_gain(gain) , m_shaping(shaping) , m_postgain(postgain) , m_ADC_mV(ADC_mV) @@ -82,8 +85,6 @@ void L1SPFilter::init_resp() // convolute with V and Y average responses ... double intrinsic_time_offset = fravg.origin / fravg.speed; - // std::cout << intrinsic_time_offset << " " << m_fine_time_offset << " " << m_coarse_time_offset << " " << - // m_gain << " " << 14.0 * units::mV/units::fC << " " << m_shaping << " " << fravg.period << std::endl; double x0 = (-intrinsic_time_offset - m_coarse_time_offset + m_fine_time_offset); double xstep = fravg.period; @@ -183,6 +184,7 @@ bool L1SPFilter::operator()(const input_pointer& in, output_pointer& out) { out = nullptr; if (!in) { + log->debug("EOS at call={}", m_count++); return true; // eos } @@ -190,8 +192,6 @@ bool L1SPFilter::operator()(const input_pointer& in, output_pointer& out) std::string sigtag = get(m_cfg, "sigtag"); std::string outtag = get(m_cfg, "outtag"); - // std::cout << smearing_vec.size() << std::endl; - int roi_pad = 0; roi_pad = get(m_cfg, "roi_pad", roi_pad); int raw_pad = 0; @@ -200,26 +200,19 @@ bool L1SPFilter::operator()(const input_pointer& in, output_pointer& out) double raw_ROI_th_nsigma = get(m_cfg, "raw_ROI_th_nsigma", 4); double raw_ROI_th_adclimit = get(m_cfg, "raw_ROI_th_adclimit", 10); - // std::cout << "Xin: " << raw_ROI_th_nsigma << " " << raw_ROI_th_adclimit << " " << overall_time_offset << " " << - // collect_time_offset << " " << roi_pad << " " << adc_l1_threshold << " " << adc_sum_threshold << " " << - // adc_sum_rescaling << " " << adc_sum_rescaling_limit << " " << l1_seg_length << " " << l1_scaling_factor << " " << - // l1_lambda << " " << l1_epsilon << " " << l1_niteration << " " << l1_decon_limit << " " << l1_resp_scale << " " << - // l1_col_scale << " " << l1_ind_scale << std::endl; init_resp(); auto adctraces = Aux::tagged_traces(in, adctag); auto sigtraces = Aux::tagged_traces(in, sigtag); if (adctraces.empty() or sigtraces.empty() or adctraces.size() != sigtraces.size()) { - std::cerr << "L1SPFilter got unexpected input: " << adctraces.size() << " ADC traces and " << sigtraces.size() - << " signal traces\n"; - THROW(RuntimeError() << errmsg{"L1SPFilter: unexpected input"}); + log->error("unexpected input: {} ADC traces, {} signal traces at call={}", + adctraces.size(), sigtraces.size(), m_count++); + raise("L1SPFilter: unexpected input"); } m_period = in->tick(); - // std::cout << m_period/units::us << std::endl; - /// here, use the ADC and signal traces to do L1SP /// put result in out_traces ITrace::vector out_traces; @@ -240,9 +233,6 @@ bool L1SPFilter::operator()(const input_pointer& in, output_pointer& out) } init_map[ch] = time_ticks; - // if (time_ticks.size()>0){ - // std::cout << ch << " " << time_ticks.size() << std::endl; - // } } // do ROI from the raw signal @@ -269,9 +259,6 @@ bool L1SPFilter::operator()(const input_pointer& in, output_pointer& out) double cut = raw_ROI_th_nsigma * sqrt((pow(mean_p1sig - mean, 2) + pow(mean_n1sig - mean, 2)) / 2.); if (cut < raw_ROI_th_adclimit) cut = raw_ROI_th_adclimit; - // if (ch==4090) - // std::cout << cut << " " << raw_pad << std::endl; - for (int qi = 0; qi < ntbins; qi++) { if (fabs(charges[qi]) > cut) { for (int qii = -raw_pad; qii != raw_pad + 1; qii++) { @@ -279,9 +266,6 @@ bool L1SPFilter::operator()(const input_pointer& in, output_pointer& out) } } } - // if (time_ticks.size()>0){ - // std::cout << ch << " " << time_ticks.size() << std::endl; - // } } // create ROIs ... @@ -332,9 +316,6 @@ bool L1SPFilter::operator()(const input_pointer& in, output_pointer& out) } if (rois_save.size() > 0) map_ch_rois[wire_index] = rois_save; - // for (auto it = rois_save.begin(); it!=rois_save.end(); it++){ - // std::cout << wire_index << " " << it->first << " " << it->second +1 << std::endl; - // } } std::map> map_ch_flag_rois; @@ -407,7 +388,6 @@ bool L1SPFilter::operator()(const input_pointer& in, output_pointer& out) } } - // std::cout << trace->channel() << std::endl; out_traces.push_back(newtrace); } @@ -457,18 +437,21 @@ bool L1SPFilter::operator()(const input_pointer& in, output_pointer& out) out_traces.at(i2) = new_trace; } - std::cerr << "L1SPFilter: frame: " << in->ident() << " " << adctag << "[" << adctraces.size() << "] + " << sigtag - << "[" << sigtraces.size() << "] --> " << outtag << "[" << out_traces.size() << "]\n"; - - // Finally, we save the traces to an output frame with tags. - - IFrame::trace_list_t tl(out_traces.size()); - std::iota(tl.begin(), tl.end(), 0); + // Finally, we save the traces to an output frame. auto sf = new Aux::SimpleFrame(in->ident(), in->time(), out_traces, in->tick()); - sf->tag_traces(outtag, tl); + if (! outtag.empty()) { + IFrame::trace_list_t tl(out_traces.size()); + std::iota(tl.begin(), tl.end(), 0); + sf->tag_traces(outtag, tl); + } out = IFrame::pointer(sf); + log->debug("call={} adctag={} sigtag={} outtag={}", m_count, adctag, sigtag, outtag); + log->debug("call={} in frame: {}", m_count, Aux::taginfo(in)); + log->debug("call={} out frame: {}", m_count, Aux::taginfo(out)); + ++m_count; + return true; } @@ -542,14 +525,6 @@ int L1SPFilter::L1_fit(std::shared_ptr& newtrace, flag_l1 = 2; } - // if (adctrace->channel() == 4079){ - // std::cout << adctrace->channel() << " " << nbin_fit << " " << start_tick << " " << end_tick << " " << temp_sum << - // " " << temp1_sum << " " << temp2_sum << " " << max_val << " " << min_val << " " << flag_l1 << std::endl; - // } - - // std::cout << temp_sum << " " << temp1_sum << " " << temp_sum/(temp1_sum*adc_sum_rescaling*1.0/nbin_fit) << " " << - // adc_ratio_threshold << " " << temp1_sum*adc_sum_rescaling*1.0/nbin_fit << " " << flag_l1 << std::endl; - if (flag_l1 == 1) { // do L1 fit ... int n_section = std::round(nbin_fit / l1_seg_length * 1.0); @@ -612,7 +587,7 @@ int L1SPFilter::L1_fit(std::shared_ptr& newtrace, l1_signal.at(j) = final_beta(j) * l1_col_scale + final_beta(nbin_fit + j) * l1_ind_scale; } int mid_bin = (smearing_vec.size() - 1) / 2; - // std::cout << smearing_vec.size() << " " << mid_bin << std::endl; + for (int j = 0; j != nbin_fit; j++) { double content = l1_signal.at(j); if (content > 0) { @@ -676,15 +651,8 @@ int L1SPFilter::L1_fit(std::shared_ptr& newtrace, l1_signal.at(k) = 0; } } - // std::cout << max_val << " " << mean_val1 << std::endl; } - // std::cout << nonzero_bins.front() << " X " << nonzero_bins.back() << std::endl; - // std::cout << ROIs.size() << std::endl; - // for (size_t i=0;i!=ROIs.size();i++){ - // std::cout << ROIs.at(i).first << " " << ROIs.at(i).second << std::endl; - // } - // finish cleaning ... } diff --git a/sigproc/src/Microboone.cxx b/sigproc/src/Microboone.cxx index 3289b3ad8..0ba34b07c 100644 --- a/sigproc/src/Microboone.cxx +++ b/sigproc/src/Microboone.cxx @@ -906,6 +906,7 @@ WireCell::Waveform::ChannelMaskMap Microboone::OneChannelNoise::apply(int ch, si // sanity check data/config match. const size_t nsiglen = signal.size(); int nmismatchlen = 0; + static bool warn_once{true}; // fixme: some channels are just bad can should be skipped. @@ -978,10 +979,11 @@ WireCell::Waveform::ChannelMaskMap Microboone::OneChannelNoise::apply(int ch, si } } - if (nmismatchlen) { + if (warn_once && nmismatchlen) { std::cerr << "OneChannelNoise: WARNING: " << nmismatchlen << " config/data mismatches. " << "#spec=" << nspec << ", #wave=" << nsiglen << ".\n" - << "\tResults may be suspect." << std::endl; + << "\tResults may be suspect. Only one warning given but there may be many suppressed, one per channel" << std::endl; + warn_once = false; } // remove the DC component diff --git a/sigproc/src/OmnibusNoiseFilter.cxx b/sigproc/src/OmnibusNoiseFilter.cxx index 5f94cb37f..ea62fbcbd 100644 --- a/sigproc/src/OmnibusNoiseFilter.cxx +++ b/sigproc/src/OmnibusNoiseFilter.cxx @@ -242,8 +242,10 @@ bool OmnibusNoiseFilter::operator()(const input_pointer& inframe, output_pointer for (size_t ind = 0; ind < itraces.size(); ++ind) { indices[ind] = ind; } - sframe->tag_traces(m_outtag, indices); - sframe->tag_frame("noisefilter"); + if (! m_outtag.empty()) { + sframe->tag_traces(m_outtag, indices); + } + sframe->tag_frame("noisefilter"); // fixme: this is unnecessary and should be removed outframe = IFrame::pointer(sframe); log->debug("call={} output frame: {}", m_count, Aux::taginfo(outframe)); diff --git a/sigproc/src/OmnibusPMTNoiseFilter.cxx b/sigproc/src/OmnibusPMTNoiseFilter.cxx index 40ccdfca7..274f32a62 100644 --- a/sigproc/src/OmnibusPMTNoiseFilter.cxx +++ b/sigproc/src/OmnibusPMTNoiseFilter.cxx @@ -224,13 +224,15 @@ bool OmnibusPMTNoiseFilter::operator()(const input_pointer& in, output_pointer& itraces.push_back(ITrace::pointer(trace)); } - IFrame::trace_list_t indices(itraces.size()); - for (size_t ind = 0; ind < itraces.size(); ++ind) { - indices[ind] = ind; + auto sframe = new Aux::SimpleFrame(in->ident(), in->time(), itraces, in->tick(), in->masks()); + if (! m_outtag.empty()) { + IFrame::trace_list_t indices(itraces.size()); + for (size_t ind = 0; ind < itraces.size(); ++ind) { + indices[ind] = ind; + } + sframe->tag_traces(m_outtag, indices); } - auto sframe = new Aux::SimpleFrame(in->ident(), in->time(), itraces, in->tick(), in->masks()); - sframe->tag_traces(m_outtag, indices); out = IFrame::pointer(sframe); return true; diff --git a/sio/inc/WireCellSio/FrameFileSource.h b/sio/inc/WireCellSio/FrameFileSource.h index 19acde7ce..b8a10693e 100644 --- a/sio/inc/WireCellSio/FrameFileSource.h +++ b/sio/inc/WireCellSio/FrameFileSource.h @@ -41,14 +41,21 @@ namespace WireCell::Sio { /** Config: "tags". - A set of tags to match against the tag portion of the .npy - file names to determine the array should be input. If the - set is empty or a tag "*" is given then all arrays with - the current frame ident are accepted. All matched arrays - of types {frame, channel, tickinfo} are interpreted as - providing tagged traces and the frame and channel arrays - are appended to the IFrame traces and channels - collections. + A set of tags selecting "framelet" array to include as tagged + traces. + + If the portion of the framelet array file name, eg + + frame__ + + is found in this set the framelet will provide tagged traces with + the tag . + + Note, the use of a that is empty "" or the special "*" is + intended for use only for the context of a frame file to indicate + untagged traces. In the context of an IFrame, it is against + convention to explicitly tag traces with either label. + */ std::vector m_tags; @@ -62,7 +69,11 @@ namespace WireCell::Sio { boost::iostreams::filtering_istream m_in; IFrame::pointer load(); - bool matches(const std::string& tag); + + // Classify a label from a framelet array name. + bool is_tagged(const std::string& tag); + bool is_untagged(const std::string& tag); + bool is_excluded(const std::string& tag); size_t m_count{0}; bool m_eos_sent{false}; diff --git a/sio/src/DepoFileSource.cxx b/sio/src/DepoFileSource.cxx index 7ac1f6a69..7d199c034 100644 --- a/sio/src/DepoFileSource.cxx +++ b/sio/src/DepoFileSource.cxx @@ -181,10 +181,13 @@ IDepoSet::pointer Sio::DepoFileSource::next() return nullptr; } - std::vector sdepos; + // build list of depos in order given + std::vector> all_sdepos; + + // First pass make all depos for (size_t ind=0; ind < ndepos; ++ind) { - auto sdepo = new Aux::SimpleDepo( + auto sdepo = std::make_shared( darr(ind, 0), // t Point(darr(ind, 2), // x darr(ind, 3), // y @@ -196,37 +199,48 @@ IDepoSet::pointer Sio::DepoFileSource::next() iarr(ind, 0), // id iarr(ind, 1)); // pdg + all_sdepos.push_back(sdepo); + } + + // second pass, save out the active and resolve the prior depos + std::vector> sdepos; + size_t npriors = 0; + size_t npriors_missing = 0; + for (size_t ind=0; ind < ndepos; ++ind) { + const auto gen = iarr(ind, 2); - if (gen > 0) { - // this depo is a prior - const size_t other = iarr(ind, 3); - if (other >= sdepos.size()) { - log->warn("call={}, prior depo {} not provided in {}", - m_count, other, sdepos.size()); - } - else { - auto idepo = IDepo::pointer(sdepo); - sdepo = nullptr; - sdepos[other]->set_prior(idepo); - } + if (!gen) { // active + sdepos.push_back(all_sdepos[ind]); + continue; + } + + // Prior + const size_t other = iarr(ind, 3); + if (other >= all_sdepos.size()) { + ++npriors_missing; + } + else { + ++npriors; + all_sdepos[other]->set_prior(all_sdepos[ind]); } - // We save both gen=0 and gen>0 as nullptrs to preserve indexing - sdepos.push_back(sdepo); } + if (npriors_missing) { + log->warn("call={}, missing {} prior depos, active={} total={}", + m_count, npriors_missing, sdepos.size(), all_sdepos.size()); + } + all_sdepos.clear(); - WireCell::IDepo::vector depos; + // convert to IDepos + WireCell::IDepo::vector idepos; for (auto sdepo: sdepos) { - if (sdepo) { - auto idepo = IDepo::pointer(sdepo); - sdepo = nullptr; - depos.push_back(idepo); - } + idepos.push_back(sdepo); } + sdepos.clear(); - log->debug("call={} loaded {} depos from ident {} stream {}", - m_count, depos.size(), ident, m_inname); + log->debug("call={} loaded {} active, {} total depos from ident {} stream {} with {} priors", + m_count, idepos.size(), ndepos, ident, m_inname, npriors); - return std::make_shared(ident, depos); + return std::make_shared(ident, idepos); } bool Sio::DepoFileSource::operator()(IDepoSet::pointer& ds) diff --git a/sio/src/FrameFileSource.cxx b/sio/src/FrameFileSource.cxx index ad7e5428d..3a54a707c 100644 --- a/sio/src/FrameFileSource.cxx +++ b/sio/src/FrameFileSource.cxx @@ -52,16 +52,13 @@ void FrameFileSource::configure(const WireCell::Configuration& cfg) m_in.clear(); input_filters(m_in, m_inname); if (m_in.size() < 1) { - THROW(ValueError() << errmsg{"FrameFileSource: unsupported inname: " + m_inname}); + raise("FrameFileSource: unsupported inname: %s", m_inname); } m_tags.clear(); for (auto jtag : cfg["tags"]) { m_tags.push_back(jtag.asString()); } - if (m_tags.empty()) { - m_tags.push_back("*"); - } m_frame_tags.clear(); for (auto jtag : cfg["frame_tags"]) { @@ -69,12 +66,26 @@ void FrameFileSource::configure(const WireCell::Configuration& cfg) } } -bool FrameFileSource::matches(const std::string& tag) +bool FrameFileSource::is_excluded(const std::string& tag) +{ + if (is_tagged(tag)) return false; + if (is_untagged(tag)) return false; + return true; +} + +bool FrameFileSource::is_untagged(const std::string& tag) +{ + return tag.empty() or tag == "*"; +} + +bool FrameFileSource::is_tagged(const std::string& tag) { + if (is_untagged(tag)) { + // avoid misinterpreting "" or "*" in m_tags + return false; + } + for (const auto& maybe : m_tags) { - if (maybe == "*") { - return true; - } if (maybe == tag) { return true; } @@ -106,6 +117,8 @@ bool FrameFileSource::read() m_cur.type = parts[0]; m_cur.tag = parts[1]; m_cur.ext = rparts[1]; + log->debug("read type={} tag={} ident={} ext={} at call={}", + m_cur.type, m_cur.tag, m_cur.ident, m_cur.ext, m_count); return true; } @@ -141,6 +154,8 @@ IFrame::pointer FrameFileSource::load() while (true) { + // Loop over arrays, falling through until we recognize the array type. + // fixme: throw on errors but return nullptr on EOF. // Read next element if current cursor is empty @@ -205,9 +220,8 @@ IFrame::pointer FrameFileSource::load() continue; } - - if (! matches(m_cur.tag)) { - log->trace("call={}, skipping unmatched tag=\"{}\" file={}", + if (is_excluded(m_cur.tag)) { + log->debug("call={}, skipping tag=\"{}\" file={}", m_count, m_cur.tag, m_inname); clear(); continue; @@ -316,7 +330,8 @@ IFrame::pointer FrameFileSource::load() auto itrace = std::make_shared(chid, tbin0, charges); all_traces.push_back(itrace); } - log->trace("call={}, add {} traces to total {}", m_count, nrows, all_traces.size()); + log->trace("call={}, add {} traces from \"{}\" to total {}", + m_count, nrows, framelet.tag, all_traces.size()); } const double time = framelets[0].tickinfo[0]; @@ -329,19 +344,25 @@ IFrame::pointer FrameFileSource::load() // Tag traces of each framelet size_t last_index=0; - for (const auto& tag : m_tags) { - auto& framelet = get_framelet(tag); + for (auto& framelet : framelets) { const size_t size = framelet.channels.size(); + + if (is_untagged(framelet.tag)) { + log->trace("call={}, {} untagged with (\"{}\")", m_count, size, framelet.tag); + last_index += size; + continue; + } + std::vector inds(size); std::iota(inds.begin(), inds.end(), last_index); last_index += size; if (framelet.summary.empty()) { - sframe->tag_traces(tag, inds); + sframe->tag_traces(framelet.tag, inds); } else { - sframe->tag_traces(tag, inds, framelet.summary); + sframe->tag_traces(framelet.tag, inds, framelet.summary); } - log->trace("call={}, tag {} traces with \"{}\"", m_count, size, tag); + log->trace("call={}, tag {} traces with \"{}\"", m_count, size, framelet.tag); } log->trace("call={}, loaded frame with {} tags, {} traces", @@ -363,13 +384,7 @@ bool FrameFileSource::operator()(IFrame::pointer& frame) return false; } - try { - frame = load(); - } - catch (IOError& err) { - log->error("call={}: {}", m_count++, err.what()); - return false; - } + frame = load(); // throws if (frame) { log->debug("call={} load frame: {}", m_count++, Aux::taginfo(frame)); diff --git a/test/docs/talks/spdir/spdir-dag-one-angle.cfg b/test/docs/talks/spdir/spdir-dag-one-angle.cfg new file mode 100644 index 000000000..3cdbad0e9 --- /dev/null +++ b/test/docs/talks/spdir/spdir-dag-one-angle.cfg @@ -0,0 +1,3 @@ +{ "theta_y_deg":[ 90 ], "theta_xz_deg":[ 0 ] } + + diff --git a/test/docs/talks/spdir/spdir-dag.dot b/test/docs/talks/spdir/spdir-dag.dot new file mode 100644 index 000000000..8923ab694 --- /dev/null +++ b/test/docs/talks/spdir/spdir-dag.dot @@ -0,0 +1,18 @@ +digraph spdir { + node[shape=box] + detlinegen drift splat sim sigproc LS4GAN calculator plotter + node[shape=ellipse] + + angles -> detlinegen -> track -> drift -> drifts -> splat -> splats + drifts -> sim -> digits + simFR -> sim + digits -> LS4GAN -> uvcgan + digits -> sigproc + spFR -> sigproc + uvcgan -> sigproc + sigproc->signals + signals -> calculator + splat -> calculator + calculator -> metrics + metrics -> plotter -> plots +} diff --git a/test/docs/talks/spdir/spdir.org b/test/docs/talks/spdir/spdir.org new file mode 100644 index 000000000..c039fe7e9 --- /dev/null +++ b/test/docs/talks/spdir/spdir.org @@ -0,0 +1,523 @@ +#+title: Wire-Cell Signal Processing Directional Metrics (spdir) +#+setupfile: ~/sync/talks/setup.org +#+LATEX_HEADER: \usepackage{siunitx} + + + +* Stuff needed to rebuild this document :noexport: + +** arXiv hates ~wget~ so manually download https://arxiv.org/pdf/1802.08709.pdf + +** Run ~spdir~ so its plots can be included here + +#+begin_src shell +if [ ! -d spdir-plots ] ; then +../../../scripts/spdir --directory spdir-plots all +fi +#+end_src + +#+RESULTS: + +This takes a while and some warnings can be expected: + +#+begin_example +OptimizeWarning: Covariance of the parameters could not be estimated +#+end_example + +Attempt to fit a metric failed. + +#+begin_example +error: (Weights sum to zero, can't be normalized) failed to calculate metrics for pln=0 ch=slice(0, 800, None) splat_filename='pdsp-splats-p0-txz82deg-ty90deg.npz' signal_filename='pdsp-signals-simhi-sphi-p0-txz82deg-ty90deg.npz' +#+end_example + +Attempt to bracket the distribution failed. + +Failures results in 100% bias/res/ineff. + +** Make the DAG visualizations + +#+begin_src shell +../../../scripts/spdir --configfile spdir-dag-one-angle.cfg --dag all > spdir-dag-one-angle.dot +dot -Tpdf -o spdir-dag-one-angle.pdf spdir-dag-one-angle.dot +dot -Tpdf -o spdir-dag.pdf spdir-dag.dot +#+end_src + +#+RESULTS: + + +** Get 2D/q1D FR plots + +#+begin_src shell +wget -q -c https://raw.githubusercontent.com/LS4GAN/toyzero/master/plots/fake-resps-diagnostic.png +wget -q -c https://raw.githubusercontent.com/LS4GAN/toyzero/master/plots/real-resps-diagnostic.png +#+end_src + +#+RESULTS: + +** UVCGAN comes back. + +The UVCGAN data comes back not in WCT frame file format so needs an "import" step done by ~wirecell-ls4gan npz-to-wct~. Besides array reformatting, UVCGAN format drops any ~tickinfo~. + +Checking ~tickinfo~ I sent digits to UVCGAN with $t_0 = -63897.76357827$ and $tbin=0$. Even though I kept working on this while UCVGAN was processing, rerunning still makes digits with this ~tickinfo~. Fix this by adding ~wirecell-util framels~ to spit out a JSON summary of the frame and update ~spdir~ so that it uses this to form a ~tickinfo~ to give to ~wirecell-ls4gan npz-to-wct~. + +** NF is bad + +It seems PDSP NF kills collection channels hit by the tracks. Probably bad-channel list. Remove it and rerun SP on both ~digits~ and ~uvcgan~. To avoid rerunning sim, use + +#+begin_example +$ ../test/scripts/spdir --configfile ../spdir-few.json --config detector=pdsp -p --touch all +$ rm pdsp-mplots* pdsp-signals* pdsp-plots* pdsp-metrics* +$ ../test/scripts/spdir --configfile ../spdir-few.json --config detector=pdsp --config uvcgan_output=uvcgan-output -p all +#+end_example + +Once that looks okay, remove the ~--configfile~ to do all angles ("few" still does 251 jobs, then 1007 more with it removed). + + + + +* ~spdir~ - SP metrics as a function of track direction + +#+begin_quote +~spdir~ provides a *standard calculator* of *efficiency, bias and resolution* of Wire-Cell *signal processing* as a function of *track direction*. +#+end_quote + +- Motivation + +- Calculation overview + +- Results on PDSP + +- Running the process + +- Future work + + +* Motivation for an ~spdir~ process + +A few needs from different areas happened to coincide: + +- \small People started to reimplement MicroBooNE SP-1 paper's "money plot" figure 30 to new detectors. + +- Production of this plot in a "standard" way for multiple detectors and across different WCT software versions was seen as a very strong regression test. + +- The final LS4GAN LDRD paper using LArTPC needs some "physics level" comparisons between nominal simulation and UVCGAN-translated and the ~spdir~ metrics are a very good match for that. + +Though each of these adds its own complexity, they are largely overlapping and so a unified, general approach has been pursued. + + +* ~spdir~ process overview + +** :B_columns: +:PROPERTIES: +:BEAMER_env: columns +:END: + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.15 +:END: + +\includegraphics[width=\textwidth]{spdir-dag.pdf} + + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.8 +:END: + +- Generate *depos* from *ideal line tracks* at *different angles*. +- Apply *drift* and *splat* with proper smearing to get "SP truth denominator". +- Apply *FR* + *sim* (signal+noise) to *drift* to get *digits*. +- Translate *nominal* digits with LS4GAN to get *uvcgan* digits. +- Apply *FR* + *sigproc* both types of digits get *signals*. +- Compare *splat* to *signals* to form *metrics* for given *plane and track angle*. +- Collect *metrics* across a configuration and make the "money plot". + +$\leftarrow$ conceptual DAG, actual DAG for one track angle $\downarrow$ + +\includegraphics[width=0.7\textwidth]{spdir-dag-one-angle.pdf} + + +* Three field responses + +** :B_columns: +:PROPERTIES: +:BEAMER_env: columns +:END: + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.8 +:END: + +\small + +| FR usage: | detector | sigproc | +| | | | +|------------------------+----------+---------| +| real nature | 3D | 2D | +| mock nature | 2D | q1D | +|------------------------+----------+---------| +| unnaturally ideal | 2D | 2D | +| very unnaturally ideal | q1D | q1D | +|------------------------+----------+---------| + + +*q1D* is *intentionally degraded* from 2D: zero all but central wire region. + +- \scriptsize FR variation still exists across the central wire region, thus "q" for "quasi". + +- Consider ~spdir~ metrics using *mock nature* as a *proxy* for *real nature*. + +- Historically we consider *unnaturally ideal* FR pairings. + +- Cover all combos of q1D/2D, $2 \times 2$ multiplier to number of jobs. + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.18 +:END: + + +\includegraphics[width=\textwidth]{real-resps-diagnostic.png} + +\includegraphics[width=\textwidth]{fake-resps-diagnostic.png} + + +* Extra requirements for LS4GAN + +LS4GAN translates sim-q1D $\to$ sim-2D$'$ as a *proxy* for translating sim-2D $\to$ det-3D$'$. + +- Manually export ADC frames from ~spdir~ to UVCGAN environment. +- UVCGAN does *cross translation* (q1D $\leftrightarrow$ 2D). +- Import translated ADC frames back to ~spdir~. +- For each translated FR, SP with either FR applied. + +Results in another $2 \times 2$ multiplier of ~spdir~ variants. + + +* Single track, no LS4GAN ~spdir~ sub-DAG + +\includegraphics[width=\textwidth]{spdir-dag-one-angle.pdf} + +#+begin_center +\tiny just for reference +#+end_center + + +* Extend MicroBooNE SP-1 paper's track direction angle convention + + +** :B_columns: +:PROPERTIES: +:BEAMER_env: columns +:END: + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.6 +:END: + +- \small MB defines angles in the "global" coordinate system. + - Coincident with the W-plane coordinate system. +- MB uses same track for metrics from all 3 planes. + - U/V's wire-track angle different than W's. +- ~spdir~ also defines track angles in each wire plane coordinate system. + - Calculates metrics using both conventions. + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.4 +:END: + +\includegraphics[width=\textwidth,page=17,clip,trim=4.5cm 10.5cm 5cm 11.5cm]{1802.08709.pdf} + +** :B_columns: +:PROPERTIES: +:BEAMER_env: columns +:END: + + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 1.0 +:END: + + +Can configure any angles, but default to MB's choice: +- $\theta_y = 90^\circ$ always (though it is possible to vary this). +- $\theta_{xz} \in \{0, 1, 3, 5, 10, 20, 30, 45, 60, 75, 80, 82, 84, 89\}^\circ$. + + +* Result plots + +Two basic types of sets of plots: + +- Metrics vs track angle in global or wire-plane coordinates. + +- Several detailed metrics plots for a given track angle and wire-plane coordinate system. + +Both have subtypes depending on sim/SP FRs +- "mock reality" (sim-2D/sp-q1D) vs "unnaturally ideal" (sim-2D/sp-2D) + +Current set of plots are all for ProtoDUNE-SP. +- "unnaturally ideal" and "global coordinates" is equivalent to MB "money plot" + +Some examples $\longrightarrow$ + + +* Unnaturally ideal: sim-2D, sp-2D + +** :B_columns: +:PROPERTIES: +:BEAMER_env: columns +:END: + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.5 +:END: + +\includegraphics[width=\textwidth]{spdir-plots/pdsp-plots-simhi-sphi-wirecoords.pdf} + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.5 +:END: + +\includegraphics[width=\textwidth]{spdir-plots/pdsp-plots-simhi-sphi-globalcoords.pdf} + + + + +* Mock reality: sim-2D, sp-q1D + +** :B_columns: +:PROPERTIES: +:BEAMER_env: columns +:END: + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.5 +:END: + +\includegraphics[width=\textwidth]{spdir-plots/pdsp-plots-simhi-splo-wirecoords.pdf} + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.5 +:END: + +\includegraphics[width=\textwidth]{spdir-plots/pdsp-plots-simhi-splo-globalcoords.pdf} + + +* Unnaturally ideal: sim-2D, sp-2D, $3^\circ$ global angles - details + +** :B_columns: +:PROPERTIES: +:BEAMER_env: columns +:END: + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.5 +:END: + +#+begin_export latex +\only<1>{\includegraphics[page=1,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-sphi-p2-txz3deg-ty90deg.pdf}} +\only<2>{\includegraphics[page=2,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-sphi-p2-txz3deg-ty90deg.pdf}} +\only<3>{\includegraphics[page=3,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-sphi-p2-txz3deg-ty90deg.pdf}} +\only<4>{\includegraphics[page=4,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-sphi-p2-txz3deg-ty90deg.pdf}} +\only<5>{\includegraphics[page=5,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-sphi-p2-txz3deg-ty90deg.pdf}} +#+end_export + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.5 +:END: + +\includegraphics[width=\textwidth]{spdir-plots/pdsp-plots-simhi-sphi-globalcoords.pdf} + + + + + +* Unnaturally ideal: sim-2D, sp-2D, $30^\circ$ global angles - details + +** :B_columns: +:PROPERTIES: +:BEAMER_env: columns +:END: + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.5 +:END: + +#+begin_export latex +\only<1>{\includegraphics[page=1,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-sphi-p2-txz30deg-ty90deg.pdf}} +\only<2>{\includegraphics[page=2,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-sphi-p2-txz30deg-ty90deg.pdf}} +\only<3>{\includegraphics[page=3,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-sphi-p2-txz30deg-ty90deg.pdf}} +\only<4>{\includegraphics[page=4,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-sphi-p2-txz30deg-ty90deg.pdf}} +\only<5>{\includegraphics[page=5,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-sphi-p2-txz30deg-ty90deg.pdf}} +#+end_export + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.5 +:END: + +\includegraphics[width=\textwidth]{spdir-plots/pdsp-plots-simhi-sphi-globalcoords.pdf} + + + + + +* Unnaturally ideal: sim-2D, sp-2D, $89^\circ$ global angles - details + +** :B_columns: +:PROPERTIES: +:BEAMER_env: columns +:END: + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.5 +:END: + +#+begin_export latex +\only<1>{\includegraphics[page=1,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-sphi-p2-txz89deg-ty90deg.pdf}} +\only<2>{\includegraphics[page=2,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-sphi-p2-txz89deg-ty90deg.pdf}} +\only<3>{\includegraphics[page=3,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-sphi-p2-txz89deg-ty90deg.pdf}} +\only<4>{\includegraphics[page=4,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-sphi-p2-txz89deg-ty90deg.pdf}} +\only<5>{\includegraphics[page=5,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-sphi-p2-txz89deg-ty90deg.pdf}} +#+end_export + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.5 +:END: + +\includegraphics[width=\textwidth]{spdir-plots/pdsp-plots-simhi-sphi-globalcoords.pdf} + + + + + +* Mock reality: sim-2D, sp-q1D (repeated) + +** :B_columns: +:PROPERTIES: +:BEAMER_env: columns +:END: + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.5 +:END: + +\includegraphics[width=\textwidth]{spdir-plots/pdsp-plots-simhi-splo-wirecoords.pdf} + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.5 +:END: + +\includegraphics[width=\textwidth]{spdir-plots/pdsp-plots-simhi-splo-globalcoords.pdf} + + +* Mock reality: sim-2D, sp-q1D, $10^\circ$ global angles - details + +** :B_columns: +:PROPERTIES: +:BEAMER_env: columns +:END: + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.5 +:END: + +#+begin_export latex +\only<1>{\includegraphics[page=1,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-splo-p2-txz10deg-ty90deg.pdf}} +\only<2>{\includegraphics[page=2,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-splo-p2-txz10deg-ty90deg.pdf}} +\only<3>{\includegraphics[page=3,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-splo-p2-txz10deg-ty90deg.pdf}} +\only<4>{\includegraphics[page=4,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-splo-p2-txz10deg-ty90deg.pdf}} +\only<5>{\includegraphics[page=5,width=\textwidth]{spdir-plots/pdsp-mplots-simhi-splo-p2-txz10deg-ty90deg.pdf}} +#+end_export + +*** :BMCOL: +:PROPERTIES: +:BEAMER_col: 0.5 +:END: + +\includegraphics[width=\textwidth]{spdir-plots/pdsp-plots-simhi-splo-globalcoords.pdf} + + + + + +* Running ~spdir~ + +\small +#+begin_example +cd /path/to/wire-cell-python +pip install -e . +pip install snakemake + +# for convenience to find "spdir" +PATH=/path/to/wire-cell-toolkit/test/scripts:$PATH + +spdir --directory spdir all +ls spdir/ | wc -l +1060 +#+end_example + +* LS4GAN translation + +LS4GAN UVCGAN-translations are optional and require external processing. + +#+begin_example +tar -cvf spdir-pdsp-digits.tar spdir/pdsp-digits-sim*.npz +#+end_example + +Send ~spdir-pdsp-digits.tar~ to Dmitrii for UVCGAN-translation and receive results, unpack and complete the ~spdir~ DAG execution. + +#+begin_example +tar -C spdir -xf translations-from-dmitrii.tar +spdir --directory spdir ls4gan +#+end_example + + +* Programs used by ~spdir~ + +~spdir~ is a [[https://snakemake.readthedocs.io/][Snakemake]] file using: + +- ~wirecell-gen detlinegen~ generate tracks vs angle (Dmitrii). +- ~wirecell-sigproc frzero~ construct degraded q1d FR. +- ~wirecell-test ssss-metrics~ write JSON metric file for each track/variant. +- ~wirecell-test plot-metrics~ plot set of JSON metric files. +- ~wirecell-test plot-ssss~ metric detailed plots +- ~wire-cell -c omnijob.jsonnet~ a generic, do-all job driven by CLI args + - requires detector support in ~cfg/layers/~ scheme (pdsp and uboone so far). + + +* Future work + +** With PDSP + +- Check sign convention with MB SP-1 metrics. +- Switch from Gaussian fit $\mu/\sigma$ to avg/RMS stats to avoid fit failures. +- Splat extra smearing is from a "Morse analysis" assuming only 2D sim/SP FR. + - Should repeat with mixed $2\times 2$ sim/SP FRs. + +** With MicroBooNE + +- Apply to ~uboone~, which has existing ~omnijob.jsonnet~ config support. +- Make quantitative comparison with MB SP-1 plot. + +** Other detectors + +- Develop WCT configuration and run on SBND, PDHD/PDVD, DUNEFD, etc. + - Likely I need help from detector/WCT experts. diff --git a/test/scripts/spdir b/test/scripts/spdir new file mode 100755 index 000000000..afe3bad03 --- /dev/null +++ b/test/scripts/spdir @@ -0,0 +1,477 @@ +#!/usr/bin/env -S snakemake --cores all --snakefile +# -*- snakemake -*- + +# +# This workflow ultimately produces "standard" plots showing per-plane SP +# metrics as a function of track directions for a given detector. +# +# CAVEAT: only canonical detector names (pdsp, uboone) supported by cfg/layers/ can be targeted. +# +# Intermidate file names follow a pattern that allows expressing a variety of combinations of processing stages. +# +# {detector}-{tier}-{modifiers}-p{plane}-txz{txz}deg-ty{ty}deg.{ext}' +# +# - detector :: canonical detector name. +# - variant :: see below +# - plane :: the number of the wire plane (0,1,2) that is targeted. +# - txz :: the theta_xz in targeted wire plane coordinates. +# - ty :: the theta_y in targeted wire plane coordinates. +# +# The "variant" expresses the path-dependent content. It conists of a "tier" and zero or more "tier modifiers". +# +# - depos :: original energy deposits of a the track at the given angles +# - drifts :: depos drifted to the response plane +# - splats :: the drifts "smeared and splatted" to the channels - FIXME: smearing assumes 2D +# - digits :: ADC waveforms, requires "-sim{lo,hi}" modifier indicating if lo (q1d) or hi (2d) FR was used. +# - uvcgan :: ADC waveforms translated by UVCGAN, requires "-tgt{lo,hi}" modifier indicating the target domain. +# - signal :: Signal waveforms from SP requires either a digit or uvcgan modifier and a "-sp{lo,hi}" modifier +# - metrics :: the metrics file for a given variant and angle. requires digits/uvcgan and signal modifiers +# - mplots :: the plot file for a given variant and angle pair. requires digits/uvcgan and signal modifiers +# - plots :: the plot file for a given variant across all angles. requires digits/uvcgan and signal modifiers +# +# Note a *-plot.pdf file holds plots over all angles for a given variant while +# *-plot-*.pdf holds plots for a given variant and angle pair. +# +# Suggested running: +# +# $ spdir --config detector=pdsp --directory=pdsp-spdir +# $ spdir --config detector=uboone --directory=uboone-spdir +# +# This script requires snakemake. Suggest to install inside a venv: +# +# pip install snakemake +# +# Note: as of 2024-02-14 an older version of pulp is required. If you get +# errors mentioning this package try: +# +# pip install 'pulp<2.8' +# +# TODO: +# +# - [ ] the smearing used by splat is tied to the FR used in the "morse" +# running. For now we use smearing determined from 2D regardless what we +# compare to. +# +# - [ ] extend support for uboone. + +# Provide hard-wired configuration defaults. User may over ride them with a +# config file or CLI config parameters. +for k,v in dict( + # The canonical detector name. Future may accept a list. + detector="pdsp", + # Matched arrays for track angles. + # Angle from Y axis to track direction + theta_y_deg=( 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90, 90 ), + # Angle from X axis to projection of track into X-Z plane + theta_xz_deg=( 0, 1, 3, 5, 10, 20, 30, 45, 60, 75, 80, 82, 84, 89 ), + # The plane numbers. Probably universal. + planes = [0,1,2], + # location of uvcgan returns. Set on CLI: --config uvcgan_output=.... + # uvcgan_output = "uvcgan-output". If this is left "None" then UVCGAN + # related subgraphs are not built. + uvcgan_output = None, + +).items(): + config.setdefault(k,v) + +import json + +# Break out some config into simple vars +thetas = list(zip(config['theta_xz_deg'], config['theta_y_deg'])) +txz = config['theta_xz_deg'] +ty = config['theta_y_deg'] +nang = len(ty) +planes = config['planes'] +gens = ['sim'] +if config.get('uvcgan_output', None): + gens=['sim','tgt'] + +def debug(*args, **kwds): + return # MAKE QUIET + import sys + sys.stderr.write(*args, **kwds) + sys.stderr.write('\n') + +class FailsafeDict(dict): + def __getitem__(self, item): + try: + return super().__getitem__(item) + except KeyError: + return "{" + str(item) + "}" +def fname(tier=None, ext='npz', **kwds): + ''' + Generate standard file name. Tier required + ''' + if tier is None: + raise ValueError("tier must be given") + + pre='{detector}-{tier}' + post='p{plane}-txz{txz}deg-ty{ty}deg.{ext}' + simfr='sim{simfr}' # nominal simulation + tgtfr='tgt{simfr}' # uvcgan translation, tgt means "target domain" + genfr='{gen}{simfr}' # wildcard match on "sim" or "tgt" + sp='sp{spfr}' + plotpost='{coords}coords.{ext}' + if tier in ('depos', 'drifts', 'splats'): + parts = [pre, post] + elif tier in ('digits', ): + parts = [pre, simfr, post] + elif tier in ('uvcgan', ): + parts = [pre, tgtfr, post] + elif tier in ('signals','metrics','mplots'): + parts = [pre, genfr, sp, post] + elif tier in ('plots',): + parts = [pre, genfr, sp, plotpost] + else: + raise ValueError(f'unknown tier: {tier}') + + pattern = '-'.join(parts) + + dat = FailsafeDict(kwds, tier=tier, ext=ext) + ret = pattern.format_map(dat) + debug(f'fname({tier=},{ext=},{kwds=})\n\t{ret}') + + return ret + +def expand_thetas(**kwds): + return [fname(txz=txz,ty=ty,**kwds) for txz,ty in thetas] + +def expand_detector(**kwds): + debug(f'expand_detector({kwds=})') + got = expand(expand_thetas(**kwds), plane=planes, detector=[config['detector']], allow_missing=True) + return got + +def expand_nodetector(**kwds): + return expand(expand_thetas(**kwds), plane=planes, detector=['{detector}']) + + +def tracking_params(w): + ''' + Make detector-specific args for detlinegen. + ''' + # MB wires are near X=0 and volume extends to +2548mm. + if w.detector == "uboone": + return "--offset '1*m,0*m,0*m'" + if w.detector == "pdsp": + return "--offset '-1*m,0*m,0*m'" + return "" + +rule tracking: + output: + depos=fname(tier='depos'), + meta=fname(tier='depos', ext='json'), + params: + offset=tracking_params + shell: ''' + wirecell-gen detlinegen \ + {params.offset} \ + --detector={wildcards.detector} \ + --plane={wildcards.plane} \ + --angle-coords=wire-plane \ + --theta_xz="{wildcards.txz}*deg" \ + --theta_y="{wildcards.ty}*deg" \ + --output_depos {output.depos} \ + --output_meta {output.meta} + ''' + +rule all_depos: + input: + expand_detector(tier='depos') + + +# PDSP period is broken. We want to import the FR file anyways, so fix the +# period along the way. For non-PDSP we effectively copy it in as-is. +def import_fr_params(w): + if w.detector == 'pdsp': + return '-P "100*ns"' + return '' + +rule import_fr: + output: + "{detector}-fields-hi.json.bz2" + params: + import_fr_params + shell: ''' + wirecell-resp condition {params} -o {output} {wildcards.detector} + ''' + +rule make_q1dfr: + input: + "{detector}-fields-hi.json.bz2" + output: + "{detector}-fields-lo.json.bz2" + shell: ''' + wirecell-sigproc frzero -n 0 -o {output} {input} + ''' + +rule all_frs: + input: + expand("{detector}-fields-{fr}.json.bz2", detector=[config['detector']], fr=['lo','hi']) + + +# Also makes splat. +# fixme: The extra smearing in splat is tied to the FR which we currently ignore. +rule drift: + input: + cfg=workflow.basedir + "/../../cfg/layers/omnijob.jsonnet", + depos=fname(tier='depos'), + fr='{detector}-fields-hi.json.bz2' + output: + drifts=fname(tier='drifts'), + splats=fname(tier='splats'), + log=fname(tier='drifts',ext='log') + shell: ''' + wire-cell -c {input.cfg} -l {output.log} -L debug \ + -A tasks=drift,splat \ + -A input={input.depos} \ + -A detector={wildcards.detector} \ + -A variant=spdir_hifr \ + --tla-code output='{{drift:"{output.drifts}",splat:"{output.splats}"}}' + ''' + +# A prototypical DAG for the drift job. Bogus file names are used as +# documentation placeholders. +rule drift_dag: + input: + workflow.basedir + "/../../cfg/layers/omnijob.jsonnet", + output: + '{detector}-drift-dag.pdf' + shell: ''' + wirecell-pgraph dotify \ + -A tasks=drift,splat \ + -A input=depos.npz \ + -A detector={wildcards.detector} \ + -A variant=spdir_hifr \ + -A output='{{drift:"drifts.npz",splat:"splats.npz"}}' \ + {input} {output} + ''' + + +rule all_drifts: + input: + expand_detector(tier='drifts') + + +rule simulate: + input: + cfg=workflow.basedir + "/../../cfg/layers/omnijob.jsonnet", + drifts=rules.drift.output.drifts, + fr='{detector}-fields-{simfr}.json.bz2' + output: + digits=fname(tier='digits'), + log=fname(tier='digits',ext='log') + wildcard_constraints: + simfr='lo|hi', + gen='sim' + shell: ''' + wire-cell -c {input.cfg} -l {output.log} -L debug \ + -A tasks=sim \ + -A input={input.drifts} \ + -A detector={wildcards.detector} \ + -A variant=spdir_{wildcards.simfr}fr \ + --tla-code output='{{sim:"{output.digits}"}}' + ''' + +rule simulate_dag: + input: + workflow.basedir + "/../../cfg/layers/omnijob.jsonnet", + output: + '{detector}-simulate-dag.pdf' + shell: ''' + wirecell-pgraph dotify \ + -A tasks=sim \ + -A input=drifts.npz \ + -A detector={wildcards.detector} \ + -A variant=spdir_hifr \ + -A output='digits.npz' \ + {input} {output} + ''' + +rule all_digits: + input: + expand(expand_detector(tier='digits'), simfr=['lo','hi']) + + +rule sigproc: + input: + cfg=workflow.basedir + "/../../cfg/layers/omnijob.jsonnet", + digits=rules.simulate.output.digits, + fr='{detector}-fields-{spfr}.json.bz2' + output: + signals=fname(tier='signals'), + log=fname(tier='signals', ext='log') + wildcard_constraints: + simfr='lo|hi', + gen="sim" + shell: ''' + wire-cell -c {input.cfg} -l {output.log} -L debug \ + -A tasks=sp \ + -A input={input.digits} \ + -A detector={wildcards.detector} \ + -A variant=spdir_{wildcards.spfr}fr \ + --tla-code output='{{sp:"{output.signals}"}}' + ''' + +# We get files back from UVCGAN as: +# +# {uvcgan_output}/tgthi/pdsp-uvcgan-tgthi-p0-txz30deg-ty90deg.npz +# +# and etc for tgtlo and need to convert them to WCT format. Part of this +# conversion involves digging out the t0 of the "digits" file that was exported +# for UVCGAN translation. +# +# CRITICAL: to be hermetic, one must freeze the output directory and all code +# while data is out for translation! +if config.get('uvcgan_output', None): + uvcgan_output=config['uvcgan_output'] + rule import_uvcgan: + input: + digits = fname(tier='digits', gen='tgt'), + uvcgan = uvcgan_output + '/tgt{simfr}/' + fname(tier='uvcgan', gen='tgt') + output: + fname(tier='uvcgan', gen='tgt') + wildcard_constraints: + simfr='lo|hi', + gen="tgt" + run: + # No way to get raw text from shell()? + jtext = '\n'.join(shell('wirecell-util framels {input.digits}', iterable=True)) + print(f'{jtext=}') + jdat = json.loads(jtext) + ti = '{t0},{tick},{tbin}'.format(**jdat) + shell(f"wirecell-ls4gan npz-to-wct --tinfo '{ti}' -n '*' -o {output} {input.uvcgan}") + + rule sigproc_uvcgan: + input: + cfg=workflow.basedir + "/../../cfg/layers/omnijob.jsonnet", + digits=fname(tier='uvcgan', gen='tgt'), + fr='{detector}-fields-{spfr}.json.bz2' + output: + signals=fname(tier='signals', gen='tgt'), + log=fname(tier='signals', gen='tgt', ext='log') + wildcard_constraints: + simfr='lo|hi', + gen="tgt" + shell: ''' + wire-cell -c {input.cfg} -l {output.log} -L debug \ + -A tasks=sp \ + -A input={input.digits} \ + -A detector={wildcards.detector} \ + -A variant=spdir_{wildcards.spfr}fr \ + --tla-code output='{{sp:"{output.signals}"}}' + ''' + + +rule sigproc_dag: + input: + workflow.basedir + "/../../cfg/layers/omnijob.jsonnet", + output: + '{detector}-sigproc-dag.pdf' + shell: ''' + wirecell-pgraph dotify \ + -A tasks=nf,sp \ + -A input=digits.npz \ + -A detector={wildcards.detector} \ + -A variant=spdir_hifr \ + -A output='{{sp:"sigproc.npz"}}' \ + {input} {output} + ''' + + +rule all_signals: + input: + expand(expand_detector(tier='signals'), gen=gens, simfr=['lo','hi'], spfr=['lo','hi']) + + +rule metrics: + input: + splats =fname(tier='splats'), + signals=fname(tier='signals'), + depos =fname(tier='depos', ext='json') + output: + metrics=fname(tier='metrics', ext='json') + shell: ''' + wirecell-test ssss-metrics \ + --output {output.metrics} \ + --params {input.depos} \ + {input.splats} {input.signals} + ''' + +rule all_metrics: + input: + expand(expand_detector(tier='metrics', ext='json'), gen=gens, simfr=['lo','hi'], spfr=['lo','hi']) + + +def plots_args(w): + 'Generate additional args for plots command' + fr = dict(lo='q1D', hi='2D') + simfr = fr[w.simfr] + spfr = fr[w.spfr] + gen = 'sim' + if w.gen == 'tgt': + gen = 'gan' + args = f"--title '{w.detector.upper()} {gen}:{simfr} SP:{spfr}'" + if w.coords == "global": + args += ' --coordinate-plane=2' + return args + +rule plots: + input: + metrics=expand_detector(tier='metrics', ext='json') + output: + metrics=fname(tier='plots', ext='pdf') + params: + args=plots_args + wildcard_constraints: + simfr='lo|hi', + gen='sim|tgt' + shell: ''' + wirecell-test plot-metrics {params.args} -o {output} {input} + ''' + + +rule all_plots: + input: + expand(fname(detector=config['detector'], tier='plots', ext='pdf'), + gen=gens, simfr=['lo','hi'], spfr=['lo','hi'], coords=['wire','global']) +rule all_plots_gan: + input: + expand(fname(detector=config['detector'], tier='plots', ext='pdf'), + gen=['tgt'], simfr=['lo','hi'], spfr=['lo','hi'], coords=['wire','global']) + + +rule metric_plots: + input: + splats =fname(tier='splats'), + signals=fname(tier='signals') + output: + plots=fname(tier='mplots', ext='pdf') + shell: ''' + wirecell-test plot-ssss --output {output.plots} {input.splats} {input.signals} + ''' + + +rule all_mplots: + input: + expand(expand_detector(tier='mplots', ext='pdf'), gen=gens, simfr=['lo','hi'], spfr=['lo','hi']) + +rule all_mplots_gan: + input: + expand(expand_detector(tier='mplots', ext='pdf'), gen=['tgt'], simfr=['lo','hi'], spfr=['lo','hi']) + +rule all: + input: + rules.all_mplots.input, rules.all_plots.input + + +rule all_gan: + input: + rules.all_mplots_gan.input, rules.all_plots_gan.input + +rule all_dags: + input: + expand(rules.drift_dag.output, detector=[config['detector']]), + expand(rules.simulate_dag.output, detector=[config['detector']]), + expand(rules.sigproc_dag.output, detector=[config['detector']]) + +