Simultaneous downsample and convolution checks

wirecell/resp/__main__.py

@@ -323,22 +323,28 @@ def plot_paths(axes, rfile, pind, n, plane, trange, frange):
 @click.option("-t", "--tick", default='64*ns',
               help="Resample the field response to have this sample period "
               "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,
+              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 "
               "figure inclusion.")
 @click.option("-o", "--output", default="lmn-fr-plots.pdf")
 def lmn_fr_plots(impact, plane, period,
                  zoom_start, zoom_window, vmm,
                  detector_name, field_file,
-                 tick, org_output, output):
+                 tick, slow_tick, slow_nticks, org_output, output):
     '''
     Make plots for LMN FR presentation.
     '''
 
     from wirecell.util import lmn
     from wirecell.resp.plots import (
         load_fr, multiply_period, fr_sig, fr_arr, eresp,
-        plot_signals, plot_ends, plot_wave_diffs)
+        plot_signals, plot_ends, plot_wave_diffs, plot_shift)
+
+    slow_tick = unitify(slow_tick)
 
     Tr = unitify(tick)
     if vmm:
@@ -352,21 +358,27 @@ 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} Ts={sigs_orig.sampling.T} {Tr=} -> {nrat=}')
+    print(f'Ns_orig={sigs_orig.sampling.N} '
+          f'Ts={sigs_orig.sampling.T} {Tr=} -> {nrat=}')
 
     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} 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}")
 
-    vrs = lmn.convolve_downsample(sigs, ers)
-    vrr = lmn.convolve_downsample(sigr, err)
+    # 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)
 
     dsigs = lmn.convolve(sigs, ers, name=r'I_{og} \otimes E_{og}')
     dsigr = lmn.convolve(sigr, err, name=r'I_{rs} \otimes E_{rs}')
@@ -462,9 +474,10 @@ def newpage(fig, name, title=''):
                               flim=(0*units.MHz, 1*units.MHz))
         newpage(fig, 'fig-cer', 'cold electronics response: $ER$')
 
-        fig, _ = plot_signals((vrs, vrr, dsigs, dsigr),
+        fig, _ = plot_signals((vrs_fast, vrs_slow),
+                              linewidth='progressive',
                               iunits='mV',
-                              tlim=(0*units.us, 100*units.us),
+                              tlim=(40*units.us, 80*units.us),
                               flim=(0*units.MHz, 1*units.MHz))
         newpage(fig, 'fig-vr', 'voltage response scd')
 
@@ -502,11 +515,10 @@ def newpage(fig, name, title=''):
               numpy.sum(numpy.abs(qdsigs.wave)), '%')
 
         # ER in fast and slow binning
-        fig, _ = plot_signals((ers, ers_ds), iunits='mV/fC',
-                              tlim=(0*units.us, 20*units.us),
-                              flim=(0*units.MHz, 1*units.MHz))
-        newpage(fig, 'fig-cer-ds',
-                'slow sampled cold electronics response: $ER$')
+        fig,_ = plot_signals((ers, ers_ss), iunits='mV/fC',
+                             tlim=(0*units.us, 20*units.us),
+                             flim=(0*units.MHz, 1*units.MHz))
+        newpage(fig, 'fig-cer-ds', 'slow sampled cold electronics response: $ER$')
 
         # FR in fast and slow
         fig,_ = plot_signals((sigs, sigs_ds, sigs_dm), iunits='femtoampere', **zoom_kwds)

wirecell/util/lmn.py

@@ -84,6 +84,7 @@ def resampling(self, N):
     def __str__(self):
         return f'N={self.N} T={self.T}'
 
+
 @dataclasses.dataclass
 class Signal:
     '''
@@ -129,7 +130,10 @@ def __init__(self, sampling, wave=None, spec=None, name=None):
             self.spec = numpy.fft.fft(self.wave)
 
     def __str__(self):
-        return f'{self.name} {self.sampling}'
+        return f'"{self.name}" {self.sampling}'
+
+    def __repr__(self):
+        return f'{self.sampling} "{self.name}"'
 
     @property
     def time_energy(self):
@@ -279,7 +283,6 @@ def rational_size(Ts, Tr, eps=1e-6):
     dT = Ts - Tr
     n = dT/egcd(Tr, dT, eps=eps)
     rn = round(n)
-    print(f"delta-n = {n}")
 
     err = abs(n - rn)
     if err > eps:
@@ -366,12 +369,12 @@ def condition(signal, Tr, eps=1e-6, name=None):
 
 def resample(signal, Nr, name=None):
     '''
-    Return a new signal of same duration that is resampled to have number of samples Nr.
+    Return a new signal of same duration that is resampled
+    to have number of samples Nr.
 
     '''
     Ns = signal.sampling.N
     resampling = signal.sampling.resampling(Nr)
-    Tr = resampling.T
 
     S = signal.spec
     R = numpy.zeros(Nr, dtype=S.dtype)
@@ -484,19 +487,30 @@ def convolution_downsample_size(Ta, Na, Tb, Nb):
     return int(Nea), int(Neb)
 
 
-def convolve_downsample(sa, sb):
+def convolve_downsample(sa, sb, name=None):
     '''
     Convolve the two signals and downsample to the slower
     '''
     if sa.sampling.T > sb.sampling.T:
         sa, sb = sb, sa
-    duration = sa.sampling.duration + sb.sampling.duration
+
+    aa = sa.sampling            # eg T=100ns
+    bb = sb.sampling            # eg T=500ns
+
+    # "total" duration must be evenly divisible by both T's!
+    duration = math.ceil(aa.duration/bb.T)*bb.T + bb.duration
     sae = resize(sa, duration)
     sbe = resize(sb, duration)
-
+    print(f'{duration=} = {sa.sampling.duration} + {sb.sampling.duration}')
+    print(f'{sa.sampling.duration/sb.sampling.T} '
+          f'{sb.sampling.duration/sa.sampling.T}')
+    print(f'{sa=}\n{sb=}\n{sae=}\n{sbe=}')
     R = sa.sampling.T / sb.sampling.T
-
-    Nr = int(R*sae.sampling.N)
+    Nrf = R*sae.sampling.N
+    Nr = int(Nrf)
+    print(f'{R=} {Nr=} {Nrf=}')
     saed = resample(sae, Nr)
-    sced = Signal(sbe.sampling, spec = saed.spec * sbe.spec)
+    print(f'{saed=}')
+    sced = Signal(sbe.sampling, spec=saed.spec * sbe.spec,
+                  name=name or f'{sa.name} (x) {sb.name}')
     return sced