diff --git a/cfg/layers/gen-mids.sh b/cfg/layers/gen-mids.sh index a8ea38ea1..4c93b5e53 100755 --- a/cfg/layers/gen-mids.sh +++ b/cfg/layers/gen-mids.sh @@ -2,17 +2,26 @@ # in-place generator for mids.jsonnet to save some typing. -cat <mids.jsonnet +me="$(realpath "$BASH_SOURCE")" +mydir="$(dirname "$me")" + +cd "$mydir" + +out="mids.jsonnet" + +cat < "$out" // This file is generated, do not edit. // Generated on $(date) by $USER on $(hostname --fqdn) - +// with $me { EOF +# this loop should keep path relative for one in mids/*/mids.jsonnet do det=$(basename $(dirname $one)) - echo " $det : import \"$one\"," >> mids.jsonnet + echo " $det : import \"$one\"," >> "$out" done -echo "}" >> mids.jsonnet +echo "}" >> "$out" + diff --git a/cfg/layers/midapi.jsonnet b/cfg/layers/midapi.jsonnet index faaf95f85..4ed0a7450 100644 --- a/cfg/layers/midapi.jsonnet +++ b/cfg/layers/midapi.jsonnet @@ -1,5 +1,10 @@ // This file defines the "mid" layer API. // +// It provides a "base class" with dummy methods that will assert if not +// overridden. +// +// Caution: due to laziness it is likely woefully incomplete. +// // 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 diff --git a/cfg/layers/mids.jsonnet b/cfg/layers/mids.jsonnet index c09002864..90ca01a4e 100644 --- a/cfg/layers/mids.jsonnet +++ b/cfg/layers/mids.jsonnet @@ -1,6 +1,6 @@ // This file is generated, do not edit. -// Generated on Tue 03 May 2022 01:14:07 PM EDT by bv on haiku.home - +// Generated on Fri Dec 8 11:53:11 AM EST 2023 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", diff --git a/cfg/layers/mids/pdsp/api/resps.jsonnet b/cfg/layers/mids/pdsp/api/resps.jsonnet new file mode 100644 index 000000000..bb25d2b19 --- /dev/null +++ b/cfg/layers/mids/pdsp/api/resps.jsonnet @@ -0,0 +1,39 @@ +function(params) { + sim: { + er: { + type: "ColdElecResponse", + name: "sim", + data: params.elec + params.ductor.binning, + }, + fr: { + type: "FieldResponse", + name: "sim", + data: { filename: params.ductor.field_file } + }, + rc: { + type: 'RCResponse', + data: params.ductor.binning { width: params.rc.width } + }, + }, + nf: { + fr: { + type: "FieldResponse", + name: "nf", + data: { filename: params.nf.field_file } + }, + }, + sp: { + er: { + type: "ColdElecResponse", + name: "sp", + data: params.elec + params.sp.binning, + }, + fr: { + type: "FieldResponse", + name: "sp", + data: { filename: params.sp.field_file } + }, + }, + +} + diff --git a/cfg/layers/mids/pdsp/api/sim.jsonnet b/cfg/layers/mids/pdsp/api/sim.jsonnet index 2edac550b..2c1ff892f 100644 --- a/cfg/layers/mids/pdsp/api/sim.jsonnet +++ b/cfg/layers/mids/pdsp/api/sim.jsonnet @@ -11,37 +11,27 @@ local wc = low.wc; local pg = low.pg; local idents = low.util.idents; -local frs = import "frs.jsonnet"; +local resps = import "resps.jsonnet"; function(services, params, options={}) { // Signal binning may be extended from nominal. local sig_binning = params.ductor.binning, - local fr = frs(params).sim, - - local cer = { - type: "ColdElecResponse", - data: params.elec + sig_binning, - }, - - local rc = { - type: 'RCResponse', - data: sig_binning { width: params.rc.width } - }, + local res = resps(params).sim, // some have more than one - local short_responses = [ cer, ], - local long_responses = [ rc, ], + local short_responses = [ res.er, ], + local long_responses = [ res.rc, ], local pirs = [{ type: 'PlaneImpactResponse', name: std.toString(plane), - uses: [fr, services.dft] + short_responses + long_responses, + uses: [res.fr, services.dft] + short_responses + long_responses, data: sig_binning { plane: plane, dft: wc.tn(services.dft), - field_response: wc.tn(fr), + field_response: wc.tn(res.fr), 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, @@ -64,15 +54,14 @@ function(services, params, options={}) { params.ductor.start_time) ]), - // API method sim.reframer - reframer :: function(anode) + reframer :: function(anode, tags=[], name=null) pg.pnode({ type: 'Reframer', - name: idents(anode), + name: if std.type(name) == "null" then idents(anode) else name, data: { anode: wc.tn(anode), - tags: [], + tags: tags, fill: 0.0, toffset: 0, tbin: params.ductor.tbin, @@ -82,7 +71,7 @@ function(services, params, options={}) { // API method sim.signal: subgraph making pure signal voltage from // depos. - signal :: function(anode) + signal :: function(anode, tags=[]) pg.pipeline([ pg.pnode({ type: 'DepoTransform', @@ -101,7 +90,7 @@ function(services, params, options={}) { nsigma: 3, }, }, nin=1, nout=1, uses = pirs + [anode, services.random, services.dft]), - $.reframer(anode)]), + $.reframer(anode, tags=tags, name=idents(anode))]), // API method sim.noise: subgraph adding noise to voltage noise :: function(anode) @@ -138,4 +127,45 @@ function(services, params, options={}) { } }, 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 89efa0b4b..71f9773bf 100644 --- a/cfg/layers/mids/pdsp/api/sp.jsonnet +++ b/cfg/layers/mids/pdsp/api/sp.jsonnet @@ -7,7 +7,7 @@ local spfilt = import "sp-filters.jsonnet"; local low = import "../../../low.jsonnet"; local wc = low.wc; -local frs = import "frs.jsonnet"; +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. @@ -19,12 +19,7 @@ function(services, params, options={}) function(anode) local fullscale = pars.digi.fullscale[1] - pars.digi.fullscale[0]; local ADC_mV_ratio = ((1 << resolution) - 1 ) / fullscale; - local fr = frs(pars).sp; - - local cer = pars.ductor.binning { - type: "ColdElecResponse", - data: pars.elec, - }; + local res = resps(params).sp; low.pg.pnode({ type: 'OmnibusSigProc', @@ -38,8 +33,8 @@ function(services, params, options={}) function(anode) anode: wc.tn(anode), dft: wc.tn(services.dft), per_chan_resp: "", - field_response: wc.tn(fr), - elecresponse: wc.tn(cer), + field_response: wc.tn(res.fr), + elecresponse: wc.tn(res.er), ftoffset: 0.0, // default 0.0 ctoffset: 1.0*wc.microsecond, // default -8.0 fft_flag: 0, // 1 is faster but higher memory, 0 is slightly slower but lower memory @@ -82,4 +77,4 @@ function(services, params, options={}) function(anode) isWrapped: false, sparse : opts.sparse, }, - }, nin=1, nout=1, uses=[anode, services.dft, fr, cer] + spfilt) + }, nin=1, nout=1, uses=[anode, services.dft, res.fr, res.er] + spfilt) diff --git a/cfg/layers/mids/pdsp/mids.jsonnet b/cfg/layers/mids/pdsp/mids.jsonnet index c4f50f578..ce660be0d 100644 --- a/cfg/layers/mids/pdsp/mids.jsonnet +++ b/cfg/layers/mids/pdsp/mids.jsonnet @@ -3,10 +3,11 @@ // 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"; function(services, variant="", options={}) api(services, variants[variant], options) diff --git a/cfg/layers/mids/pdsp/variants/nominal.jsonnet b/cfg/layers/mids/pdsp/variants/nominal.jsonnet index e4b5d0969..67cddac98 100644 --- a/cfg/layers/mids/pdsp/variants/nominal.jsonnet +++ b/cfg/layers/mids/pdsp/variants/nominal.jsonnet @@ -1,13 +1,28 @@ -// This is the base for *pdsp* variant parameter objects. - -// Do NOT place any "if" in this file! If a new variant is needed, -// make an empty, new file in this directory and either import/inherit -// from this structure overiding only the new values. Do not -// copy-paste. Refactor this file if needed to achieve breif -// variants. +// This is the file cfg/layers/mids/pdsp/variant/nominal.jsonnet. +// If this file is found at any other path, I will delete it without notice. + +// This is the BASE configuration for *pdsp* variant parameter objects. It +// yields an object of a schema expected throughout cfg/layers/mids/pdsp/. + +// Do NOT change this file. + +// If the nominal parameter do not suit your needs, follow these steps: +// +// 1. Pick a single, short name that describes your needs. +// +// 2. Start an EMPTY file in cfg/layers/mids/pdsp/variant/.jsonnet +// +// 3. Import the nominal.jsonnet object +// +// 4. Yield the nominal object with any values overridden. +// +// 5. Add this new object to the variants object in mids.jsonnet. local wc = import "wirecell.jsonnet"; +local detectors = import "detectors.jsonnet"; +local mydet = detectors.pdsp; + { // Define the LAr properties. In principle, this is NOT subject // to variance though in principle it could be (eg, top/bottom of @@ -54,7 +69,8 @@ local wc = import "wirecell.jsonnet"; // The exhaustive list of the location of every single wire // (or strip) segment. - wires_file: "protodune-wires-larsoft-v4.json.bz2", + // 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, @@ -162,7 +178,8 @@ local wc = import "wirecell.jsonnet"; // The ductor transforms drifted depos to currents ductor: { - field_file: "dune-garfield-1d565.json.bz2", + # 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 @@ -198,10 +215,23 @@ local wc = import "wirecell.jsonnet"; 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: "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 @@ -274,7 +304,8 @@ local wc = import "wirecell.jsonnet"; // Imaging paramter pack img : { // For now we use MicroBooNE's - "charge_error_file": "microboone-charge-error.json.bz2", + #"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, diff --git a/cfg/layers/mids/pdsp/variants/simple.jsonnet b/cfg/layers/mids/pdsp/variants/simple.jsonnet new file mode 100644 index 000000000..8ca144938 --- /dev/null +++ b/cfg/layers/mids/pdsp/variants/simple.jsonnet @@ -0,0 +1,321 @@ +// 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 new file mode 100644 index 000000000..8d6cef948 --- /dev/null +++ b/cfg/layers/mids/pdsp/variants/ssss.jsonnet @@ -0,0 +1,48 @@ +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 + { + splat: super.splat + { + + // Run wirecell-gen morse-* to find these numbers that match the extra + // spread the sigproc induces. + "smear_long": [ + 2.691862363980221, + 2.6750200122535057, + 2.7137567141154055 + ], + "smear_tran": [ + 0.7377218875719689, + 0.7157764520393882, + 0.13980698710556544 + ] + } +}; + + +{ + // Used for test-morse-pdsp + morse_nominal: override, + + // Used for test-ssss-pdsp + ssss_nominal: override, + ssss_smeared: smeared, + + // Used for test-spdir + spdir: smeared +} diff --git a/cfg/layers/mids/uboone/api/sim.jsonnet b/cfg/layers/mids/uboone/api/sim.jsonnet index 2d761b9cf..11db507e0 100644 --- a/cfg/layers/mids/uboone/api/sim.jsonnet +++ b/cfg/layers/mids/uboone/api/sim.jsonnet @@ -98,6 +98,13 @@ function(services, params, options={}) { local pirs = make_pirs(index); pg.pipeline([ pg.pnode({ + // despite the name, just calls whatever its "drifter" is + type: 'DepoSetDrifter', + name: broken_detector[index].name, + data: { + drifter: 'WireBoundedDepos:' + broken_detector[index].name, + } + }, nin=1, nout=1, uses=[pg.pnode({ type: 'WireBoundedDepos', name: broken_detector[index].name, data: { @@ -105,15 +112,7 @@ function(services, params, options={}) { regions: broken_detector[index].regions, mode: broken_detector[index].mode, }, - }, nin=1, nout=1, uses = [anode]), - pg.pnode({ - type: 'Bagger', - name: broken_detector[index].name, - data: { - gate: [params.ductor.start_time, - params.ductor.start_time+params.ductor.readout_time], - }, - }, nin=1, nout=1), + }, nin=1, nout=1, uses = [anode])]), transforms(anode, index, pirs)[params.ductor.transform], ]), @@ -145,10 +144,9 @@ function(services, params, options={}) { local ubsigtags = ['ubsig%d'%n for n in [0,1,2]], // API method sim.signal: subgraph making pure signal voltage from - // depos. Note, you MUST provide a stream of individial depos, - // not a stream of deposets! + // depos. signal :: function(anode) pg.pipeline([ - pg.fan.pipe('DepoFanout', [make_branch(anode, index) for index in [0,1,2]], + pg.fan.pipe('DepoSetFanout', [make_branch(anode, index) for index in [0,1,2]], 'FrameFanin', 'ubsigraph', ubsigtags), pg.pnode({ type: 'Reframer', diff --git a/cfg/wirecell.jsonnet b/cfg/wirecell.jsonnet index 7adf97b08..bc0ed59c5 100644 --- a/cfg/wirecell.jsonnet +++ b/cfg/wirecell.jsonnet @@ -330,6 +330,16 @@ unique_list(l):: std.foldl($.unique_helper, l, []), + // Return an array. If l is array, return it. If string, split it, if object return field names + listify(l, d=',') :: + local t = std.type(l); + if t == "string" then + std.split(l, d) + else if t == "object" then + std.objectValues(l) + else + l, + // Round a floating point to nearest integer. It's a bit weird to // go through a format/parse. Maybe there's a better way? roundToInt(x):: std.parseInt("%d" % (x+0.5)), diff --git a/gen/docs/depofluxsplat-summary.org b/gen/docs/depofluxsplat-summary.org new file mode 100644 index 000000000..ed598982c --- /dev/null +++ b/gen/docs/depofluxsplat-summary.org @@ -0,0 +1,19 @@ + +#+caption: Sensitive volume diagonal and ideal line track endpoints. +#+name: tab:diagonal-endpoints +| x (mm) | y (mm) | z (mm) | +|--------|--------|--------| +| -3578.36| 76.1| 3.3497| +| -1.5875| 6060.0| 2303.37| + + +#+caption: Depo t/x ranges: +#+name: tab:depo-tx-ranges +| | min | max | units | +|---------|------|-----|-------| +| depos t | 0.00|0.024 | us | +| drift t | -57.50|2227.381 | us | +| depos x | -3578.36|-1.59 | mm | +| drift x | -3487.89|-3487.89 | mm | + + diff --git a/gen/docs/depofluxsplat.org b/gen/docs/depofluxsplat.org new file mode 100644 index 000000000..deea0479c --- /dev/null +++ b/gen/docs/depofluxsplat.org @@ -0,0 +1,128 @@ +#+title: DepoFluxSplat +#+subtitle: Converting Drifted Depos to Signal-like waveforms +#+LATEX_HEADER: \usepackage[margin=1.0in]{geometry} +#+options: ':t +#+BIND: org-latex-prefer-user-labels t + +* Overview + +The ~DepoFluxSplat~ component transforms depos from an input ~IDepoSet~ to an output ~IFrame~. The transform involves adding in quadrature additional Gaussian extent to any extent the depos may posses (eg, as gained during drifting) and essentially binning the result onto a grid formed by the wires of each plane and time sampling. The resulting waveforms can be considered "true signals" and may be meaningfully compared to the "reconstructed signals" that result in the same ~IDepoSet~ passing through the simulation (~DepoTransform~) and signal processing (~OmnibusSigProc~). + +See also ~DepoFluxWriter~ component that resides in LArSoft's ~larwirecell~. It performs a similar transform but outputs to the ~art::Event~ a collection of ~SimChannel::IDE~ objects. ~DepoFluxWriter~ operates on multiple ~AnodePlane~ together while ~DepoFluxSplat~ operates on a per-anode context. In WCT there is also an older ~DepoSplat~ that operates on a stream of individual ~IDepo~. This component is deprecated by ~DepoFluxSplat~. + +* Interpreting the output + +The ~DepoFluxSplat~ can produce a *sparse* or a *dense* frame. A sparse frame will consist of a set of distinct traces for each depo. These traces contribute a portion of the true signal due to their associated depo. This means multiple traces may be found in the frame for any given channel ident. A dense frame will accumulate these contributes so that no more than one trace will be present for any given channel ident. + +Whether sparse or dense, the resulting true signal waveforms represent a measure of the "true" flux of ionization electrons through a particular *measurement plane* (MP) that is perpendicular to the nominal drift direction. The MP is identified with the *collection plane* (CP) by default. Options are described below to define another location for the MP. + +Another logical plane that is of importance here is the *response plane* (RP). It is at the RP where the drift simulation (~Drifter~) piles up depos prior to these depos being output and passed to the detector response simulation (~DepoTransform~) in the full simulation. These same depos are typically input to ~DepoTransform~ which will further drift them from the RP to the MP assuming a constant electric field. + +At the MP, the distribution of depo charge is binned in drift time and spatially in the MP. In time, the flux is integrated over small sample periods ("ticks"). In space, the flux is integrated over the region near all the electrodes that are connected to a common electronics channel. A region is taken to run the length of an electrode (wire or strip) and extend on either side of its center line by one half pitch. + +* Nominal time + +The *nominal time* is defined such that it measures when the depo charge distribution passes through the MP. It is the sum of the following times: + +- absolute kinematics time :: of the initial kinematics objects (eg "Geant4 time"). +- interaction time :: to create a particular charged particle +- tracking time :: to reach a given step of that particle track that produces a depo. +- drifting time :: to drift that depo to the RP. +- transport time :: to transport that depo to the MP. + +Thus, the *nominal time* measures when a portion of a depo passes through the *measurement plane*. This time is absolute, physical and correct given the physics models. The nominal time is independent from any man-made times. In particular it is unrelated to any "trigger time". + +* Acceptance window + +An *acceptance window* governs what depos will be considered and is reflected in the *frame time*. It is defined in the *nominal time* clock. As such, it must be made to bracket the expected arrival times to the MP for all depos that the user wishes to include. In defining the *acceptance window* the user must carefully consider how *nominal time* is defined above + +For example, if Geant4 is given an "event start time" of T0 and assuming interaction and tracking times are neglected, one may pick a window starting from T0 with a duration of one maximum drift time. This will exactly capture all depos throughout the drift volume. + +The acceptance window is configured using: + +- ~tick~ :: The sample time over which to integrate depo flux into time bins. +- ~window_start~ :: The start of the acceptance window in *nominal time* measured at the MR. +- ~window_duration~ :: The duration of the acceptance window. + +The actual acceptance window will truncate ~window_duration~ so that it is an integral number of ~tick~. By default, the ~window_start~ will be used to set the frame reference time. + +* Incorrecting time offsets + +To provide users with questionable flexibility, ~DepoFluxSplat~ accepts two types of arbitrary time offsets. These are applied by arbitrarily adjusting either the frame reference time or the ~tbin~ times of traces. These time offsets are not considered when applying the *acceptance window*. The configuration parameters are: + +- ~reference_time~ :: an absolute time *subtracted* from the *nominal time* in setting the frame time. Default is 0. + +- ~time_offsets~ :: a 3-array providing a relative time *added* to the *nominal time* on a per-plane basis. Default is empty. If defined, the time offsets are converted to integer number of ticks and added to the ~tbin~ value of all traces on a per-plane basis. + +* Additional smearing + +~DepoFluxSplat~ can apply additional smearing to depos independently in the longitudinal (time) and the transverse direction of each plane with the following options: + +- ~smear_long~ :: a unitless, scalar value giving the number of ticks over which to smear. It forms a second Gaussian sigma that is added in quadrature with the original longitudinal sigma of the depo. + +- ~smear_tran~ :: is a unitless, scalar or 3-array giving number of pitches over which to smear. If scalar, the same smearing is applied regardless of the plane. + +- ~nsigma~ :: the cut-off of the Gaussian extent of a depo in units of number-of-sigma. + +* Anode and field response + +- ~anode~ :: names an ~IAnodePlane~ giving the context for the instance of ~DepoFluxSplat~. + +- ~field_response~ :: names an ~IFieldResponse~ giving the field response data to assume. Only the ~origin~ and ~speed~ values are used to calculate the *transport time* defined above. + +* Example + +The test ~test/test/test-ssss-pdsp.bats~ uses the ~pdsp~ "simple" detector variant (not the "nominal" which meant to be close to reality) to perform splat, sim and sigproc (ssss). + +The sensitive volume bounding box for anode 0 face 0 is described in table [[tab:diagonal-endpoints]]. +The two X coordinated are on the "anode" and "cathode" planes. +Consider an ideal line track of ionization charge connecting these two corners. +The range in the time and X coordinates for the original depos dotting this line and after drifting are showin in table [[tab:depo-tx-ranges]]. + +#+include: depofluxsplat-summary.org + +Before drifting, depos are "piled up" in time due to the particle moving at the speed of light and their X coordinates span the full drift distance as the track is diagonal through the entire drift volume. After drifting, the depos are "piled up" on the response plane and their times are spread due to differing drift lengths. The negative times are due to original depos between the response and anode planes being "backed up" to the respone plane. + + +Figure [[fig:comp]] visualizes the ionized electron charge distribution in the three views as determined by the ~DepoFluxSplat~ ("splat") and full simulation and signal processing reconstruction ("signal"). Figure [[fig:zoom]] shows the regions at the beginning and ending of the track for the V-plane view for splat and signal and their difference. + +#+attr_latex: +#+caption: Comparison of the diagonal line track for ~DepoFluxSplat~ (top, "splat") and full simulation and signal processing (bottom, "signal"). The images show channel vs time and the inset plots show activity integrated along each axis. +#+name: fig:comp +[[file:depofluxsplat/plots-nosmear0000.png]] + +#+caption: The splat and signal activity near the two ends of the track and their pixel-wise difference. +#+name: fig:zoom +#+begin_center +#+attr_latex: :width 0.4\textwidth :center +[[file:depofluxsplat/plots-nosmear0003.png]] +#+attr_latex: :width 0.4\textwidth :center +[[file:depofluxsplat/plots-nosmear0004.png]] +#+end_center + +The substantial difference indicates a mismatch between the two. +Figure [[fig:wf]] selects particular channels from these two regions and show their time domain waveforms. + +#+caption: The splat and signal activity near the end points of the track for select channels. +#+name: fig:wf +#+begin_center +#+attr_latex: :width 0.4\textwidth :center +[[file:depofluxsplat/plots-nosmear0005.png]] +#+attr_latex: :width 0.4\textwidth :center +[[file:depofluxsplat/plots-nosmear0006.png]] +#+end_center + + +These clearly show differences. Some difference is expected. The splat shape comes from integrating regions of purely Gaussian-distributed charge from the depos while simulation applies a convolution and sum over the field response which a very non-Gaussian kernel. Furthermore, the signal processing applies a low-pass filter which adds width to the resulting signal. The fact that their integrals do not match is not understood. + +The ~DepoFluxSplat~ allows two arbitrary Gaussian widths (transverse and longitudinal) to be added in quadrature to those carried by the drifted depos. The figure [[fig:wfss]] illustrate an approximation of what these additional 2D Gaussians can do to the waveforms by directly convolving a Gaussian along only the time dimension (not channel). The "smear" value is the Gaussian sigma measured in units of ticks. The "scale" value is multiplicative. + + +#+caption: The splat and signal activity near the end points of the track for select channels with smearing and scaling applied. +#+name: fig:wfss +#+begin_center +#+attr_latex: :width 0.4\textwidth :center +[[file:depofluxsplat/plots-smear0005.png]] +#+attr_latex: :width 0.4\textwidth :center +[[file:depofluxsplat/plots-smear0006.png]] +#+end_center diff --git a/gen/inc/WireCellGen/DepoFluxSplat.h b/gen/inc/WireCellGen/DepoFluxSplat.h new file mode 100644 index 000000000..ee37b4de5 --- /dev/null +++ b/gen/inc/WireCellGen/DepoFluxSplat.h @@ -0,0 +1,132 @@ +/** DepoFluxSplat is an approximate combination DepoTransform and OmnibusSigProc. + + DepoFluxSplat transforms an IDepoSet into an IFrame that represents a "true + signal" comparable to the output of signal processing. It does this by + adding configurable amount of Gaussian smearing to the drifted depos to + account for the smearing accrued by the convolution/deconvolution cycle of + sim+sigproc and then bins the resulting distribution onto channels via their + wires. + + DepoFluxSplat operates on the context of a single anode and will write + either a sparse or dense frame. When the frame is sparse, each depo results + in a family of traces that will in general overlap in tick and channel with + the family of traces from nearby depos. + + Differences between DepoFluxSplat and DepoSplat: + - DepoFluxSplat is general to all detectors (no hard-wired magic numbers) + - DepoFluxSplat consumes IDepoSet instead of IDepo + + Differences between DepoFluxSplat and DepoFluxWriter (from larwirecell) + - DepoFluxSplat is "pure" Wire-Cell, no interace with art::Event. + - DepoFluxSplat outputs IFrame instead of SimChannel + - DepoFluxSplat operates in the context of a single AnodePlane. + - DepoFluxSplat drops the misguided handling of "energy". + + DepoFluxSplat internals are largely copy-pasted from DepoFluxWriter. It is + hoped that DepoFluxWriter will be refactored to replace its internals with + instances of DepoFluxSplat. + +*/ +#ifndef WIRECELLGEN_DEPOFLUXSPLAT +#define WIRECELLGEN_DEPOFLUXSPLAT + +#include "WireCellAux/Logger.h" +#include "WireCellIface/IDepoFramer.h" +#include "WireCellIface/IConfigurable.h" +#include "WireCellIface/IAnodePlane.h" +#include "WireCellIface/IAnodeFace.h" + +namespace WireCell::Gen { + + class DepoFluxSplat : public Aux::Logger, + public IDepoFramer, + public IConfigurable + { + public: + DepoFluxSplat(); + ~DepoFluxSplat(); + + virtual bool operator()(const input_pointer& in, output_pointer& out); + + virtual void configure(const WireCell::Configuration& cfg); + virtual Configuration default_configuration() const; + + private: + + /// Configuration: + /// + /// anode - the type/name for the IAnodePlane. + /// + IAnodePlane::pointer m_anode; + + /// field_response - name of an IFieldResponse from which drift_speed + /// (speed) and response_plane (origin) is taken. These two parameters + /// may be explicitly overridden (see next). If both are overridden + /// then field_response parameter may be left unspecified. + + /// drift_speed, response_plane - give explicit drift speed and + /// response plane distance. Default is unspecified and the values + /// provided by the field response are used. + double m_speed{0}, m_origin{0}; + + /// tick - the sample time over which to integrate depo flux + /// into time bins. + /// + /// window_start - the start of an acceptance window in NOMINAL + /// TIME. + /// + /// window_duration - the duration of the acceptance window. + Binning m_tbins; + + /// sparse - bool, if true save a sparse IFrame, else dense. Default is + /// false. Sparse frames can be very large an only useful for special + /// studies that require the many-to-many mapping between depos and + /// trace contributions. + bool m_sparse{false}; + + /// nsigma - number of sigma at which to truncate a depo Gaussian + double m_nsigma{3.0}; + + /// referenced_time - An arbitrary time that is SUBTRACTED from the + /// NOMINAL TIME of the window start when setting the output frame time. + double m_reftime{0}; + + /// time_offsets - An arbitrary time that is ADDED to NOMINAL TIMES for + /// each plane. If set, this offset is reflected in the "tbin" value of + /// traces for a given plane. + std::vector m_tick_offsets; + + /// smear_tran - Extra smearing applied to the depo in the + /// transverse direction of each plane. This given in units of + /// pitch (ie, smear_tran=2.0 would additionally smear over 2 + /// pitch distances) and it is added in quadrature with the + /// Gaussian sigma of the depo. If zero (default) no extra + /// smearing is applied. + /// + /// A single scalar number may be provided which will add a common + /// smearing for all planes. Or, a per-plane smearing may be provided. + std::vector m_smear_tran; + + /// smear_long - Extra smearing applied to the depo in the + /// longitudinal direction. This given in units of tick (ie, + /// smear_long=2.0 would additionally smear across two ticks) + /// and added in quadrature with the Gaussian sigma of the + /// depo. If zero (default) no extra smearing is applied. + /// + /// A single scalar number may be provided which will add a common + /// smearing for all planes. Or, a per-plane smearing may be provided. + std::vector m_smear_long; + + + // internal, return nullptr if depo is not in anode plane. o.w. return + // anode face. + IAnodeFace::pointer find_face(const WireCell::IDepo::pointer& depo); + + // Count executions for debug log. + size_t m_count{0}; + + }; + +} + +#endif diff --git a/gen/inc/WireCellGen/DepoSplat.h b/gen/inc/WireCellGen/DepoSplat.h index 730d9e3d2..0ee016306 100644 --- a/gen/inc/WireCellGen/DepoSplat.h +++ b/gen/inc/WireCellGen/DepoSplat.h @@ -3,8 +3,7 @@ * It is approximately equivalent to combined simulation and sigproc * using the same response. * - * FIXME: A new DepoSetSplat needs to be written as an IDepoFramer to - * avoid the high cost of sending individual depos to N APAs. + * See DepoSetSplat for a version that input depo sets. */ #ifndef WIRECELLGEN_DEPOSPLAT diff --git a/gen/inc/WireCellGen/DuctorFramer.h b/gen/inc/WireCellGen/DuctorFramer.h new file mode 100644 index 000000000..8adcf625c --- /dev/null +++ b/gen/inc/WireCellGen/DuctorFramer.h @@ -0,0 +1,51 @@ +/** + Apply a ductor across a depo set to produce a frame. + + This delegates to two sub DFP nodes: + + - an IDuctor is fed a stream of IDepos from the input IDepoSet and produces + zero or more IFrames. + + - if the set of IFrame produced over the IDepoSet has more than one entry, + the set is fed an IFrameFanin in order to merge them into a single IFrame. + Note the IFrameFanin must be able to handle DYNAMIC MULTIPLICITY. In the + most likely case of using the concrete FrameFanin class for this, the + multiplicity should be set to 0. + */ + +#ifndef WIRECELLGEN_DUCTORFRAMER +#define WIRECELLGEN_DUCTORFRAMER + +#include "WireCellIface/IDepoFramer.h" +#include "WireCellIface/IConfigurable.h" +#include "WireCellIface/IFrameFanin.h" +#include "WireCellIface/IDuctor.h" +#include "WireCellAux/Logger.h" + +namespace WireCell::Gen { + + class DuctorFramer : public Aux::Logger, + public IDepoFramer, + public IConfigurable + { + public: + + DuctorFramer(); + virtual ~DuctorFramer(); + + virtual bool operator()(const input_pointer& in, output_pointer& out); + + virtual void configure(const WireCell::Configuration& cfg); + virtual WireCell::Configuration default_configuration() const; + + private: + + IDuctor::pointer m_ductor{nullptr}; + IFrameFanin::pointer m_fanin{nullptr}; + size_t m_count{0}; + }; +} + +#endif + + diff --git a/gen/inc/WireCellGen/FrameFanin.h b/gen/inc/WireCellGen/FrameFanin.h index b8fcaa01f..a05eb0aed 100644 --- a/gen/inc/WireCellGen/FrameFanin.h +++ b/gen/inc/WireCellGen/FrameFanin.h @@ -1,3 +1,23 @@ +/** Fanin N frames and output one frame. + + The FrameFanin has two operational modes. + + Fixed multiplicity mode is set by a non-zero "multiplicity" value. This + mode is appropriate for FrameFanin to operate in a DFP graph. + + Dynamic multiplicity mode is set by a zero "multiplicity" value. This mode + is relevant for using FrameFanin as a tool inside another DFP node which may + generate a variable set of frames that need to be merged into one. In this + mode it will accept any size input vector of frames. Arrays of trace tags + and tag rules may still be supplied. The arrays are mapped to the vector of + input frames similarly to the case of fixed multiplicity. However, any + frames at indices in the vector equal to or larger than the size of an array + will use the last element in that array. + + See also comments around the configuration parameters. + + */ + #ifndef WIRECELL_GEN_FRAMEFANIN #define WIRECELL_GEN_FRAMEFANIN @@ -30,7 +50,10 @@ namespace WireCell { virtual WireCell::Configuration default_configuration() const; private: - size_t m_multiplicity; + + std::string tag(size_t port) const; + + size_t m_multiplicity{0}; std::vector m_tags; tagrules::Context m_ft; int m_count{0}; diff --git a/gen/inc/WireCellGen/Reframer.h b/gen/inc/WireCellGen/Reframer.h index 2b0a0bf44..965287aa9 100644 --- a/gen/inc/WireCellGen/Reframer.h +++ b/gen/inc/WireCellGen/Reframer.h @@ -58,14 +58,54 @@ namespace WireCell { ITrace::vector process_one(const ITrace::vector& itraces); /// no summary version IAnodePlane::pointer m_anode; - // Consider traces with these tags. No tags mean all traces. + /// Configure: tags + /// + /// Will reframe each set of tagged traces separately, carrying + /// forward the tag to the output. + /// + /// If empty, then all traces in the frame are reframed regardless + /// of any trace tagging. std::vector m_input_tags; - // If not empty, apply this tag to output frame - std::string m_frame_tag; - - double m_toffset, m_fill; - int m_tbin, m_nticks; + /// Configure: ignore_tags + /// + /// If input_tags is empty and yet tagged traces are found in the + /// frame it signifies the user likely made an error and a warning + /// is printed. Set this to true to supress the warning. A true + /// value is inconsistent with a non-empty "tags" array. + bool m_ignore_tags{false}; + + /// Configure: frame_tag + /// + /// If set, apply as a frame tag on the output frame. + std::string m_frame_tag{""}; + + /// Configure: toffset + /// + /// The output frame reference time is nominally set to the the + /// input frame reference time and extended by "tbin" worth of + /// "tick" in the output traces. The "toffset" time will be + /// arbitrarily added to that nominal time. + double m_toffset{0}; + + /// Configure: fill + /// + /// The value to assign to output trace samples that are not + /// represented in the input frame. + double m_fill{0}; + + /// Configure: tbin + /// + /// The starting time bin (tick) of the traces. This effectively + /// truncates the output frame array. + int m_tbin{0}; + + /// Configure: nticks + /// + /// The size of the output frame array along the time dimension. + int m_nticks{0}; + + // count calls for more useful log messages. size_t m_count{0}; }; } // namespace Gen diff --git a/gen/src/DepoFluxSplat.cxx b/gen/src/DepoFluxSplat.cxx new file mode 100644 index 000000000..837f1c5bd --- /dev/null +++ b/gen/src/DepoFluxSplat.cxx @@ -0,0 +1,407 @@ +#include "WireCellGen/DepoFluxSplat.h" + +#include "WireCellIface/IFieldResponse.h" + +#include "WireCellAux/SimpleFrame.h" +#include "WireCellAux/SimpleTrace.h" + +#include "WireCellUtil/Range.h" +#include "WireCellUtil/Array.h" +#include "WireCellUtil/Configuration.h" +#include "WireCellUtil/NamedFactory.h" +#include "WireCellUtil/Units.h" + +#include "WireCellGen/GaussianDiffusion.h" + +#include // std::plus + +WIRECELL_FACTORY(DepoFluxSplat, + WireCell::Gen::DepoFluxSplat, + WireCell::INamed, + WireCell::IConfigurable, + WireCell::IDepoFramer) + +using namespace WireCell; +using namespace WireCell::Aux; +using WireCell::Range::irange; + + +Gen::DepoFluxSplat::DepoFluxSplat() + : Aux::Logger("DepoFluxSplat", "gen") +{ +} + + +Gen::DepoFluxSplat::~DepoFluxSplat() +{ +} + + +WireCell::Configuration Gen::DepoFluxSplat::default_configuration() const +{ + Configuration cfg; + + // Accept array of strings or single string + cfg["anode"] = Json::arrayValue; + cfg["field_response"] = Json::stringValue; + + // time binning + const double default_tick = 0.5 * units::us; + cfg["tick"] = default_tick; + cfg["window_start"] = 0; + cfg["window_duration"] = 8096 * default_tick; + + cfg["sparse"] = m_sparse; + + cfg["nsigma"] = m_nsigma; + + cfg["reference_time"] = 0; + cfg["time_offsets"] = Json::arrayValue; + + cfg["smear_long"] = 0.0; + cfg["smear_tran"] = 0.0; + + return cfg; +} + +static +std::vector get_n(const Configuration& cfg, size_t n=3) +{ + if (cfg.isDouble()) { + return std::vector(n, cfg.asDouble()); + } + if (cfg.size() == n) { + std::vector ret; + for (auto const& one : cfg) { + ret.push_back(one.asDouble()); + } + return ret; + } + return std::vector(n, 0); +} + +void Gen::DepoFluxSplat::configure(const WireCell::Configuration& cfg) +{ + // For response plane info + auto jfrname = cfg["field_response"]; + if (! jfrname.isNull()) { + auto ifr = Factory::find_tn(jfrname.asString()); + const auto& fr = ifr->field_response(); + m_speed = fr.speed; + m_origin = fr.origin; + } + auto jspeed = cfg["drift_speed"]; + if (! jspeed.isNull()) { + m_speed = jspeed.asDouble(); + } + auto jorigin = cfg["response_plane"]; + if (! jorigin.isNull()) { + m_origin = jorigin.asDouble(); + } + + // Anode plane for down-selecting depos. + std::string anode_tn = cfg["anode"].asString(); + m_anode = Factory::find_tn(anode_tn); + + // Acceptance window + const double wtick = get(cfg, "tick", 0.5 * units::us); + const double wstart = cfg["window_start"].asDouble(); + const double wduration = cfg["window_duration"].asDouble(); + const int nwbins = wduration / wtick; + m_tbins = Binning(nwbins, wstart, wstart + nwbins * wtick); + + // Frame form. + m_sparse = get(cfg, "sparse", m_sparse); + + // Gaussian cut-off. + m_nsigma = get(cfg, "nsigma", m_nsigma); + + // Additional smearing. + m_smear_long = get_n(cfg["smear_long"]); + m_smear_tran = get_n(cfg["smear_tran"]); + + // Arbitrary time subtracted from window_start when setting frame time + m_reftime = get(cfg, "reference_time", m_reftime); + + // Arbitrary time added to tbin of traces on per plane basis. + m_tick_offsets.clear(); + m_tick_offsets.resize(3,0); + auto jto = cfg["time_offsets"]; + if (jto.isArray()) { + if (jto.size() == 3) { + for (int ind = 0; ind < 3; ++ind) { + m_tick_offsets[ind] = jto[ind].asDouble() / wtick; + } + } + else if (!jto.empty()) { + THROW(ValueError() << errmsg{"DepoFluxSplat: time_offsets must be empty or be a 3-array"}); + } + } + log->debug("speed={} mm/us, origin={} mm, tbins: {} {}us ticks: [{},{}]us, reftime={} us, tick offsets=({},{},{})", + m_speed / (units::mm/units::us), m_origin/units::mm, + m_tbins.nbins(), m_tbins.binsize()/units::us, + m_tbins.min()/units::us, m_tbins.max()/units::us, + m_reftime/units::us, + m_tick_offsets[0],m_tick_offsets[1],m_tick_offsets[2]); + +} + +IAnodeFace::pointer Gen::DepoFluxSplat::find_face(const IDepo::pointer& depo) +{ + for (auto face : m_anode->faces()) { + auto bb = face->sensitive(); + if (bb.inside(depo->pos())) { return face; } + } + return nullptr; +} + +// Return intersection of two half open ranges. +using intrange_t = std::pair; + +// Return intersection of half-open ranges r1 and r2. If r1 is fully before r2, +// return [r2.first,r2.first] and if fully after return [r2.second.r2.second]. +static intrange_t intersect(intrange_t const& r1, intrange_t const& r2) +{ + if (r1.second <= r2.first) { + return std::make_pair(r2.first, r2.first); + } + if (r1.first >= r2.second) { + return std::make_pair(r2.second, r2.second); + } + return std::make_pair(std::max(r1.first, r2.first), + std::min(r1.second, r2.second)); +} + +// A base class for sparse vs dense accumulation +struct Accumulator { + + virtual ~Accumulator() {}; + // Accumulate one depo's patch + virtual void add(int chid, int tbin, const std::vector& charge) = 0; + virtual IFrame::pointer frame(int ident, double time, double tick) = 0; + virtual size_t ntraces() const = 0; +}; + + +// Accumulate in a sparse way. +struct SparseAccumulator : public Accumulator { + ITrace::vector itraces; + + virtual ~SparseAccumulator() {} + virtual void add(int chid, int tbin, const std::vector& charge) { + itraces.push_back(std::make_shared(chid,tbin,charge)); + } + virtual IFrame::pointer frame(int ident, double time, double tick) { + return std::make_shared(ident, time, itraces, tick); + } + virtual size_t ntraces() const { return itraces.size(); } +}; + +static intrange_t make_intrange(int beg, size_t siz) +{ + return std::make_pair(beg, beg+siz); +} +static intrange_t union_intrange(intrange_t const& r1, intrange_t const& r2) +{ + return std::make_pair(std::min(r1.first, r2.first), + std::min(r1.second, r2.second)); +} + +// Accumulate in a dense way. +struct DenseAccumulator : public Accumulator { + using SharedSimpleTrace = std::shared_ptr; + std::map traces; // map chid->trace + virtual void add(int chid, int tbin, const std::vector& charge) { + auto tp = traces[chid]; + if (!tp) { // first seen + traces[chid] = std::make_shared(chid,tbin,charge); + return; + } + const auto& oldcharge = tp->charge(); + auto have = make_intrange(tp->tbin(), oldcharge.size()); + auto want = make_intrange(tbin, charge.size()); + auto need = union_intrange(have, want); + if (need.first < have.first || need.second > have.second) { + std::vector newcharge(need.second-need.first, 0); + std::copy(oldcharge.begin(), oldcharge.end(), + newcharge.begin() + have.first-need.first); + std::transform(charge.begin(), charge.end(), + newcharge.begin() + want.first-need.first, + newcharge.begin() + want.first-need.first, + std::plus()); + tp->charge() = newcharge; + } + } + + virtual IFrame::pointer frame(int ident, double time, double tick) { + ITrace::vector itraces; + for (const auto& [chid, tp] : traces) { + itraces.push_back(tp); + } + return std::make_shared(ident, time, itraces, tick); + } + virtual size_t ntraces() const { return traces.size(); } +}; + + + +bool Gen::DepoFluxSplat::operator()(const input_pointer& in, output_pointer& out) +{ + out = nullptr; + if (!in) { + log->debug("EOS at {}", m_count); + ++m_count; + return true; // EOS + } + + std::unique_ptr accum; + if (m_sparse) { + accum = std::make_unique(); + } + else { + accum = std::make_unique(); + } + + size_t ndepos_seen=0; + size_t ndepos_skipped=0; + size_t nplanes_skipped=0; + for (const auto& depo : *in->depos()) { + if (!depo) { + ++ndepos_skipped; + continue; + } + + auto face = find_face(depo); + if (!face) { + ++ndepos_skipped; + continue; + } + + ++ndepos_seen; + + // Depo is at response plane. Find its time at the collection + // plane assuming it were to continue along a uniform field. + // After this, all times are nominal up until we add arbitrary + // time offsets at output. + const double nominal_depo_time = depo->time() + m_origin / m_speed; + + // Tabulate depo flux for wire regions from each plane + for (auto plane : face->planes()) { + + int iplane = plane->planeid().index(); + if (iplane < 0) { + ++nplanes_skipped; + continue; + } + + // Allow for extra smear in time unique to each plane. + double sigma_L = depo->extent_long(); // [length] + const double smear_long = m_smear_long[iplane]; + if (smear_long > 0) { + const double extra = smear_long * m_tbins.binsize() * m_speed; + sigma_L = sqrt(sigma_L * sigma_L + extra * extra); + } + Gen::GausDesc time_desc(nominal_depo_time, sigma_L / m_speed); + + // Check if patch is fully outside time binning + { + double nmin_sigma = time_desc.distance(m_tbins.min()); + double nmax_sigma = time_desc.distance(m_tbins.max()); + + double eff_nsigma = depo->extent_long() > 0 ? m_nsigma : 0; + if (nmin_sigma > eff_nsigma || nmax_sigma < -eff_nsigma) { + ++ndepos_skipped; + continue; + } + } + + const Pimpos* pimpos = plane->pimpos(); + auto& wires = plane->wires(); + auto wbins = pimpos->region_binning(); // wire binning + + double sigma_T = depo->extent_tran(); + const double smear_tran = m_smear_tran[iplane]; + if (smear_tran > 0) { + const double extra = smear_tran * wbins.binsize(); + sigma_T = sqrt(sigma_T * sigma_T + extra * extra); + } + + const double center_pitch = pimpos->distance(depo->pos()); + Gen::GausDesc pitch_desc(center_pitch, sigma_T); + { + double nmin_sigma = pitch_desc.distance(wbins.min()); + double nmax_sigma = pitch_desc.distance(wbins.max()); + + double eff_nsigma = depo->extent_tran() > 0 ? m_nsigma : 0; + if (nmin_sigma > eff_nsigma || nmax_sigma < -eff_nsigma) { + // We are more than "N-sigma" outside the wire plane. + break; + } + } + + // The heavy lifting + Gen::GaussianDiffusion gd(depo, time_desc, pitch_desc); + gd.set_sampling(m_tbins, wbins, m_nsigma, 0, 1); + + // Transfer depo's patch to itraces + const auto patch = gd.patch(); // 2D array + + + // The absolute pitch bin for the first row of the patch array. + const int pbin0 = gd.poffset_bin(); + + // Absolute pitch bin range truncated to fit wire plane. + const auto p_range = intersect({pbin0, pbin0 + patch.rows()}, + {0, wbins.nbins()}); + if (p_range.first == p_range.second) { + ++nplanes_skipped; + continue; + } + + // The absolute tick bin for the first column of the patch array + const int tbin0 = gd.toffset_bin(); + + // Absolute tick bin range truncated to fit time window. + const auto t_range = intersect({tbin0, tbin0 + patch.cols()}, + {0, m_tbins.nbins()}); + if (t_range.first == t_range.second) { + ++nplanes_skipped; + continue; + } + + // Iterate over the valid wires covered by the patch + for (int pbin : irange(p_range)) { + auto iwire = wires[pbin]; + const int chid = iwire->channel(); + + const int ncharges = t_range.second - t_range.first; + std::vector charge(ncharges); + + const int prel = pbin - pbin0; + const int trel = t_range.first - tbin0; + + // truly cursed + Eigen::VectorXf::Map(&charge[0], ncharges) = patch.row(prel).segment(trel, charge.size()); + + // Differing conventions exist for the sign of the charge of + // "number of electrons" in the depo. Force positive signal. + std::transform(charge.begin(), charge.end(), charge.begin(), + static_cast(&std::abs)); + + accum->add(chid, t_range.first + m_tick_offsets[iplane], charge); + } // wires + } // plane + } // depos + + out = accum->frame(in->ident(), m_tbins.min() - m_reftime, m_tbins.binsize()); + log->debug("splat {} ndepos={}/{}/[{}] ntraces={}", + out->ident(), ndepos_seen, in->depos()->size(), nplanes_skipped, accum->ntraces()); + ++m_count; + return true; +} + + +// Local Variables: +// mode: c++ +// c-basic-offset: 4 +// End: diff --git a/gen/src/DepoSetDrifter.cxx b/gen/src/DepoSetDrifter.cxx index 6eda0dbf9..5d4fd6cc6 100644 --- a/gen/src/DepoSetDrifter.cxx +++ b/gen/src/DepoSetDrifter.cxx @@ -64,7 +64,7 @@ bool DepoSetDrifter::operator()(const input_pointer& in, output_pointer& out) // The EOS comes through all_depos.pop_back(); - log->debug("call={} drifted ndepos={} in={} out={}", m_count, all_depos.size(), charge_in, charge_out); + log->debug("call={} drifted ndepos={} Qout={} ({}%)", m_count, all_depos.size(), charge_out, 100.0*charge_out/charge_in); out = std::make_shared(m_count, all_depos); ++m_count; return true; diff --git a/gen/src/DepoTransform.cxx b/gen/src/DepoTransform.cxx index 629798f84..6537ff29a 100644 --- a/gen/src/DepoTransform.cxx +++ b/gen/src/DepoTransform.cxx @@ -94,6 +94,10 @@ void Gen::DepoTransform::configure(const WireCell::Configuration& cfg) m_drift_speed = get(cfg, "drift_speed", m_drift_speed); m_frame_count = get(cfg, "first_frame_number", m_frame_count); + log->debug("tick={} us, start={} us, readin={} us, drift_speed={} mm/us", + m_tick/units::us, m_start_time/units::us, + m_readout_time/units::us, m_drift_speed/(units::mm/units::us)); + auto jpirs = cfg["pirs"]; if (jpirs.isNull() or jpirs.empty()) { std::string msg = "must configure with some plane impact response components"; @@ -118,11 +122,11 @@ WireCell::Configuration Gen::DepoTransform::default_configuration() const put(cfg, "fluctuate", false); /// The open a gate. This is actually a "readin" time measured at - /// the input ("reference") plane. + /// the input ("response") plane. put(cfg, "start_time", m_start_time); /// The time span for each readout. This is actually a "readin" - /// time span measured at the input ("reference") plane. + /// time span measured at the input ("response") plane. put(cfg, "readout_time", m_readout_time); /// The sample period @@ -202,6 +206,9 @@ bool Gen::DepoTransform::operator()(const input_pointer& in, output_pointer& out auto trace = make_shared(chid, tbin, charge); traces.push_back(trace); } + // fixme: use SPDLOG_LOGGER_DEBUG + log->debug("plane={} face={} depos={} total traces={}", + iplane, face->ident(), face_depos.size(), traces.size()); } } diff --git a/gen/src/Drifter.cxx b/gen/src/Drifter.cxx index 7b07bf20e..090871968 100644 --- a/gen/src/Drifter.cxx +++ b/gen/src/Drifter.cxx @@ -111,7 +111,12 @@ void Gen::Drifter::configure(const WireCell::Configuration& cfg) THROW(ValueError() << errmsg{"no xregions given"}); } for (auto jone : jxregions) { - m_xregions.push_back(Xregion(jone)); + Xregion xr(jone); + m_xregions.push_back(xr); + log->debug("xregion: anode: {} mm, response: {} mm, cathode: {} mm", + xr.anode / units::mm, + xr.response / units::mm, + xr.cathode / units::mm); } log->debug("time offset: {} ms, drift speed: {} mm/us", m_toffset / units::ms, m_speed / (units::mm / units::us)); @@ -184,7 +189,8 @@ bool Gen::Drifter::insert(const input_pointer& depo) dT = sqrt(2.0 * m_DT * dt + dT * dT); } - auto newdepo = make_shared(depo->time() + direction * dt + m_toffset, pos, Qf, depo, dL, dT); + auto newdepo = make_shared(depo->time() + direction * dt + m_toffset, + pos, Qf, depo, dL, dT, depo->id()); xrit->depos.insert(newdepo); return true; } @@ -195,8 +201,6 @@ bool by_time(const IDepo::pointer& lhs, const IDepo::pointer& rhs) { return lhs- void Gen::Drifter::flush(output_queue& outq) { for (auto& xr : m_xregions) { - log->debug("xregion: anode: {} mm, response: {} mm, cathode: {} mm, flushing {}", xr.anode / units::mm, - xr.response / units::mm, xr.cathode / units::mm, xr.depos.size()); outq.insert(outq.end(), xr.depos.begin(), xr.depos.end()); xr.depos.clear(); } diff --git a/gen/src/DuctorFramer.cxx b/gen/src/DuctorFramer.cxx new file mode 100644 index 000000000..57ca4bfa1 --- /dev/null +++ b/gen/src/DuctorFramer.cxx @@ -0,0 +1,79 @@ +#include "WireCellGen/DuctorFramer.h" +#include "WireCellUtil/NamedFactory.h" + +WIRECELL_FACTORY(DuctorFramer, + WireCell::Gen::DuctorFramer, + WireCell::INamed, + WireCell::IDepoFramer, + WireCell::IConfigurable) + +using namespace WireCell; + +Gen::DuctorFramer::DuctorFramer() + : Aux::Logger("DuctorFramer", "gen") +{ +} + +Gen::DuctorFramer::~DuctorFramer() +{ +} + +WireCell::Configuration Gen::DuctorFramer::default_configuration() const +{ + Configuration cfg; + // The IDuctor that does the transformation + cfg["ductor"] = Json::nullValue; + // The IFrameFanin that merges frames in the case that the IDuctor returns + // more than one over a given depo set. + cfg["fanin"] = Json::nullValue; + return cfg; +} + +void Gen::DuctorFramer::configure(const WireCell::Configuration& cfg) +{ + auto ductor_tn = get(cfg, "ductor", ""); + m_ductor = Factory::find_tn(ductor_tn); + auto fanin_tn = get(cfg, "fanin", ""); + m_fanin = Factory::find_tn(fanin_tn); +} + + +bool Gen::DuctorFramer::operator()(const input_pointer& in, output_pointer& out) +{ + out=nullptr; + if (!in) { + log->debug("EOS at call={}", m_count); + ++m_count; + return true; + } + + IDepo::vector in_depos(in->depos()->begin(), in->depos()->end()); + in_depos.push_back(nullptr); // input EOS + + IFrameFanin::input_vector all_frames; + + for (auto idepo : in_depos) { + IDuctor::output_queue more; + (*m_ductor)(idepo, more); + all_frames.insert(all_frames.end(), more.begin(), more.end()); + } + + // The EOS should come through + if (all_frames.back()) { + raise("ductor did not pass through an EOS"); + } + all_frames.pop_back(); + + const IFrameFanin::input_vector& call_frames = all_frames; + IFrameFanin::output_pointer op; + bool ok = (*m_fanin)(call_frames, op); + if (ok) { + out = op; + } + else { + log->warn("DuctorFramer fanin failed"); + } + ++m_count; + return true; +} + diff --git a/gen/src/FrameFanin.cxx b/gen/src/FrameFanin.cxx index a91b1d19b..7823839ab 100644 --- a/gen/src/FrameFanin.cxx +++ b/gen/src/FrameFanin.cxx @@ -14,7 +14,6 @@ using namespace WireCell; Gen::FrameFanin::FrameFanin(size_t multiplicity) : Aux::Logger("FrameFanin", "glue") - , m_multiplicity(multiplicity) { } Gen::FrameFanin::~FrameFanin() {} @@ -22,6 +21,9 @@ Gen::FrameFanin::~FrameFanin() {} WireCell::Configuration Gen::FrameFanin::default_configuration() const { Configuration cfg; + + // The number of input ports. If set to zero the fanin will have DYNAMIC + // MULTIPLICITY and simply accept what it is given. cfg["multiplicity"] = (int) m_multiplicity; // A non-null entry in this array is taken as a string and used to @@ -29,6 +31,10 @@ WireCell::Configuration Gen::FrameFanin::default_configuration() const // they are placed to the output frame. This creates a tagged // subset of output traces corresponding to the entire set of // traces in the input frame from the associated port. + // + // When DYNAMIC MULTIPLICITY is used (multiplicity=0) the last tag in this + // array is applied to the all ports numbered greater or equal to the size + // of the tag array. cfg["tags"] = Json::arrayValue; // Tag rules are an array, one element per input port. Each @@ -36,24 +42,45 @@ WireCell::Configuration Gen::FrameFanin::default_configuration() const // their values are an object keyed by a regular expression // (regex) and with values that are a single tag or an array of // tags. See util/test_tagrules for examples. + // + // When DYNAMIC MULTIPLICITY (multiplicity=0) the last tag rule in this + // array is applied to all ports numbered greater or equal to the size of + // the tag rules array. cfg["tag_rules"] = Json::arrayValue; + // See also comments at top of FrameFanin.h + return cfg; } +std::string Gen::FrameFanin::tag(size_t port) const +{ + if (m_tags.empty()) return ""; + if (port < m_tags.size()) return m_tags.at(port); + return m_tags.back(); +} + void Gen::FrameFanin::configure(const WireCell::Configuration& cfg) { int m = get(cfg, "multiplicity", (int) m_multiplicity); - if (m <= 0) { + if (m < 0) { log->critical("illegal multiplicity: {}", m); THROW(ValueError() << errmsg{"FrameFanin multiplicity must be positive"}); } + if (m == 0) { + log->debug("FrameFanin using dynamic multiplicity"); + } + m_multiplicity = m; - m_tags.resize(m); - // Tag entire input frame worth of traces in the output frame. - auto jtags = cfg["tags"]; - for (int ind = 0; ind < m; ++ind) { - m_tags[ind] = convert(jtags[ind], ""); + // If tag is not empty string it will be used to tag the traces from the + // corresponding input port. + for (const auto& jtag : cfg["tags"]) { + if (jtag.isNull()) { + m_tags.push_back(""); + } + else { + m_tags.push_back(jtag.asString()); + } } // frame/trace tagging at the port level @@ -62,6 +89,9 @@ void Gen::FrameFanin::configure(const WireCell::Configuration& cfg) std::vector Gen::FrameFanin::input_types() { + if (!m_multiplicity) { + return std::vector(); + } const std::string tname = std::string(typeid(input_type).name()); std::vector ret(m_multiplicity, tname); return ret; @@ -83,12 +113,13 @@ bool Gen::FrameFanin::operator()(const input_vector& invec, output_pointer& out) return true; } - if (invec.size() != m_multiplicity) { - log->critical("input vector size={} my multiplicity={}", invec.size(), m_multiplicity); + const size_t multiplicity = invec.size(); + if (m_multiplicity && multiplicity != m_multiplicity) { + log->critical("input vector size={} my multiplicity={}", multiplicity, m_multiplicity); THROW(ValueError() << errmsg{"input vector size mismatch"}); } - std::vector by_port(m_multiplicity); + std::vector by_port(multiplicity); std::vector > stash; @@ -103,7 +134,7 @@ bool Gen::FrameFanin::operator()(const input_vector& invec, output_pointer& out) // creating an empty name map; only want combine the two maskmaps from each APA into a single maskmap, NOT merge masks std::map empty_name_map; - for (size_t iport = 0; iport < m_multiplicity; ++iport) { + for (size_t iport = 0; iport < multiplicity; ++iport) { const size_t trace_offset = out_traces.size(); const auto& fr = invec[iport]; @@ -120,7 +151,8 @@ bool Gen::FrameFanin::operator()(const input_vector& invec, output_pointer& out) // tag based on user rules. auto fintags = fr->frame_tags(); fintags.push_back(""); - auto fo = m_ft.transform(iport, "frame", fintags); + const size_t n = m_ft.nrules("frame"); + auto fo = m_ft.transform(iport >= n ? n-1 : iport, "frame", fintags); fouttags.insert(fo.begin(), fo.end()); } @@ -128,7 +160,8 @@ bool Gen::FrameFanin::operator()(const input_vector& invec, output_pointer& out) // collect transformed trace tags, any trace summary and their // offset trace indeices. annoying_factor *= 10; for (auto inttag : fr->trace_tags()) { - tagrules::tagset_t touttags = m_ft.transform(iport, "trace", inttag); + const size_t n = m_ft.nrules("frame"); + tagrules::tagset_t touttags = m_ft.transform(iport >= n ? n-1 : iport, "trace", inttag); if (touttags.empty()) { continue; } @@ -144,7 +177,7 @@ bool Gen::FrameFanin::operator()(const input_vector& invec, output_pointer& out) } }; - if (!m_tags[iport].empty()) { + if (tag(iport).size()) { IFrame::trace_list_t tl(traces->size()); std::iota(tl.begin(), tl.end(), trace_offset); by_port[iport] = tl; @@ -155,9 +188,9 @@ bool Gen::FrameFanin::operator()(const input_vector& invec, output_pointer& out) } auto sf = new Aux::SimpleFrame(one->ident(), one->time(), out_traces, one->tick(), out_cmm); - for (size_t iport = 0; iport < m_multiplicity; ++iport) { - if (m_tags[iport].size()) { - sf->tag_traces(m_tags[iport], by_port[iport]); + for (size_t iport = 0; iport < multiplicity; ++iport) { + if (tag(iport).size()) { + sf->tag_traces(tag(iport), by_port[iport]); } } diff --git a/gen/src/Reframer.cxx b/gen/src/Reframer.cxx index 8454a1bc9..047f15abd 100644 --- a/gen/src/Reframer.cxx +++ b/gen/src/Reframer.cxx @@ -21,10 +21,6 @@ using WireCell::Aux::SimpleFrame; Gen::Reframer::Reframer() : Aux::Logger("Reframer", "gen") - , m_toffset(0.0) - , m_fill(0.0) - , m_tbin(0) - , m_nticks(0) { } @@ -39,11 +35,12 @@ WireCell::Configuration Gen::Reframer::default_configuration() const // tags to find input traces/frames cfg["tags"] = Json::arrayValue; // tag to apply to output frame - cfg["frame_tag"] = ""; + cfg["frame_tag"] = m_frame_tag; cfg["tbin"] = m_tbin; cfg["nticks"] = m_nticks; cfg["toffset"] = m_toffset; cfg["fill"] = m_fill; + cfg["ignore_tags"] = m_ignore_tags; return cfg; } @@ -58,11 +55,15 @@ void Gen::Reframer::configure(const WireCell::Configuration& cfg) for (auto jtag : cfg["tags"]) { m_input_tags.push_back(jtag.asString()); } - m_frame_tag = get(cfg, "frame_tag", ""); + m_frame_tag = get(cfg, "frame_tag", m_frame_tag); m_toffset = get(cfg, "toffset", m_toffset); m_tbin = get(cfg, "tbin", m_tbin); m_fill = get(cfg, "fill", m_fill); m_nticks = get(cfg, "nticks", m_nticks); + m_ignore_tags = get(cfg, "ignore_tags", m_ignore_tags); + if ( m_ignore_tags && m_input_tags.size() ) { + raise("providing and ignoring tags is not consistent"); + } } std::pair Gen::Reframer::process_one(const ITrace::vector& itraces, const IFrame::trace_summary_t& isummary) { @@ -149,6 +150,16 @@ bool Gen::Reframer::operator()(const input_pointer& inframe, output_pointer& out std::unordered_map< std::string, IFrame::trace_summary_t> tag_summary; if (m_input_tags.empty()) { + auto ttags = inframe->trace_tags(); + if (! ttags.empty()) { + std::stringstream ss; + for (const auto& tag : ttags) { + ss << " " << tag; + } + if (! m_ignore_tags) { + log->warn("will combine traces from {} trace tags in frame:{}", ttags.size(), ss.str()); + } + } out_traces = process_one(*(inframe->traces())); } else { diff --git a/iface/inc/WireCellIface/IFrameFanin.h b/iface/inc/WireCellIface/IFrameFanin.h index eb270663c..fda82ea7f 100644 --- a/iface/inc/WireCellIface/IFrameFanin.h +++ b/iface/inc/WireCellIface/IFrameFanin.h @@ -10,9 +10,11 @@ namespace WireCell { * inputs is left to the implementation. */ class IFrameFanin : public IFaninNode { - public: + public: virtual ~IFrameFanin(); + typedef std::shared_ptr pointer; + virtual std::string signature() { return typeid(IFrameFanin).name(); } // Subclass must implement: diff --git a/pgraph/inc/WireCellPgraph/Node.h b/pgraph/inc/WireCellPgraph/Node.h index 08dc622cf..cbd05c310 100644 --- a/pgraph/inc/WireCellPgraph/Node.h +++ b/pgraph/inc/WireCellPgraph/Node.h @@ -58,6 +58,22 @@ namespace WireCell { return true; } + std::vector disconnected_ports() + { + std::vector ret; + for (auto& p : input_ports()) { + if (!p.edge()) { + ret.push_back(p); + } + } + for (auto& p : output_ports()) { + if (!p.edge()) { + ret.push_back(p); + } + } + return ret; + } + protected: // Concrete class should fill during construction PortList m_ports[Port::ntypes]; diff --git a/pgraph/inc/WireCellPgraph/Port.h b/pgraph/inc/WireCellPgraph/Port.h index bd5920973..b8071298d 100644 --- a/pgraph/inc/WireCellPgraph/Port.h +++ b/pgraph/inc/WireCellPgraph/Port.h @@ -32,18 +32,18 @@ namespace WireCell { Port(Node* node, Type type, std::string signature, std::string name = ""); - bool isinput(); - bool isoutput(); - Edge edge(); + bool isinput() const; + bool isoutput() const; + Edge edge() const; // Connect an edge, returning any previous one. Edge plug(Edge edge); // return edge queue size or 0 if no edge has been plugged - size_t size(); + size_t size() const; // Return true if queue is empty or no edge has been plugged. - bool empty(); + bool empty() const; // Get the next data. By default this pops the data off // the queue. To "peek" at the data, pas false. @@ -54,9 +54,12 @@ namespace WireCell { // Get back the associated Node. Node* node(); + const Node* node() const; - const std::string& name(); - const std::string& signature(); + const std::string& name() const; + const std::string& signature() const; + + std::string ident() const; private: Node* m_node; // node to which port belongs diff --git a/pgraph/src/Graph.cxx b/pgraph/src/Graph.cxx index 26d300e9f..9375f2ced 100644 --- a/pgraph/src/Graph.cxx +++ b/pgraph/src/Graph.cxx @@ -39,8 +39,9 @@ bool Graph::connect(Node* tail, Node* head, size_t tpind, size_t hpind) m_edges_forward[tail].push_back(head); m_edges_backward[head].push_back(tail); - SPDLOG_LOGGER_TRACE(l, "connect {}:({}:{}) -> {}({}:{})", tail->ident(), demangle(tport.signature()), tpind, - head->ident(), demangle(hport.signature()), hpind); + // was at in TRACE macro, maybe deserves a bit more prominence + l->debug("connect {}:({}:{}) -> {}({}:{})", tail->ident(), demangle(tport.signature()), tpind, + head->ident(), demangle(hport.signature()), hpind); return true; } @@ -154,12 +155,17 @@ bool Graph::call_node(Node* node) bool Graph::connected() { + bool okay = true; for (auto n : m_nodes) { - if (!n->connected()) { - return false; + auto bad = n->disconnected_ports(); + if (bad.empty()) continue; + okay = false; + l->warn("disconnected node: {}", n->ident()); + for (const auto& p : bad) { + l->warn("\tport: {}", p.ident()); } } - return true; + return okay; } void Graph::print_timers() const diff --git a/pgraph/src/Pgrapher.cxx b/pgraph/src/Pgrapher.cxx index 1056213c9..4016a87f6 100644 --- a/pgraph/src/Pgrapher.cxx +++ b/pgraph/src/Pgrapher.cxx @@ -24,13 +24,28 @@ static std::pair get_node(WireCell::Configuration // We should NOT be the one creating this component. auto nptr = WireCell::Factory::find_maybe_tn(node); if (!nptr) { - THROW(ValueError() << errmsg{"failed to get node"}); + raise("failed to get node \"%s\"", node); } int port = get(jone, "port", 0); return std::make_pair(nptr, port); } +static +std::string port_to_string(WireCell::Configuration node) +{ + std::stringstream ss; + ss << node["node"].asString() << ":" << node["port"].asString(); + return ss.str(); +} + +static +std::string edge_to_string(WireCell::Configuration tail, + WireCell::Configuration head) +{ + return port_to_string(tail) + " -> " + port_to_string(head); +} + void Pgrapher::configure(const WireCell::Configuration& cfg) { @@ -45,12 +60,12 @@ void Pgrapher::configure(const WireCell::Configuration& cfg) bool ok = m_graph.connect(fac(tail.first), fac(head.first), tail.second, head.second); if (!ok) { log->critical("failed to connect edge: {}", jedge); - THROW(ValueError() << errmsg{"failed to connect edge"}); + raise("failed to connect edge %s", edge_to_string(jedge["tail"], jedge["head"])); } } if (!m_graph.connected()) { log->critical("graph not fully connected"); - THROW(ValueError() << errmsg{"graph not fully connected"}); + raise("graph not fully connected"); } } diff --git a/pgraph/src/Port.cxx b/pgraph/src/Port.cxx index 9c83cfba4..1ee814a93 100644 --- a/pgraph/src/Port.cxx +++ b/pgraph/src/Port.cxx @@ -1,9 +1,13 @@ #include "WireCellPgraph/Port.h" +#include "WireCellUtil/Type.h" + +#include #include using namespace std; using namespace WireCell::Pgraph; +using WireCell::demangle; Port::Port(Node* node, Type type, std::string signature, std::string name) : m_node(node) @@ -14,10 +18,10 @@ Port::Port(Node* node, Type type, std::string signature, std::string name) { } -bool Port::isinput() { return m_type == Port::input; } -bool Port::isoutput() { return m_type == Port::output; } +bool Port::isinput() const { return m_type == Port::input; } +bool Port::isoutput() const { return m_type == Port::output; } -Edge Port::edge() { return m_edge; } +Edge Port::edge() const { return m_edge; } // Connect an edge, returning any previous one. Edge Port::plug(Edge edge) @@ -28,7 +32,7 @@ Edge Port::plug(Edge edge) } // return edge queue size or 0 if no edge has been plugged -size_t Port::size() +size_t Port::size() const { if (!m_edge) { return 0; @@ -37,7 +41,7 @@ size_t Port::size() } // Return true if queue is empty or no edge has been plugged. -bool Port::empty() +bool Port::empty() const { if (!m_edge or m_edge->empty()) { return true; @@ -77,5 +81,21 @@ void Port::put(Data& data) m_edge->push_back(data); } -const std::string& Port::name() { return m_name; } -const std::string& Port::signature() { return m_sig; } +const std::string& Port::name() const { return m_name; } +const std::string& Port::signature() const { return m_sig; } + +std::string Port::ident() const +{ + std::stringstream ss; + ss << ""; + return ss.str(); +} diff --git a/sigproc/src/OmnibusSigProc.cxx b/sigproc/src/OmnibusSigProc.cxx index f66223da0..4175de20a 100644 --- a/sigproc/src/OmnibusSigProc.cxx +++ b/sigproc/src/OmnibusSigProc.cxx @@ -733,6 +733,7 @@ void OmnibusSigProc::init_overall_response(IFrame::pointer frame) m_pad_nticks = m_fft_nticks - m_nticks; } + // Fixme: this should be moved into configure() auto ifr = Factory::find_tn(m_field_response); // Get full, "fine-grained" field responses defined at impact // positions. diff --git a/test/test/spdir-metric-gen-tests.sh b/test/test/spdir-metric-gen-tests.sh new file mode 100755 index 000000000..e1a7be988 --- /dev/null +++ b/test/test/spdir-metric-gen-tests.sh @@ -0,0 +1,13 @@ +#!/bin/bash + +# This is a tad kludgey but we want a main BATS file to run the "spdir-metric" +# test separately for each detector so waf will run each in parallel. But, we +# do not want to edit each one just to place the detector name where needed. +# So, we let sed do the editing. + +tdir="$(dirname "$(realpath "$BASH_SOURCE")")" + +for det in pdsp uboone +do + sed -e s/@DETECTOR@/$det/g < $tdir/spdir-metric.tmpl > $tdir/test-spdir-metric-${det}.bats +done diff --git a/test/test/spdir-metric-pdsp.bats b/test/test/spdir-metric-pdsp.bats new file mode 100644 index 000000000..4d8421f9d --- /dev/null +++ b/test/test/spdir-metric-pdsp.bats @@ -0,0 +1,5 @@ +#!/usr/bin/env bats +bats_load_library "wct-bats.sh" +@test "spdir metric $DETECTOR generate depos" { + echo "spdir metric $DETECTOR generate depos" 1>&3 +} diff --git a/test/test/spdir-metric.jsonnet b/test/test/spdir-metric.jsonnet new file mode 100644 index 000000000..9a2f35c41 --- /dev/null +++ b/test/test/spdir-metric.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.sim.signal(anode), + mid.sim.noise(anode), + mid.sim.digitizer(anode), + ], + sig: [ + mid.sim.signal(anode), + ], + noi: [ + mid.sim.noise(anode), + ], + dig: [ + mid.sim.digitizer(anode), + ], + nf: [ + mid.sigproc.nf(anode), + ], + sp: [ + mid.sigproc.sp(anode), + ] + }, + pre_sink(stage) :: + if stage == "splat" + then [ mid.sim.splat(anode) ] + else [], + + reframer(stage) :: + local reframers = { + splat: [mid.sim.reframer(anode, name=outputs[stage])], + sp: [mid.sim.reframer(anode, name=outputs[stage], tags=["gauss0"])], + }; + 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 stage == "drift" + 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 mid = high.mid(detector, variant); + 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/test/test/spdir-metric.org b/test/test/spdir-metric.org new file mode 100644 index 000000000..a4e29a893 --- /dev/null +++ b/test/test/spdir-metric.org @@ -0,0 +1,40 @@ +#+title: The SP directional metric test +#+setupfile: ../../docs/setup-topic.org + +* Overview + +The full "spdir-metric" test runs the following chain: + + 1) *Generate* ideal line tracks at different angles in the detector volume. (~wirecell-*~) + 2) *Simulate* noise and ADC waveform from the tracks. (~wire-cell~) + 3) *Translate* apply an optional ADC waveform transform. (external) + 4) *Filter* apply WCT noise filtering. (~wire-cell~) + 5) *Signals* apply WCT signal processing. (~wire-cell~) + 6) *Metrics* for each track, calculate and output the metrics. (~wirecell-*~) + 7) *Plots* consume all metrics and make standardized plots. (~wirecell-*~) + +Comments: + +- The main command used to execute each step is given in parenthesis. + +- The *translate* step is supported by the infrastructure for this test but is exercised outside the scope of WCT testing. With no translate stage, *simulate*, *filter* and *signals* stages can be run in a single call to ~wire-cell~. + +- The output of the *metrics* stage will be stored into the WCT test data repository to allow for historical tests. + +* Meta + + The source for this document is in org markup. You may read that form directly or you may use Emacs to /export/ the org content to HTML with ~C-c C-e h h~ or to PDF with ~C-c C-e l p~. The org content may also be /tangled/ with ~C-c C-v t~ to produce the individual BATS test files which of which executes the stages for a given supported detector. +Depending on where you are reading this document you may find these links provide it in the different formats: [[file:spdir-metric.html][HTML]], [[file:spdir-metric.pdf][PDF]]. + +* The WCT configuration + +A single top level ~wire-cell~ configuration file configures how to run the *simulate*, *filter* and *signals* stages. It requires specifying the detector, input and output file names and optional list of tasks. Eg, + +#+begin_example +wire-cell \ + -A detector=pdsp \ + -A input=depos.npz \ + -A output=sig.npz \ + test/test/spdir-metric.jsonnet +#+end_example + diff --git a/test/test/spdir-metric.tmpl b/test/test/spdir-metric.tmpl new file mode 100644 index 000000000..ca15b7e2d --- /dev/null +++ b/test/test/spdir-metric.tmpl @@ -0,0 +1,158 @@ +#!/usr/bin/env bats +# -*- shell-script -*- + +# Caution: if this file has .bats extension it is fully generated. +# Any edits are subject to loss. +# Instead, edit spdir-metric.tmpl and run spdir-metric-gen-tests.sh. + +bats_load_library "wct-bats.sh" + +DETECTOR=@DETECTOR@ + +# The detector variant. +VARIANT="spdir" + +function gen_fname () { + local tier="$1" ; shift + local num="$1" ; shift + local ext="${1:-npz}" + printf "%s-%s-%02.0f.%s" "$DETECTOR" "$tier" "$num" "$ext" +} + + +# The MicroBooNE conventional directions are defined in X-axis of figure 30 of +# the the MB SP paper with angle definitions in figure 29. +# +# Two big caveats in the MB convention: +# +# 1) It relies on U/V wire direction being mirror-symmetric about W. +# +# 2) It samples U/V wire direction space differently than W. Event when the +# "primed" (U/V) angle theta_xz' matches the W angle theta_xz, the theta_y' and +# theta_y differ. And so it is not trivial to make U/V vs W comparisons. +mb_theta_xz_deg=( 0 1 3 5 10 20 30 45 60 75 80 82 84 89 ) +mb_theta_y_deg=( 90 90 90 90 90 90 90 90 90 90 90 90 90 90 ) + +# In future, perhaps we pick others. For now, accept MB's +theta_xz_deg=( ${mb_theta_xz_deg[@]} ) +theta_y_deg=( ${mb_theta_y_deg[@]} ) +ndirs="${#theta_xz_deg[@]}" +directions=$(seq 0 $(( $ndirs - 1 ))) + +function generate_depos () +{ + local num="$1" ; shift + local theta_xz=${theta_xz_deg[$num]}'*deg' + local theta_y=${theta_y_deg[$num]}'*deg' + + echo "generate_depos: $num $theta_xz $theta_y" + + local fname="$(gen_fname depos $num)" + local pname="$(gen_fname params $num json)" + + # note: MB only select based on collection plane=2 + wcpy gen detlinegen \ + -d $DETECTOR --apa 0 --plane 2 --angle-coords=wire-plane \ + --theta_xz="$theta_xz" --theta_y="$theta_y" \ + -o "$fname" -m "$pname" + + # avoid empty zip of 22 bytes + file_larger_than "$fname" 23 +} + + +@test "spdir metric @DETECTOR@ generate depos" { + cd_tmp file + + for num in $directions + do + generate_depos $num + done +} + + +@test "spdir metric @DETECTOR@ sim+nf+sp one job graph viz" { + cd_tmp file + local cfg_file="$(relative_path spdir-metric.jsonnet)" + local num="${directions[0]}" + local depos="$(gen_fname depos "$num")" + # save out the intermediate drifted depos for input to metric + local drifts="$(gen_fname drifts "$num")" + local splats="$(gen_fname splats "$num")" + local signals="$(gen_fname signals "$num")" + local pdf="$(gen_fname graph "$num" pdf)" + local log="$(gen_fname wct "$num" log)" + wirecell-pgraph dotify \ + -A input="$depos" \ + -A output="{drift:\"$drifts\",splat:\"$splats\",sp:\"$signals\"}" \ + -A detector=$DETECTOR \ + -A variant="$VARIANT" \ + "$cfg_file" "$pdf" +} + +@test "spdir metric @DETECTOR@ sim+nf+sp run job" { + cd_tmp file + local cfg_file="$(relative_path spdir-metric.jsonnet)" + for num in ${directions[*]} + do + local depos="$(gen_fname depos "$num")" + # save out the intermediate drifted depos for input to metric + local drifts="$(gen_fname drifts "$num")" + local splats="$(gen_fname splats "$num")" + local signals="$(gen_fname signals "$num")" + local log="$(gen_fname wct "$num" log)" + wire-cell -c "$cfg_file" \ + -l "$log" -L debug \ + -A input="$depos" \ + --tla-code output="{drift:\"$drifts\",splat:\"$splats\",sp:\"$signals\"}" \ + -A detector=$DETECTOR \ + -A variant="$VARIANT" + + done +} + + +# FIXME +# +# We do not yet have a command to generate a metric from each frame file. +# When we do, fix this function. For now, it does some bogosity. +# +# This function should consume params.json +function calcualte_metrics () +{ + local jname="$(gen_fname metrics $num json)" + cat < $jname +{ + "num": $num + "theta": $num, + "phi": $num, + "nchan": 0, + "sigtot": 0, + "eff": [], + "bias": [], + "res": [] +} +EOF + # This metrics.json should produce rich but summary data. Eg, the + # eff/bias/res should be arrays the hit channels. +} +# FIXME + +@test "spdir metric @DETECTOR@ metrics" { + cd_tmp file + for num in ${directions[*]} + do + calcualte_metrics "$num" + done + +} + + +# FIXME: we also need a canonical plotter that consumes all metric.json. +@test "spdir metric @DETECTOR@ plots" { + cd_tmp file + for num in ${directions[*]} + do + ls -l "$(gen_fname metrics $num json)" >> "$(gen_fname plots all txt)" + done +} diff --git a/test/test/ssss-pdsp.py b/test/test/ssss-pdsp.py new file mode 100755 index 000000000..8c9a0cc46 --- /dev/null +++ b/test/test/ssss-pdsp.py @@ -0,0 +1,597 @@ +#!/usr/bin/env python +import numpy +import numpy.ma as ma + +import scipy +import matplotlib.pyplot as plt +import click +from pathlib import Path +from wirecell import units +from wirecell.util.cli import log, context +from wirecell.util.plottools import pages +from wirecell.util.peaks import ( + baseline_noise, select_activity, + find1d as find_peaks, + gauss as gauss_func +) + +@context("ssss-pdsp") +def cli(ctx): + ''' + The ssss-pdsp test + ''' + pass + + +@cli.command("log-summary") +@click.argument("log") +def depo_summary(log): + ''' + Print summary of wct log + ''' + want = " face:0 with 3 planes and sensvol: " + for line in open(log).readlines(): + if want in line: + n = line.find(": [(") + len(want) + off = len(want) + line.find(want) + 1 + parts = [p.strip() for p in line[off:-2].split(" --> ")] + p1 = numpy.array(list(map(float, parts[0][1:-1].split(" "))))/units.mm + p2 = numpy.array(list(map(float, parts[1][1:-1].split(" "))))/units.mm + print(f''' +#+caption: Sensitive volume diagonal and ideal line track endpoints. +#+name: tab:diagonal-endpoints +| x (mm) | y (mm) | z (mm) | +|--------|--------|--------| +| {p1[0]}| {p1[1]}| {p1[2]}| +| {p2[0]}| {p2[1]}| {p2[2]}| + ''') + return + +def load_depos(fname): + fp = numpy.load(fname) + d = fp["depo_data_0"] + i = fp["depo_info_0"] + ind = i[:,2] == 0 + return d[ind,:] + +@cli.command("depo-summary") +@click.argument("depos") +@click.argument("drift") +def depo_summary(depos, drift): + ''' + Print summary of depos + ''' + d1 = load_depos(depos) + d2 = load_depos(drift) + das = [d1,d2] + + tmins = numpy.array([numpy.min(d[:,0]) for d in das])/units.us + # print(f'{tmins=}') + + tmaxs = numpy.array([numpy.max(d[:,0]) for d in das])/units.us + # print(f'{tmaxs=}') + + xmins = numpy.array([numpy.min(d[:,2]) for d in das])/units.mm + xmaxs = numpy.array([numpy.max(d[:,2]) for d in das])/units.mm + + # [t/x][min/max][depos/drift] + + # copy-paste into org and hit tab + print(f''' +#+caption: Depo t/x ranges: +#+name: tab:depo-tx-ranges +| | min | max | units | +|---------|------|-----|-------| +| depos t | {tmins[0]:.2f}|{tmaxs[0]:.3f} | us | +| drift t | {tmins[1]:.2f}|{tmaxs[1]:.3f} | us | +| depos x | {xmins[0]:.2f}|{xmaxs[0]:.2f} | mm | +| drift x | {xmins[1]:.2f}|{xmaxs[1]:.2f} | mm | + ''') + +def load_frame(fname, tunit=units.us): + ''' + Load a frame with time values in explicit units. + ''' + fp = numpy.load(fname) + f = fp["frame_*_0"] + t = fp["tickinfo_*_0"] + c = fp["channels_*_0"] + + c2 = numpy.array(c) + numpy.sort(c2) + assert numpy.all(c == c2) + + cmin = numpy.min(c) + cmax = numpy.max(c) + nch = cmax-cmin+1 + ff = numpy.zeros((nch, f.shape[1]), dtype=f.dtype) + for irow, ch in enumerate(c): + ff[cmin-ch] = f[irow] + + t0 = t[0]/tunit + tick = t[1]/tunit + tf = t0 + f.shape[1]*tick + + # edge values + extent=(t0, tf+tick, cmin, cmax+1) + origin = "lower" # lower flips putting [0,0] at bottom + ff = numpy.flip(ff, axis=0) + return (ff, extent, origin) + +@cli.command("frame-summary") +@click.argument("files", nargs=-1) +def frame_summary(files): + data = dict() + for fname in files: + p = Path(fname) + data[p.stem] = load_frame(fname) + + + rowpat = '| {name} | {start} | {duration:.0f} | {npos} | {vmin:.1f} | {vmax:.1f} |' + rows = [ + '#+caption: Frame times in microsecond.', + '#+name: tab:frame times', + '| frame | start | duration | n>0 | min | max |', + ] + for name, (f,e,o) in sorted(data.items()): + p = dict( + name = name, + start = e[0], + duration = e[1]-e[0], + npos = numpy.sum(f > 0), + vmin = numpy.min(f), + vmax = numpy.max(f), + ) + row = rowpat.format(**p) + rows.append(row) + + print('\n'.join(rows)) + +from matplotlib.gridspec import GridSpec, GridSpecFromSubplotSpec + +def dump_frame(name,f): + print(name,f.shape,numpy.sum(f)) + +def plot_ticks(ax, t0, tf, drift, f1, f2, channel_ranges, speed = 1.565*units.mm/units.us): + nticks = f1.shape[1] + times = numpy.linspace(t0,tf,nticks+1, endpoint=True) + deltat = times[1]-times[0] + print(f'time binning: [{t0:.1f},{tf:1f}] us x {nticks}') + dump_frame('splat',f1) + dump_frame('signal',f2) + + print(f'warning: using speed of {speed/(units.mm/units.us)} mm/us') + ts = (drift[:,0] + 10*units.cm / speed)/units.us + qs = -1*drift[:,1] + + sigs = drift[:,5] / speed / units.us + nsig=5 + nspread = numpy.array((nsig/deltat)*sigs, dtype=int) + + gt = numpy.zeros(nticks) + for t,q,n in zip(ts,qs,nspread): + tbin = int((t-t0)/deltat) + gt[tbin-n:tbin+n+1] += q/(2*n+1) + + # h = numpy.histogram(ts, bins=times, weights=qs) + # ax.plot(h[1][:-1], h[0], label='depo charge') + + for p,chans in zip("UVW",channel_ranges): + val1 = f1[chans,:].sum(axis=0) + val2 = f2[chans,:].sum(axis=0) + ax.plot(times[:-1], val1, label=p+' splat') + ax.plot(times[:-1], val2, label=p+' signal') + val1 = f1.sum(axis=0) + val2 = f2.sum(axis=0) + ax.plot(times[:-1], val2/3.0, label='total signal / 3') + ax.plot(times[:-1], gt, label='depo charge') + ax.plot(times[:-1], val1/3.0, label='total splat / 3') + + + + # ax.set_xlim((1400,1600)) + # ax.set_xlim((0,200)) + ax.legend(loc='right') + + + +def plot_frame(gs, f, e, o="lower", channel_ranges=None, which="splat", tit=""): + t0,tf,c0,cf=e + + gs = GridSpecFromSubplotSpec(2,2, subplot_spec=gs, + height_ratios = [5,1], width_ratios = [6,1]) + + + fax = plt.subplot(gs[0,0]) + tax = plt.subplot(gs[1,0], sharex=fax) + cax = plt.subplot(gs[0,1], sharey=fax) + + cax.set_xlabel(which) + fax.set_ylabel("channel") + if which=="signal": + tax.set_xlabel("time [us]") + + if tit: + plt.title(tit) + plt.setp(fax.get_xticklabels(), visible=False) + plt.setp(cax.get_yticklabels(), visible=False) + if which=="splat": + plt.setp(tax.get_xticklabels(), visible=False) + + im = fax.imshow(f, extent=e, origin=o, aspect='auto', vmax=500, cmap='hot_r') + + tval = f.sum(axis=0) + t = numpy.linspace(e[0],e[1],f.shape[1]+1,endpoint=True) + tax.plot(t[:-1], tval) + if channel_ranges: + for p,chans in zip("UVW",channel_ranges): + val = f[chans,:].sum(axis=0) + c1 = chans.start + c2 = chans.stop + print(c1,c2,numpy.sum(val)) + tax.plot(t[:-1], val, label=p) + fax.plot([t0,tf], [c1,c1]) + fax.text(t0 + 0.1*(tf-t0), c1 + 0.5*(c2-c1), p) + fax.plot([t0,tf], [c2-1,c2-1]) + tax.legend() + + cval = f.sum(axis=1) + c = numpy.linspace(e[2],e[3],f.shape[0]+1,endpoint=True) + cax.plot(cval, c[:-1]) + + return im + +def plot_zoom(f1, f2, ts, cs, tit, o="lower"): + + e = (ts[0]*0.5, ts[1]*0.5, cs[0], cs[1]) + ts = slice(*ts) + cs = slice(*cs) + + f1 = f1[cs, ts] + f2 = f2[cs, ts] + + fig, axes = plt.subplots(3,1, sharex=True, sharey=True) + plt.suptitle(tit) + + vmax=2000 + args=dict(extent=e, origin=o, aspect='auto', vmin=-vmax, vmax=vmax, cmap='Spectral') + + im0 = axes[0].imshow(f1, **args) + im1 = axes[1].imshow(f2, **args) + im2 = axes[2].imshow(f1-f2, **args) + axes[2].set_xlabel('time [us]') + fig.subplots_adjust(right=0.85) + cbar_ax = fig.add_axes([0.88, 0.15, 0.04, 0.7]) + fig.colorbar(im1, cax=cbar_ax) + +def smear_tick(f1, sigma=3.0, nsigma=6): + print(f'smear: tot sig {numpy.sum(f1)}') + gx = numpy.arange(-sigma*nsigma, sigma*nsigma, 1) + if gx.size%2: + gx = gx[:-1] + gauss = numpy.exp(-0.5*(gx/sigma)**2) + gauss /= numpy.sum(gauss) + smear = numpy.zeros_like(f1[0]) + hngx = gx.size//2 + smear[:hngx] = gauss[hngx:] + smear[-hngx:] = gauss[:hngx] + nticks = f1[0].size + totals = numpy.sum(f1, axis=1) + for irow in range(f1.shape[0]): + wave = f1[irow] + tot = numpy.sum(wave) + wave = scipy.signal.fftconvolve(wave, smear)[:nticks] + f1[irow] = wave * tot / numpy.sum(wave) + return f1 + + +def plot_channels(f1, f2, chans, t0, t1, tit=""): + x = numpy.arange(t0, t1) + + if tit: + plt.title(tit) + for c in chans: + w1 = f1[c][t0:t1] + w2 = f2[c][t0:t1] + tot1 = numpy.sum(w1) + tot2 = numpy.sum(w2) + plt.plot(x, w1, label=f'ch{c} q={tot1:.0f} splat') + plt.plot(x, w2, label=f'ch{c} q={tot2:.0f} signal') + plt.legend() + plt.xlabel("ticks") + + +@cli.command("plots") +@click.option("--channel-ranges", default="0,800,1600,2560", + help="comma-separated list of channel idents defining ranges") +@click.option("--smear", default=0.0, help="Apply post-hoc smearing across tick dimension") +@click.option("--scale", default=0.0, help="Apply post-hoc multiplicative scale to waveforms") +@click.option("-o",'--output', default='plots.pdf') +@click.argument("depos") +@click.argument("drift") +@click.argument("splat") +@click.argument("signal") +def plots(channel_ranges, smear, scale, depos, drift, splat, signal, output, **kwds): + ''' + Make plots comparing {depos,drift,signal,splat}.npz + ''' + # dumb sanity check user gave us files in the right order + assert "depos" in depos + assert "drift" in drift + assert "splat" in splat + assert "signal" in signal + + + if channel_ranges: + channel_ranges = list(map(int,channel_ranges.split(","))) + channel_ranges = list(zip(channel_ranges[:-1], channel_ranges[1:])) + + fig = plt.figure() + + # d1 = load_depos(depos) + d2 = load_depos(drift) + f1,e1,o1 = load_frame(splat) + print(f'extent={e1}') + if smear: + f1 = smear_tick(f1, smear) + if scale: + f1 *= scale + + f2,e2,o2 = load_frame(signal) + + with pages(output) as out: + + # plt.suptitle(f'{smear=:.1f} {scale=:.1f}') + pgs = GridSpec(1,2, figure=fig, width_ratios = [7,0.2]) + gs = GridSpecFromSubplotSpec(2, 1, pgs[0,0]) + im1 = plot_frame(gs[0], f1, e1, o1, channel_ranges, which="splat") + im2 = plot_frame(gs[1], f2, e2, o2, channel_ranges, which="signal") + fig.colorbar(im2, cax=plt.subplot(pgs[0,1])) + plt.tight_layout() + out.savefig() + + plt.clf() + + fig, ax = plt.subplots(1,1, sharex=True, sharey=True) + plot_ticks(ax, e1[0], e1[1], d2, f1, f2, channel_ranges) + out.savefig() + + plt.clf() + plot_zoom(f1, f2, [0,400], [1400,1600], + f"splat - signal difference, begin of track, V-plane") + out.savefig() + + plt.clf() + plot_zoom(f1, f2, [4400, 4800], [1100, 1300], + f"splat - signal difference, end of track, V-plane") + out.savefig() + + plt.clf() + plot_zoom(f1, f2, [100,200], [1525,1560], + tit=f"splat - signal difference, begin of track, V-plane") + out.savefig() + + plt.clf() + plot_zoom(f1, f2, [4600, 4800], [1190, 1230], + tit=f"splat - signal difference, end of track, V-plane") + out.savefig() + + plt.clf() + plot_channels(f1, f2, [1540, 1530, 1520, 1510], 100, 350, + tit=f'V-plane start') + out.savefig() + + plt.clf() + plot_channels(f1, f2, [1200, 1210, 1220, 1230], 4550, 4750, + tit=f'V-plane end') + out.savefig() + +def relbias(a,b): + rb = numpy.zeros_like(a) + ok = b>0 + rb[ok] = a[ok]/b[ok] - 1 + return rb + +@cli.command("plot-ssss") +@click.option("--channel-ranges", default="0,800,1600,2560", + help="comma-separated list of channel idents defining ranges") +@click.option("--nsigma", default=3.0, + help="Relative threshold on signal in units of number of sigma of noise width") +@click.option("-o",'--output', default='plots.pdf') +@click.argument("depos") +@click.argument("drift") +@click.argument("splat") +@click.argument("signal") +def plot_ssss(channel_ranges, depos, drift, splat, signal, + nsigma, output, **kwds): + ''' + Plot splat and sim+SP signals. + ''' + # dumb sanity check user gave us files in the right order + assert "depos" in depos + assert "drift" in drift + assert "splat" in splat + assert "signal" in signal + + if channel_ranges: + channel_ranges = list(map(int,channel_ranges.split(","))) + channel_ranges = [slice(*cr) for cr in zip(channel_ranges[:-1], channel_ranges[1:])] + + fig = plt.figure() + + # d1 = load_depos(depos) + d2 = load_depos(drift) + + f1,e1,o1 = load_frame(splat) + f2,e2,o2 = load_frame(signal) + + tick_us = (e2[1] - e2[0])/(f1.shape[1]+1) + print(f'{tick_us=}') + + + with pages(output) as out: + + pgs = GridSpec(1,2, figure=fig, width_ratios = [7,0.2]) + gs = GridSpecFromSubplotSpec(2, 1, pgs[0,0]) + im1 = plot_frame(gs[0], f1, e1, o1, channel_ranges, which="splat") + im2 = plot_frame(gs[1], f2, e2, o2, channel_ranges, which="signal") + fig.colorbar(im2, cax=plt.subplot(pgs[0,1])) + plt.tight_layout() + out.savefig() + + byplane = list() + + # Per channel range plots. + for pln, ch in enumerate(channel_ranges): + letter = "UVW"[pln] + + spl = select_activity(f1, ch, nsigma) + sig = select_activity(f2, ch, nsigma) + print(f'selection: {ch} {spl.selection.shape} {sig.selection.shape}') + + # Find the bbox that bounds the biggest splat object. + biggest = spl.plats.sort_by("sums")[-1] + bbox = spl.plats.bboxes[biggest] + print(f'{bbox=}') + + byplane.append((spl, sig, bbox)) + + spl_act = spl.thresholded[bbox] + sig_act = sig.thresholded[bbox] + + print('shapes',spl_act.shape, sig_act.shape) + + + # bias of first w.r.t. second + bias1 = relbias(sig_act, spl_act) + bias2 = relbias(spl_act, sig_act) + + plt.clf() + fig, axes = plt.subplots(nrows=2, ncols=2, sharex=True, sharey=True) + plt.suptitle(f'{letter}-plane') + args=dict(aspect='auto') + im1 = axes[0,0].imshow(sig_act, **args) + fig.colorbar(im1, ax=axes[0,0]) + im2 = axes[0,1].imshow(spl_act, **args) + fig.colorbar(im2, ax=axes[0,1]) + + args = dict(args, cmap='jet', vmin=-50, vmax=50) + + im3 = axes[1,0].imshow(100*bias1, **args) + fig.colorbar(im3, ax=axes[1,0]) + + im4 = axes[1,1].imshow(100*bias2, **args) + fig.colorbar(im4, ax=axes[1,1]) + + axes[0,0].set_title(f'signal {nsigma=}') + axes[0,1].set_title(f'splat {nsigma=}') + + axes[1,0].set_title(f'splat/signal - 1 [%]') + axes[1,1].set_title(f'signal/splat - 1 [%]') + + chan_tit = 'chans (rel)' + tick_tit = 'ticks (rel)' + axes[0,0].set_ylabel(chan_tit) + axes[1,0].set_ylabel(chan_tit) + axes[1,0].set_xlabel(tick_tit) + axes[1,1].set_xlabel(tick_tit) + + fig.subplots_adjust(right=0.85) + + plt.tight_layout() + out.savefig() + + plt.clf() + plt.suptitle('By channel') + fig, axes = plt.subplots(nrows=2, ncols=3, sharey="row") + for pln, (spl, sig, bbox) in enumerate(byplane): + + spl_qch = numpy.sum(spl.activity[bbox], axis=1) + spl_tot = numpy.sum(spl_qch) + sig_qch = numpy.sum(sig.activity[bbox], axis=1) + sig_tot = numpy.sum(sig_qch) + + print('ch tots:', spl_qch.shape, sig_qch.shape) + + ### MD SP paper definitions: + # Charge bias - the mean fractional difference between the total + # deconvolved charge and the true charge. In terms of a track (line + # charge), the mean value of the distribution of each wire’s + # integrated deconvolved charges from a range of wires will be used + # to calculate the charge bias with respect to the true charge + # within one wire pitch. + # + # Charge resolution - the standard deviation of the total + # deconvolved charge relative to the mean deconvolved charge. In + # analogy to the definition of charge bias, in terms of a track + # (line charge), the distribution of each wire’s integrated + # deconvolved charges from a range of wires will apply. + # + # Charge inefficiency - specifically associated with tracks (line + # charges), the fraction of the wires which have ZERO deconvolved + # charge (no ROI found). Note that these wires will not be involved + # in the calculation of charge bias or charge resolution. + #### + + # However, while a channel that has no signal is "inefficient" it is + # possible for a channel to be "over efficient" by having signal + # where there is no splat. This has as much to do with (smeared) + # "splat" being an imperfect "true signal" as it has anything to do + # with "real" sim+SP signals. So, we interpret MB's "ROI" above to + # mean channels that have both signal and splat and then inefficient + # channels have only splat and over efficient channels have only + # signal. + + # either-or, exclude channels where both are zero + eor = numpy.logical_or (spl_qch > 0, sig_qch > 0) + # both are nonzero + both = numpy.logical_and(spl_qch > 0, sig_qch > 0) + nosig = numpy.logical_and(spl_qch > 0, sig_qch == 0) + nospl = numpy.logical_and(spl_qch == 0, sig_qch > 0) + + neor = numpy.sum(eor) + # inefficiency percent + nnosig = numpy.sum(nosig) + ineff = 100*nnosig/neor + # over efficiency percent + nnospl = numpy.sum(nospl) + oveff = 100*nnospl/neor + print(f'{neor=} {nnosig=} {nnospl=}') + + nbins = 50 + reldiff = 100.0*(spl_qch[both] - sig_qch[both])/spl_qch[both], + bln = baseline_noise(reldiff, nbins, nbins//2) + counts, edges = bln.hist + + # we plot: + # 1. The total charge on each channel for signal and splat + # 2. The histogram of their difference + + letter = "UVW"[pln] + + ax1,ax2 = axes[:,pln] + + ax1.plot(sig_qch, label='signal') + ax1.plot(spl_qch, label='splat') + ax1.set_xlabel('chans (rel)') + ax1.set_ylabel('electrons') + ax1.set_title(f'{letter} ineff={ineff:.1f}%') + ax1.legend() + + ax2.step(edges[:-1], counts, label='data') + model = gauss_func(edges[:-1], bln.A, bln.mu, bln.sigma) + ax2.plot(edges[:-1], model, label='fit') + ax2.set_title(f'mu={bln.mu:.2f}%\nsig={bln.sigma:.2f}%') + ax2.set_xlabel('difference [%]') + ax2.set_ylabel('counts') + ax2.legend() + + plt.suptitle('(splat - signal) / splat') + plt.tight_layout() + out.savefig() + +def main(): + cli() + +if '__main__' == __name__: + main() diff --git a/test/test/test-morse-pdsp.bats b/test/test/test-morse-pdsp.bats new file mode 100644 index 000000000..8a3e40a02 --- /dev/null +++ b/test/test/test-morse-pdsp.bats @@ -0,0 +1,47 @@ +#!/usr/bin/env bats + +bats_load_library "wct-bats.sh" + +setup_file () { + cd_tmp + echo 'depo-morse --wires pdsp --output depos.npz' > depo-morse-args-tmp + mv_if_diff depo-morse-args-tmp depo-morse-args +} + +@test "morse gen depos" { + cd_tmp file + read -a args < depo-morse-args + + run_idempotently -s depo-morse-args -t depos.npz -- \ + wirecell-gen "${args[@]}" + file_larger_than depos.npz 23 +} + +@test "morse run wct" { + cd_tmp file + local cfg_file="$(relative_path spdir-metric.jsonnet)" + + run_idempotently -s "$cfg_file" -s depos.npz -t wct.log -t drift.npz -t splat.npz -t digit.npz -t signal.npz -- \ + wire-cell \ + -l wct.log -L debug \ + -A input="depos.npz" \ + --tla-code output="{drift:\"drift.npz\",splat:\"splat.npz\",sim:\"digit.npz\",sp:\"signal.npz\"}" \ + -A detector=pdsp \ + -A variant=morse_nominal \ + -A tasks="drift,splat,sim,sp" \ + "$cfg_file" +} + +@test "morse make plots" { + cd_tmp file + + run_idempotently -s signal.npz -t peaks.json -- \ + wirecell-gen morse-summary -o peaks.json signal.npz + + run_idempotently -s peaks.json -t splat-params.json -- \ + wirecell-gen morse-splat -o splat-params.json peaks.json + + run_idempotently -s signal.npz -t peaks.pdf -- \ + wirecell-gen morse-plots -o peaks.pdf signal.npz + +} diff --git a/test/test/test-spdir-metric-pdsp.bats b/test/test/test-spdir-metric-pdsp.bats new file mode 100644 index 000000000..e27ca0e29 --- /dev/null +++ b/test/test/test-spdir-metric-pdsp.bats @@ -0,0 +1,158 @@ +#!/usr/bin/env bats +# -*- shell-script -*- + +# Caution: if this file has .bats extension it is fully generated. +# Any edits are subject to loss. +# Instead, edit spdir-metric.tmpl and run spdir-metric-gen-tests.sh. + +bats_load_library "wct-bats.sh" + +DETECTOR=pdsp + +# The detector variant. +VARIANT="spdir" + +function gen_fname () { + local tier="$1" ; shift + local num="$1" ; shift + local ext="${1:-npz}" + printf "%s-%s-%02.0f.%s" "$DETECTOR" "$tier" "$num" "$ext" +} + + +# The MicroBooNE conventional directions are defined in X-axis of figure 30 of +# the the MB SP paper with angle definitions in figure 29. +# +# Two big caveats in the MB convention: +# +# 1) It relies on U/V wire direction being mirror-symmetric about W. +# +# 2) It samples U/V wire direction space differently than W. Event when the +# "primed" (U/V) angle theta_xz' matches the W angle theta_xz, the theta_y' and +# theta_y differ. And so it is not trivial to make U/V vs W comparisons. +mb_theta_xz_deg=( 0 1 3 5 10 20 30 45 60 75 80 82 84 89 ) +mb_theta_y_deg=( 90 90 90 90 90 90 90 90 90 90 90 90 90 90 ) + +# In future, perhaps we pick others. For now, accept MB's +theta_xz_deg=( ${mb_theta_xz_deg[@]} ) +theta_y_deg=( ${mb_theta_y_deg[@]} ) +ndirs="${#theta_xz_deg[@]}" +directions=$(seq 0 $(( $ndirs - 1 ))) + +function generate_depos () +{ + local num="$1" ; shift + local theta_xz=${theta_xz_deg[$num]}'*deg' + local theta_y=${theta_y_deg[$num]}'*deg' + + echo "generate_depos: $num $theta_xz $theta_y" + + local fname="$(gen_fname depos $num)" + local pname="$(gen_fname params $num json)" + + # note: MB only select based on collection plane=2 + wcpy gen detlinegen \ + -d $DETECTOR --apa 0 --plane 2 --angle-coords=wire-plane \ + --theta_xz="$theta_xz" --theta_y="$theta_y" \ + -o "$fname" -m "$pname" + + # avoid empty zip of 22 bytes + file_larger_than "$fname" 23 +} + + +@test "spdir metric pdsp generate depos" { + cd_tmp file + + for num in $directions + do + generate_depos $num + done +} + + +@test "spdir metric pdsp sim+nf+sp one job graph viz" { + cd_tmp file + local cfg_file="$(relative_path spdir-metric.jsonnet)" + local num="${directions[0]}" + local depos="$(gen_fname depos "$num")" + # save out the intermediate drifted depos for input to metric + local drifts="$(gen_fname drifts "$num")" + local splats="$(gen_fname splats "$num")" + local signals="$(gen_fname signals "$num")" + local pdf="$(gen_fname graph "$num" pdf)" + local log="$(gen_fname wct "$num" log)" + wirecell-pgraph dotify \ + -A input="$depos" \ + -A output="{drift:\"$drifts\",splat:\"$splats\",sp:\"$signals\"}" \ + -A detector=$DETECTOR \ + -A variant="$VARIANT" \ + "$cfg_file" "$pdf" +} + +@test "spdir metric pdsp sim+nf+sp run job" { + cd_tmp file + local cfg_file="$(relative_path spdir-metric.jsonnet)" + for num in ${directions[*]} + do + local depos="$(gen_fname depos "$num")" + # save out the intermediate drifted depos for input to metric + local drifts="$(gen_fname drifts "$num")" + local splats="$(gen_fname splats "$num")" + local signals="$(gen_fname signals "$num")" + local log="$(gen_fname wct "$num" log)" + wire-cell -c "$cfg_file" \ + -l "$log" -L debug \ + -A input="$depos" \ + --tla-code output="{drift:\"$drifts\",splat:\"$splats\",sp:\"$signals\"}" \ + -A detector=$DETECTOR \ + -A variant="$VARIANT" + + done +} + + +# FIXME +# +# We do not yet have a command to generate a metric from each frame file. +# When we do, fix this function. For now, it does some bogosity. +# +# This function should consume params.json +function calcualte_metrics () +{ + local jname="$(gen_fname metrics $num json)" + cat < $jname +{ + "num": $num + "theta": $num, + "phi": $num, + "nchan": 0, + "sigtot": 0, + "eff": [], + "bias": [], + "res": [] +} +EOF + # This metrics.json should produce rich but summary data. Eg, the + # eff/bias/res should be arrays the hit channels. +} +# FIXME + +@test "spdir metric pdsp metrics" { + cd_tmp file + for num in ${directions[*]} + do + calcualte_metrics "$num" + done + +} + + +# FIXME: we also need a canonical plotter that consumes all metric.json. +@test "spdir metric pdsp plots" { + cd_tmp file + for num in ${directions[*]} + do + ls -l "$(gen_fname metrics $num json)" >> "$(gen_fname plots all txt)" + done +} diff --git a/test/test/test-spdir-metric-uboone.bats b/test/test/test-spdir-metric-uboone.bats new file mode 100644 index 000000000..7636173e3 --- /dev/null +++ b/test/test/test-spdir-metric-uboone.bats @@ -0,0 +1,158 @@ +#!/usr/bin/env bats +# -*- shell-script -*- + +# Caution: if this file has .bats extension it is fully generated. +# Any edits are subject to loss. +# Instead, edit spdir-metric.tmpl and run spdir-metric-gen-tests.sh. + +bats_load_library "wct-bats.sh" + +DETECTOR=uboone + +# The detector variant. +VARIANT="spdir" + +function gen_fname () { + local tier="$1" ; shift + local num="$1" ; shift + local ext="${1:-npz}" + printf "%s-%s-%02.0f.%s" "$DETECTOR" "$tier" "$num" "$ext" +} + + +# The MicroBooNE conventional directions are defined in X-axis of figure 30 of +# the the MB SP paper with angle definitions in figure 29. +# +# Two big caveats in the MB convention: +# +# 1) It relies on U/V wire direction being mirror-symmetric about W. +# +# 2) It samples U/V wire direction space differently than W. Event when the +# "primed" (U/V) angle theta_xz' matches the W angle theta_xz, the theta_y' and +# theta_y differ. And so it is not trivial to make U/V vs W comparisons. +mb_theta_xz_deg=( 0 1 3 5 10 20 30 45 60 75 80 82 84 89 ) +mb_theta_y_deg=( 90 90 90 90 90 90 90 90 90 90 90 90 90 90 ) + +# In future, perhaps we pick others. For now, accept MB's +theta_xz_deg=( ${mb_theta_xz_deg[@]} ) +theta_y_deg=( ${mb_theta_y_deg[@]} ) +ndirs="${#theta_xz_deg[@]}" +directions=$(seq 0 $(( $ndirs - 1 ))) + +function generate_depos () +{ + local num="$1" ; shift + local theta_xz=${theta_xz_deg[$num]}'*deg' + local theta_y=${theta_y_deg[$num]}'*deg' + + echo "generate_depos: $num $theta_xz $theta_y" + + local fname="$(gen_fname depos $num)" + local pname="$(gen_fname params $num json)" + + # note: MB only select based on collection plane=2 + wcpy gen detlinegen \ + -d $DETECTOR --apa 0 --plane 2 --angle-coords=wire-plane \ + --theta_xz="$theta_xz" --theta_y="$theta_y" \ + -o "$fname" -m "$pname" + + # avoid empty zip of 22 bytes + file_larger_than "$fname" 23 +} + + +@test "spdir metric uboone generate depos" { + cd_tmp file + + for num in $directions + do + generate_depos $num + done +} + + +@test "spdir metric uboone sim+nf+sp one job graph viz" { + cd_tmp file + local cfg_file="$(relative_path spdir-metric.jsonnet)" + local num="${directions[0]}" + local depos="$(gen_fname depos "$num")" + # save out the intermediate drifted depos for input to metric + local drifts="$(gen_fname drifts "$num")" + local splats="$(gen_fname splats "$num")" + local signals="$(gen_fname signals "$num")" + local pdf="$(gen_fname graph "$num" pdf)" + local log="$(gen_fname wct "$num" log)" + wirecell-pgraph dotify \ + -A input="$depos" \ + -A output="{drift:\"$drifts\",splat:\"$splats\",sp:\"$signals\"}" \ + -A detector=$DETECTOR \ + -A variant="$VARIANT" \ + "$cfg_file" "$pdf" +} + +@test "spdir metric uboone sim+nf+sp run job" { + cd_tmp file + local cfg_file="$(relative_path spdir-metric.jsonnet)" + for num in ${directions[*]} + do + local depos="$(gen_fname depos "$num")" + # save out the intermediate drifted depos for input to metric + local drifts="$(gen_fname drifts "$num")" + local splats="$(gen_fname splats "$num")" + local signals="$(gen_fname signals "$num")" + local log="$(gen_fname wct "$num" log)" + wire-cell -c "$cfg_file" \ + -l "$log" -L debug \ + -A input="$depos" \ + --tla-code output="{drift:\"$drifts\",splat:\"$splats\",sp:\"$signals\"}" \ + -A detector=$DETECTOR \ + -A variant="$VARIANT" + + done +} + + +# FIXME +# +# We do not yet have a command to generate a metric from each frame file. +# When we do, fix this function. For now, it does some bogosity. +# +# This function should consume params.json +function calcualte_metrics () +{ + local jname="$(gen_fname metrics $num json)" + cat < $jname +{ + "num": $num + "theta": $num, + "phi": $num, + "nchan": 0, + "sigtot": 0, + "eff": [], + "bias": [], + "res": [] +} +EOF + # This metrics.json should produce rich but summary data. Eg, the + # eff/bias/res should be arrays the hit channels. +} +# FIXME + +@test "spdir metric uboone metrics" { + cd_tmp file + for num in ${directions[*]} + do + calcualte_metrics "$num" + done + +} + + +# FIXME: we also need a canonical plotter that consumes all metric.json. +@test "spdir metric uboone plots" { + cd_tmp file + for num in ${directions[*]} + do + ls -l "$(gen_fname metrics $num json)" >> "$(gen_fname plots all txt)" + done +} diff --git a/test/test/test-sss-pdsp.bats b/test/test/test-sss-pdsp.bats new file mode 100644 index 000000000..7b039befd --- /dev/null +++ b/test/test/test-sss-pdsp.bats @@ -0,0 +1,74 @@ +#!/usr/bin/env bats + +bats_load_library "wct-bats.sh" + +function gen_fname () { + local tier="$1" ; shift + local num="$1" ; shift + local ext="${1:-npz}" + printf "%s-%s-%02.0f.%s" "pdsp" "$tier" "$num" "$ext" +} + +function generate_depos () { + local num="1" + local fname="$(gen_fname depos $num)" + local pname="$(gen_fname params $num json)" + + # APA0, Face 0 is approx at x=[0,-3.6m], y=[0,6m], z=[0,2.6m] + # local corner='-1*m,3*m,1.0*m' + local corner='0,0,0' + # local diagonal='-1*m,0*m,0.6*m' + local diagonal='-3.6*m,6*m,2.6*m' + + + # wcpy gen depo-lines --track-info "$pname" -T 0.0 -c "$corner" -d "$diagonal" -o "$fname" --seed $num + wcpy gen depo-point -n 50000 -p '-1*m,3*m,1.3*m' -o "$fname" + touch "$pname" + + # avoid empty zip of 22 bytes + file_larger_than "$fname" 23 +} + +@test "sss generate depos" { + cd_tmp file + generate_depos +} + + +@test "sss run job" { + cd_tmp file + + local cfg_file="$(relative_path spdir-metric.jsonnet)" + local num=1 + local depos="$(gen_fname depos "$num")" + # save out the intermediate drifted depos for input to metric + local drifts="$(gen_fname drifts "$num")" + local splats="$(gen_fname splats "$num")" + local signals="$(gen_fname signals "$num")" + local log="$(gen_fname wct "$num" log)" + + # sort of pointless until depo-lines is also run idempotently + run_idempotently -s "$depos" -t "$drifts" -t "$splats" -t "$signals" -- \ + wire-cell -c "$cfg_file" \ + -l "$log" -L debug \ + -A variant=simple \ + -A input="$depos" \ + --tla-code output="{drift:\"$drifts\",splat:\"$splats\",sp:\"$signals\"}" \ + -A detector=pdsp \ + -A tasks="drift,splat,sim,sp" + +} + +@test "sss plots" { + cd_tmp file + local num=1 + + for thing in splats signals + do + + wirecell-plot frame -t '*' -n wave \ + -o "$(gen_fname "$thing" "$num" pdf)" \ + --single "$(gen_fname "$thing" "$num")" + done + +} diff --git a/test/test/test-ssss-pdsp.bats b/test/test/test-ssss-pdsp.bats new file mode 100644 index 000000000..4d6bb4013 --- /dev/null +++ b/test/test/test-ssss-pdsp.bats @@ -0,0 +1,134 @@ +#!/usr/bin/env bats + + +bats_load_library "wct-bats.sh" + + +setup_file () { + cd_tmp + + 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 + +} + +@test "ssss pdsp depo line" { + cd_tmp file + + read -a args < depo-line-args + + run_idempotently -s depo-line-args -t depos.npz -- \ + wirecell-gen "${args[@]}" + [[ -s depos.npz ]] +} + +@test "ssss pdsp do splats" { + cd_tmp file + local cfg_file="$(relative_path spdir-metric.jsonnet)" + + for what in nominal smeared + do + run_idempotently -s "$cfg_file" -t "splat-$what.log" -t "splat-$what.npz" -- \ + wire-cell \ + -l "splat-$what.log" -L debug \ + -A input="depos.npz" \ + --tla-code output="{splat:\"splat-$what.npz\"}" \ + -A detector=pdsp \ + -A variant="ssss_$what" \ + -A tasks="drift,splat" \ + "$cfg_file" + done +} + +@test "ssss pdsp gen cfg" { + cd_tmp file + local cfg_file="$(relative_path spdir-metric.jsonnet)" + + run_idempotently -s "$cfg_file" -t dfp.json -- \ + wcsonnet \ + -A input="depos.npz" \ + --tla-code output="{drift:\"drift.npz\",sp:\"signal.npz\"}" \ + -A detector=pdsp \ + -A variant="ssss_nominal" \ + -A tasks="drift,sim,sp" \ + -o dfp.json \ + "$cfg_file" +} + +@test "ssss pdsp graph viz" { + cd_tmp file + + run_idempotently -s dfp.json -t dfp.pdf -- \ + wirecell-pgraph dotify dfp.json dfp.pdf + +} + +@test "ssss pdsp full wct" { + cd_tmp file + + run_idempotently -s dfp.json -t wct.log -t drift.npz -t splat.npz -t signal.npz -- \ + wire-cell -l wct.log -L debug dfp.json + +} + +@test "ssss pdsp log summary" { + cd_tmp file + + # {"anode":-3578.36,"cathode":-1.587500000000091,"response":-3487.8875} + jq -c '.[]|select (.type|contains("AnodePlane"))|.data.faces[0]' < dfp.json > xregions.json + + # [(-3578.36 76.1 3.3497) --> (-1.5875 6060 2303.37)] + grep sensvol wct.log | grep 'face:0'|head -1|sed 's/^.*sensvol: //g' > sensvol.txt + + # we use these numbers for the line to start with. + + # fixme: check that we get the numbers we expect +} + +@test "ssss pdsp plots" { + cd_tmp file + + + run_idempotently -s depos.npz -t depos-qxz.pdf -- \ + wirecell-gen plot-depos --single -p qxz \ + -o depos-qxz.pdf depos.npz + + run_idempotently -s drift.npz -t drift-qzt.pdf -- \ + wirecell-gen plot-depos --single -p qzt \ + -o drift-qzt.pdf drift.npz + + local helper="$(relative_path ssss-pdsp.py)" + + for what in nominal smeared + do + + for ext in png pdf + do + + run_idempotently -s "$helper" \ + -s depos.npz -s drift.npz -s "splat-$what.npz" -s signal.npz \ + -t "plots-$what.$ext" -- \ + python "$helper" plots \ + --output="plots-$what.$ext" \ + --channel-ranges '0,800,1600,2560' \ + {depos,drift,"splat-$what",signal}.npz + + done + done +} + +@test "ssss pdsp summary" { + cd_tmp file + + local helper="$(relative_path ssss-pdsp.py)" + python "$helper" log-summary wct.log > summary.org + python "$helper" depo-summary depos.npz drift.npz >> summary.org + +} diff --git a/test/wct-bats.sh b/test/wct-bats.sh index b6cd66f42..7b7650317 100755 --- a/test/wct-bats.sh +++ b/test/wct-bats.sh @@ -167,6 +167,7 @@ function check () { info "RUNNING: $*" local status run "$@" + if [ -n "$output" ] ; then info "OUTPUT:\n$output" else @@ -1091,6 +1092,29 @@ function download_file () { } +# mv_if_diff +# +# Move to if they differ. +# +# This can be useful to initate a chain of run_idempotently steps. +mv_if_diff () { + local src="$1" ; shift + local dst="$1" ; shift + + if [ ! -f "$dst" ] ; then + echo "dst does not exist $dst moving $src" 1>&3 + mv "$src" "$dst" + return + fi + if [ -n "$( diff "$src" "$dst" )" ] ; then + echo "differ: $dst and $src" 1>&3 + mv "$src" "$dst" + return + fi + rm -f "$src" +} + + function help_one () { local func="$1" ; shift diff --git a/util/inc/WireCellUtil/TagRules.h b/util/inc/WireCellUtil/TagRules.h index 53cbcdf8e..4742a6ab4 100644 --- a/util/inc/WireCellUtil/TagRules.h +++ b/util/inc/WireCellUtil/TagRules.h @@ -68,7 +68,15 @@ namespace WireCell { class Context { std::unordered_map > m_rulesets; - public: + public: + + // Return the number of rulesets under a given rule name. + size_t nrules(const std::string& name) const { + auto it = m_rulesets.find(name); + if (it == m_rulesets.end()) return 0; + return it->second.size(); + } + // This should be an array of objects each keyed by a // context name and with values providing a ruleset. void configure(const Configuration& cfg); diff --git a/wct-deps.png b/wct-deps.png deleted file mode 100644 index cf9c0657b..000000000 Binary files a/wct-deps.png and /dev/null differ