Miniscule fixes for LMN related

wirecell/gen/noise.py

@@ -8,26 +8,27 @@
 
 import numpy
 from math import sqrt, pi, ceil, floor
+from wirecell.lmn import hermitian_mirror
 
 def rayleigh(x,sigma=1):
     s2 = sigma**2
     return (x/s2)*numpy.exp(-0.5*x**2/s2)
 
-def hermitian_mirror(spec):
-    hm = numpy.array(spec)
-
-    # nyquist bin index if size is even, else bin just below Nyquist edge.
-    halfsize = hm.size//2
-
-    # zero freq must be real
-    hm[0] = numpy.abs(hm[0])
-    hm[1:halfsize] = hm[1:halfsize]
-    if 0 == hm.size%2:          # even with Nyquist bin
-        hm[halfsize] = numpy.abs(hm[halfsize])
-        hm[halfsize+1:] = numpy.conjugate(hm[halfsize-1:0:-1])
-    else:
-        hm[halfsize+1:] = numpy.conjugate(hm[halfsize:0:-1])
-    return hm
+# def hermitian_mirror(spec):
+#     hm = numpy.array(spec)
+
+#     # nyquist bin index if size is even, else bin just below Nyquist edge.
+#     halfsize = hm.size//2
+
+#     # zero freq must be real
+#     hm[0] = numpy.abs(hm[0])
+#     hm[1:halfsize] = hm[1:halfsize]
+#     if 0 == hm.size%2:          # even with Nyquist bin
+#         hm[halfsize] = numpy.abs(hm[halfsize])
+#         hm[halfsize+1:] = numpy.conjugate(hm[halfsize-1:0:-1])
+#     else:
+#         hm[halfsize+1:] = numpy.conjugate(hm[halfsize:0:-1])
+#     return hm
 
 def fictional(freqs, rel=0.1):
     r = rayleigh(freqs, freqs[-1]*rel)

wirecell/plot/__main__.py

@@ -345,9 +345,11 @@ def plot_slice(ax, slc):
               help="limit y-axis range in raw numbers, default is auto range")
 @click.option("-f", "--frange", default=None, type=str,
               help="limit frequency range, eg '0,100*kHz'")
+@click.option("--drawstyle", default="steps-mid", type=str,
+              help="draw style to use")
 @image_output
 @click.argument("frame_files", nargs=-1)
-def channels(output, channel, trange, frange, yrange, frame_files, **kwds):
+def channels(output, channel, trange, frange, yrange, drawstyle, frame_files, **kwds):
     '''
     Plot channels from multiple frame files.
 
@@ -378,8 +380,13 @@ def channels(output, channel, trange, frange, yrange, frame_files, **kwds):
 
     frames = {ff: list(load_frames(ff)) for ff in frame_files}
 
-    drawstyle = 'steps-mid'
-
+    progressive = False
+    if drawstyle == "progressive":
+        progressive = True
+        drawstyle = None
+    if not drawstyle:
+        drawstyle = None
+
     with output as out:
 
         for chan in channels:
@@ -389,19 +396,19 @@ def channels(output, channel, trange, frange, yrange, frame_files, **kwds):
 
             for ind, (fname, frs) in enumerate(frames.items()):
                 stem = Path(fname).stem
-                # linewidth = len(frames) - ind
-                linewidth = 1
+                if progressive:
+                    linewidth = len(frames) - ind
+                else:
+                    linewidth = 1
 
                 for fr in frs:
 
                     wave = fr.waves(chan)
                     axes[0].set_title("waveforms")
-                    print(f'{linewidth=}')
                     axes[0].plot(fr.times/units.us, wave, linewidth=linewidth, drawstyle=drawstyle)
                     if trange:
                         axes[0].set_xlim(trange[0]/units.us, trange[1]/units.us)
                     if yrange:
-                        print(yrange)
                         axes[0].set_ylim(*yrange)
                     else:
                         plottools.rescaley(axes[0], fr.times/units.us, wave,
@@ -419,7 +426,6 @@ def channels(output, channel, trange, frange, yrange, frame_files, **kwds):
                         axes[1].set_xlim(0, fr.freqs_MHz[fr.nticks//2])
                     axes[1].set_xlabel("frequency [MHz]")
                     axes[1].legend(loc='upper right')
-                    print(fr.nticks, fr.period/units.ns, fr.duration/units.us)
 
 
             if not out.single:

wirecell/resp/__main__.py

@@ -325,7 +325,7 @@ def plot_paths(axes, rfile, pind, n, plane, trange, frange):
               "with units, eg '64*ns'")
 @click.option("--slow-tick", default='500*ns',
               help="Sample period of ADC for nominal ER")
-@click.option("--slow-nticks", default=200,
+@click.option("--slow-nticks", default=200, type=int,
               help="Number of ADC ticks for nominal ER")
 @click.option("-O", "--org-output", default=None,
               help="Generate an org mode file with macros expanding to "
@@ -341,7 +341,7 @@ def lmn_fr_plots(impact, plane, period,
 
     from wirecell.util import lmn
     from wirecell.resp.plots import (
-        load_fr, multiply_period, fr_sig, fr_arr, eresp,
+        load_fr, multiply_period, fr_sig, fr_arr, eresp, wct_pir_resample,
         plot_signals, plot_ends, plot_wave_diffs, plot_shift)
 
     slow_tick = unitify(slow_tick)
@@ -358,27 +358,32 @@ def lmn_fr_plots(impact, plane, period,
         FRs.period = unitify(period)
     nrat = lmn.rational_size(FRs.period, Tr)
     sigs_orig = fr_sig(FRs, "I_{og}", impact, plane)
-    print(f'Ns_orig={sigs_orig.sampling.N} '
-          f'Ts={sigs_orig.sampling.T} {Tr=} -> {nrat=}')
+    # print(f'Ns_orig={sigs_orig.sampling.N} '
+    #       f'Ts={sigs_orig.sampling.T} {Tr=} -> {nrat=}')
+
+    sig_pir = wct_pir_resample(sigs_orig, slow_tick, slow_nticks,
+                               name='I_{pir}')
 
     sigs = lmn.rational(sigs_orig, Tr)
     arrs = fr_arr(FRs)
 
     FRr = res.resample(FRs, Tr)
     sigr = fr_sig(FRr, "I_{rs}", impact, plane)
     arrr = fr_arr(FRr)
-    print(f'Ns={sigs.sampling.N} Ts={sigs.sampling.T} '
-          f'Nr={sigr.sampling.N} Tr={sigr.sampling.T}')
+    # print(f'Ns={sigs.sampling.N} Ts={sigs.sampling.T} '
+    #       f'Nr={sigr.sampling.N} Tr={sigr.sampling.T}')
 
     ers = eresp(sigs.sampling, "E_{og}")
     err = eresp(sigr.sampling, "E_{rs}")
 
     # check simultaneous downsample+convolution
     slow_sampling = lmn.Sampling(T=slow_tick, N=slow_nticks)
     ers_slow = eresp(slow_sampling, "E_{og,slow}")
-    vrs_slow = lmn.convolve_downsample(sigs, ers_slow)
-    # should also work when samplings match
-    vrs_fast = lmn.convolve_downsample(sigs, ers)
+
+    # check voltage.  Needs tick multiplied to FRxER
+    # vrs_slow = lmn.convolve_downsample(sigs, ers_slow, name='V_{cd}')
+    vrs_pir  = multiply_period(lmn.convolve(sig_pir, ers_slow), name="V_{pir}")
+    vrs_fast = multiply_period(lmn.convolve(sigs, ers), name="V_{fast}")
 
     dsigs = lmn.convolve(sigs, ers, name=r'I_{og} \otimes E_{og}')
     dsigr = lmn.convolve(sigr, err, name=r'I_{rs} \otimes E_{rs}')
@@ -409,7 +414,8 @@ def lmn_fr_plots(impact, plane, period,
     # convolve then downsample/decimate
     qdsigs_cds = lmn.interpolate(dtsigs, sim_tick,name='V_{cds}')
     qdsigs_cdm = lmn.decimate(dtsigs, decimation, name='V_{cdm}')
-    print(f'{dtsigs.wave.shape=} {qdsigs_cds.wave.shape=} {qdsigs_dsc.wave.shape=} {qdsigs_cdm.wave.shape=} {qdsigs_dmc.wave.shape=}')
+    # print(f'{dtsigs.wave.shape=} {qdsigs_cds.wave.shape=} '
+    #       f'{qdsigs_dsc.wave.shape=} {qdsigs_cdm.wave.shape=} {qdsigs_dmc.wave.shape=}')
 
     ## end building data arrays
 
@@ -453,33 +459,33 @@ def newpage(fig, name, title=''):
         frtit = 'current field response: $FR$ ' \
             f'({detector_name} plane {plane} impact {impact})'
 
-        fig, _ = plot_signals((sigs, sigr), iunits='femtoampere', **full_range)
-        newpage(fig, 'fig-fr', frtit)
+        # fig, _ = plot_signals((sigs, sigr), iunits='femtoampere', **full_range)
+        # newpage(fig, 'fig-fr', frtit)
 
-        fig, _ = plot_signals((sigs, sigr), iunits='femtoampere', **zoom_kwds)
+        fig, _ = plot_signals((sigs, sigr, sig_pir), iunits='femtoampere', **zoom_kwds)
         newpage(fig, 'fig-fr-zoom', frtit + ' zoom')
 
-        ilim = 0.05*units.femtoampere
-        fig, _ = plot_signals((sigs, sigr), iunits='femtoampere',
-                              ilim=(-ilim, ilim), **front_range)
-        newpage(fig, 'fig-fr-front', frtit)
+        # ilim = 0.05*units.femtoampere
+        # fig, _ = plot_signals((sigs, sigr), iunits='femtoampere',
+        #                       ilim=(-ilim, ilim), **front_range)
+        # newpage(fig, 'fig-fr-front', frtit)
 
-        fig, _ = plot_ends((arrs, arrr), (sigs.name, sigr.name),
-                           iunits='femtoampere')
-        newpage(fig, 'fig-ends', 'current response ends')
+        # fig, _ = plot_ends((arrs, arrr), (sigs.name, sigr.name),
+        #                    iunits='femtoampere')
+        # newpage(fig, 'fig-ends', 'current response ends')
 
-        fig, _ = plot_signals((ers, err),
+        fig, _ = plot_signals((ers, err, ers_slow),
                               iunits='mV/fC',
                               tlim=(0*units.us, 20*units.us),
                               flim=(0*units.MHz, 1*units.MHz))
         newpage(fig, 'fig-cer', 'cold electronics response: $ER$')
 
-        fig, _ = plot_signals((vrs_fast, vrs_slow),
-                              linewidth='progressive',
-                              iunits='mV',
-                              tlim=(40*units.us, 80*units.us),
-                              flim=(0*units.MHz, 1*units.MHz))
-        newpage(fig, 'fig-vr', 'voltage response scd')
+        # fig, _ = plot_signals((vrs_fast, vrs_pir),
+        #                       linewidth='progressive',
+        #                       iunits='mV',
+        #                       tlim=(40*units.us, 80*units.us),
+        #                       flim=(0*units.MHz, 1*units.MHz))
+        # newpage(fig, 'fig-vr', 'voltage response scd')
 
         fig, _ = plot_signals((dsigs, dsigr, dsigr2),
                               iunits='femtoampere*mV/fC',
@@ -511,8 +517,8 @@ def newpage(fig, name, title=''):
         axes[0].set_ylim(-tiny, tiny)
         newpage(fig, 'fig-v-back', f'voltage response ({detector_name} '
                 f'plane {plane} impact {impact}, zoom back)')
-        print('q=', 100*numpy.sum(qdsigs.wave) /
-              numpy.sum(numpy.abs(qdsigs.wave)), '%')
+        # print('q=', 100*numpy.sum(qdsigs.wave) /
+        #       numpy.sum(numpy.abs(qdsigs.wave)), '%')
 
         # ER in fast and slow binning
         fig,_ = plot_signals((ers, ers_ss), iunits='mV/fC',

wirecell/resp/plots.py

@@ -1,4 +1,5 @@
 import numpy
+from math import floor
 import matplotlib.pyplot as plt
 from wirecell import units
 from wirecell.util import lmn
@@ -39,6 +40,45 @@ def fr_arr(fr, plane=0):
     return ret
 
 
+def spectrum_resize(spec, newsize, norm):
+    '''
+    Replicate the function of this name from PIR.
+    '''
+    oldsize = spec.size
+    if oldsize == newsize:
+        return spec
+    oldhalf = int(1 + floor(oldsize / 2))
+    newhalf = int(1 + floor(newsize / 2))
+
+    half = min(oldhalf, newhalf)
+    ret = numpy.zeros(newsize, dtype=spec.dtype)
+    ret[:half] = spec[:half]
+    ret = norm * lmn.hermitian_mirror(ret);
+    return ret
+
+def wct_pir_resample(sig, tick, short_length, name=None):
+    '''
+    Replicate the resampling done in PIR.
+    '''
+    rawresp_tick = sig.sampling.T
+    rawresp_ntick = sig.sampling.N
+    wave = numpy.array(sig.wave);
+    wave.resize( int(short_length * tick / rawresp_tick ) )
+
+    # unlike WCT PIR, we keep this an interpolation instead of directly jumping
+    # from FR to T*FR.  And, note the norm uses original size as the
+    # wave.resize() adds no power.
+    norm = short_length/rawresp_ntick
+    spec = spectrum_resize(numpy.fft.fft(wave), short_length, norm)
+
+    ret = lmn.Signal(lmn.Sampling(T=tick, N=short_length), spec=spec,
+                     name=name or sig.name);
+    print(f'{sig} {ret}')
+    return ret
+
+
+
+
 from wirecell.sigproc.response import electronics
 def eresp(ss, name="coldelec", gain=14*units.mV/units.fC, shaping=2*units.us):
     '''

wirecell/util/__main__.py

@@ -21,6 +21,37 @@ def cli(ctx):
     pass
 
 
+@cli.command("convdown")
+@click.option("-n", "--sampling-number", default=None, type=int,
+              help="Original number of samples.")
+@click.option("-t", "--sampling-period", default=None, type=str,
+              help="Original sample period,eg '100*ns'")
+@click.option("-N", "--resampling-number", default=None, type=int,
+              help="Resampled number of samples.")
+@click.option("-T", "--resampling-period", default=None, type=str,
+              help="Resampled sample period, eg '64*ns'")
+@click.option("-e", "--error", default=1e-6,
+              help="Precision by which integer and "
+              "rationality conditions are judged")
+def cmd_convdown(sampling_number, sampling_period, resampling_period, resampling_number, error):
+    '''
+    Calculate numbers for "simultaneous convolution and downsample" method.
+
+    Eg:
+
+        $ wirecell-util convdown -n 625 -t '100*ns' -N 200 -T "500*ns"
+        (Ta=100.0 ns, Na=625) -> 1625, (Tb=500.0 ns, Nb=200) -> 325
+
+    '''
+    from wirecell.util import lmn
+    Ta = unitify(sampling_period)
+    Na = sampling_number
+    Tb = unitify(resampling_period)
+    Nb = resampling_number
+    Nea, Neb = lmn.convolution_downsample_size(Ta, Na, Tb, Nb)
+    print(f'(Ta={Ta/units.ns:.1f} ns, {Na=}) -> {Nea}, (Tb={Tb/units.ns:.1f} ns, {Nb=}) -> {Neb}')
+
+
 @cli.command("lmn")
 @click.option("-n", "--sampling-number", default=None, type=int,
               help="Original number of samples.")