diff --git a/wirecell/resp/__main__.py b/wirecell/resp/__main__.py index 4ae2910..7a7a29f 100644 --- a/wirecell/resp/__main__.py +++ b/wirecell/resp/__main__.py @@ -323,6 +323,10 @@ 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.") @@ -330,7 +334,7 @@ def plot_paths(axes, rfile, pind, n, plane, trange, frange): 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. ''' @@ -338,7 +342,9 @@ 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, - 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,7 +358,8 @@ 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) @@ -360,13 +367,18 @@ def lmn_fr_plots(impact, plane, period, 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) diff --git a/wirecell/util/lmn.py b/wirecell/util/lmn.py index 870750b..b8fe995 100644 --- a/wirecell/util/lmn.py +++ b/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