More plots checking LMN

wirecell/resp/__main__.py

@@ -596,6 +596,103 @@ def newpage(fig, name, title=''):
             oo.write('\n')
 
 
+@cli.command("lmn-pdsp-plots")
+@click.argument("output")
+def lmn_pdsp_plots(output):
+    '''
+    Generate PDF file with plots illustrating LMN on PDSP 
+    '''
+    from wirecell.util import lmn
+    import wirecell.resp.resample as res
+    from wirecell.resp.plots import load_fr, eresp, plot_paths, multiply_period
+    from wirecell.resp.util import fr2sigs
+    from wirecell.util.plottools import pages
+
+    FRs = [load_fr("pdsp")]
+    FRs[0].period = 100*units.ns  # fix inaccuracies
+    FRs.append(res.resample(FRs[0], 64*units.ns))
+
+    names = ["top", "bot"]
+    ss_fast = [lmn.Sampling(T=fr.period,
+                            N=fr.planes[0].paths[0].current.size)
+                 for fr in FRs]
+    er_fast = [eresp(ss, name="ER_{%s}" % name)
+               for ss, name in zip(ss_fast, names)]
+    fr_fast = [fr2sigs(fr) for fr in FRs]
+
+    def visit_lol(lol, func):
+        if isinstance(lol, list):
+            return [visit_lol(one, func) for one in lol]
+        return func(lol)
+
+    qr_fast = [visit_lol(frs,
+                         lambda sig: multiply_period(lmn.convolve(sig, er)))
+               for er, frs in zip(er_fast, fr_fast)]
+
+    downsample_factors = (5, 8)
+    qr_slow = [visit_lol(one, lambda sig:
+                         lmn.interpolate(sig, sig.sampling.T*dsf))
+               for one, dsf in zip(qr_fast, downsample_factors)]
+
+    Tslow = qr_slow[0][0][0].sampling.T
+    ss_slow = [qr[0][0].sampling for qr in qr_slow]
+    qr_slow_rs = visit_lol(qr_slow[1], lambda sig:
+                           lmn.interpolate(sig, Tslow))
+    qr_slow.append(qr_slow_rs)
+    names.append('res')
+
+
+    # ss_slow = [dr[0][0].sampling for dr in dr_slow]
+    # for one in ss_slow:
+    #     print('slow',one)
+
+    # dr_slow_rs = visit_lol(dr_slow[1], lambda sig:
+    #                        lmn.interpolate(sig, ss_slow[0].T))
+    # dr_slow.append(dr_slow_rs)
+    # names.append('res')
+
+    # qr_slow = visit_lol(dr_slow, lambda sig: multiply_period(sig))
+
+    # qr_diff = [
+    #     visit_lol(list(zip(pa,pb)), lambda sigs:
+    #               lmn.Signal(sigs[0].sampling,
+    #                          wave=sigs[0].wave-sigs[1].wave[:ss_slow[0].N]))
+    #            for pa, pb in zip(qr_slow[0], qr_slow[-1])]
+
+    with pages(output) as out:
+        def page(name):
+            print(name)
+            plt.suptitle(name)
+            plt.tight_layout()
+            out.savefig()
+            plt.close()
+
+        for iplane, letter in enumerate("UVW"):
+
+            def trio(kind, responses):
+                indices = [0,1]
+                if len(responses) == 3:
+                    indices = [2,0,1]
+                for ind in indices:
+                    name = names[ind]
+                    rs = responses[ind]
+                    rplane = rs[iplane]
+                    plot_paths(rplane)
+                    page(f'{letter} {name} ${kind}$ {rplane[0].sampling}')
+
+
+
+            trio('FR', fr_fast)
+            trio('T \cdot FR \circledast ER', qr_fast)
+            trio('T \cdot FR \circledast ER', qr_slow)
+            # trio('FRxER', dr_slow)
+            # trio('T.FRxER', qr_slow)
+
+            # plot_paths(qr_diff[iplane])
+            # page(f'T.FRxER diff {letter}')
+
+
+
 def main():
     cli(obj=dict())
 

wirecell/resp/plots.py

@@ -56,13 +56,14 @@ def spectrum_resize(spec, newsize, norm):
     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 = numpy.array(sig.wave)
     wave.resize( int(short_length * tick / rawresp_tick ) )
 
     # unlike WCT PIR, we keep this an interpolation instead of directly jumping
@@ -72,10 +73,10 @@ def wct_pir_resample(sig, tick, short_length, name=None):
     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}')
+                     name=name or sig.name)
+    # print(f'{sig=} {ret=}')
     return ret
-    
+
 
 
 
@@ -88,13 +89,15 @@ def eresp(ss, name="coldelec", gain=14*units.mV/units.fC, shaping=2*units.us):
     eresp = electronics(times, gain, shaping)
     return lmn.Signal(ss, wave = eresp, name=name)
 
+
 def label(sig):
     '''
     Return a legend label for signal.
     '''
     ss = sig.sampling
     return f'${sig.name},\ N={ss.N},\ T={ss.T/units.ns:.0f}\ ns$'
 
+
 def plot_signals(sigs, tlim=None, flim=None, tunits="us", funits="MHz",
                  iunits="femtoampere", ilim=None,
                  drawstyle='steps-mid', linewidth = 1,
@@ -322,3 +325,48 @@ def load_fr(detector):
     if isinstance(got, list):
         return got[0]
     return got
+
+
+def plot_paths(rplane,
+               tlim=(55*units.us, 70*units.us), flim=(0, 1*units.MHz),
+               pixels=True):
+    '''
+    Plot path response signals in one plane
+    '''
+    tustr = "us"
+    fustr = "MHz"
+    tuval = getattr(units, tustr)
+    fuval = getattr(units, fustr)
+
+    fig, axes = plt.subplots(ncols=2, nrows=2, sharex='col')
+
+    npaths = len(rplane)
+    imps = numpy.linspace(0, npaths, npaths, endpoint=False)
+    ss = rplane[0].sampling
+    t_lrtb = (0, ss.duration/tuval, 0, npaths)
+    f_lrtb = (0, ss.F/fuval, 0, npaths)
+
+    wave = lmn.sigs2arr(rplane, "wave")
+    spec = numpy.abs(lmn.sigs2arr(rplane, "spec"))
+
+    tix, fix = axes[0]
+    # pcolormesh is slower than imshow
+    if pixels:
+        kwds = dict(aspect='auto', interpolation='none')
+        tix.imshow(wave, extent=t_lrtb, cmap="seismic", **kwds)
+        fix.imshow(spec, extent=f_lrtb, **kwds)
+    else:
+        kwds = dict()
+        tt_mg, ti_mg = numpy.meshgrid(ss.times/tuval, imps)
+        ff_mg, fi_mg = numpy.meshgrid(ss.freqs/fuval, imps)
+        tix.pcolormesh(tt_mg, ti_mg, wave, cmap="seismic", **kwds)
+        fix.pcolormesh(ff_mg, fi_mg, spec, **kwds)
+    tix.set_xlim(tlim[0]/tuval, tlim[1]/tuval)
+    fix.set_xlim(flim[0]/fuval, flim[1]/fuval)
+
+    tpx, fpx = axes[1]
+    for irow, sig in enumerate(rplane):
+        tpx.plot(ss.times/tuval, sig.wave)
+        fpx.plot(ss.freqs/fuval, numpy.abs(sig.spec))
+
+    return fig, axes

wirecell/resp/resample.py

@@ -43,6 +43,7 @@ def resample(fr, Tr, Ts=None, **kwds):
                                     plane.pitch))
     return FieldResponse(planes, fr.axis, fr.origin, fr.tstart, Tr, fr.speed)
 
+
 def rolloff(fr, fstart=0.9):
     '''
     Return an FR with a linear roll-off applied.

wirecell/resp/util.py

@@ -0,0 +1,20 @@
+from wirecell.sigproc.response.arrays import pr2array
+from wirecell.util import lmn
+
+def pr2sigs(pr, period):
+    '''
+    Return list of signals for path responses in plane response pr.
+    '''
+    return [
+        lmn.Signal(lmn.Sampling(T=period, N=imp.current.size),
+                   wave=imp.current, name=str(imp.pitchpos))
+        for imp in pr.paths
+    ]
+
+
+def fr2sigs(fr):
+    '''
+    Return list of list of response signals for each plane.
+    '''
+    return [pr2sigs(pr, fr.period) for pr in fr.planes]
+

wirecell/sigproc/response/arrays.py

@@ -41,18 +41,18 @@ def pr2array(pr, nimperwire = 6, nbinsperwire = 10):
             p1 = pr.paths[ibin0 + ind+0].pitchpos
             p2 = pr.paths[ibin0 + ind+1].pitchpos
             pm = 0.5*(p1+p2)
-            
-            obin = iwire * nbinsperwire + ind;
+
+            obin = iwire * nbinsperwire + ind
 
             res[obin] = m
             pitches[obin] = pm
 
     res = res + numpy.flipud(res)
     pitches = pitches - numpy.flip(pitches)
-        
+
     # for path in pr.paths:
     #     print ("%.3f mm"%(path.pitchpos/units.mm))
-    return res,pitches
+    return res, pitches
 
 
 def toarray(fr, nimperwire=6, nbinsperwire=10):
@@ -162,7 +162,7 @@ def coldelec(fra, fr, gain, shaping):
     # for ind, pr in enumerate(fr.planes):
     #     fra['resp%d' % pr.planeid] = responses[ind]
 
-    return fra;
+    return fra
 
 
 def fr2arrays(fr, gain=None, shaping=None):

wirecell/util/lmn.py

@@ -247,6 +247,13 @@ def interval_multiply(self, other, name=None):
                       name=name or self.name)
 
 
+def sigs2arr(sigs, which='wave'):
+    '''
+    Given list of consistent sigs, return 2D array of wave or spec
+    '''
+    return numpy.vstack([getattr(sig, which) for sig in sigs])
+
+
 def bezout(a, b, eps=1e-6):
     '''Greated common divisor and Bezout coefficients.
 
@@ -558,3 +565,4 @@ def convolve_downsample(sa, sb, name=None):
     sced = Signal(sbe.sampling, spec=saed.spec * sbe.spec,
                   name=name or f'{sa.name} (x) {sb.name}')
     return sced
+

wirecell/util/test/test_cont_conv.py

@@ -0,0 +1,4 @@
+#!/usr/bin/env pytest
+'''
+Test method for continuous convolution.
+'''