diff --git a/wirecell/gen/noise.py b/wirecell/gen/noise.py index ffab0e8..e6b490f 100644 --- a/wirecell/gen/noise.py +++ b/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) diff --git a/wirecell/plot/__main__.py b/wirecell/plot/__main__.py index 5071c19..21717ca 100644 --- a/wirecell/plot/__main__.py +++ b/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: diff --git a/wirecell/resp/__main__.py b/wirecell/resp/__main__.py index 7a7a29f..17a77e8 100644 --- a/wirecell/resp/__main__.py +++ b/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,8 +358,11 @@ 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) @@ -367,8 +370,8 @@ 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} ' - 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}") @@ -376,9 +379,11 @@ def lmn_fr_plots(impact, plane, period, # 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', diff --git a/wirecell/resp/plots.py b/wirecell/resp/plots.py index 35472f5..31a5cd4 100644 --- a/wirecell/resp/plots.py +++ b/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): ''' diff --git a/wirecell/util/__main__.py b/wirecell/util/__main__.py index 2642ec4..5734c98 100644 --- a/wirecell/util/__main__.py +++ b/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.") diff --git a/wirecell/util/lmn.py b/wirecell/util/lmn.py index b8fe995..f8ad392 100644 --- a/wirecell/util/lmn.py +++ b/wirecell/util/lmn.py @@ -12,6 +12,28 @@ import matplotlib.pyplot as plt from wirecell.util.cli import debug +def hermitian_mirror(spec): + ''' + Return a Hermitian-symmetric version of spec. + + The spec should be full size and the first half is Hermitian-reflected to + the second half, respecting the zero and Nyquist bin (if exists). + ''' + 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 + @dataclasses.dataclass class Sampling: ''' @@ -205,6 +227,24 @@ def resize(self, Ns, time_padding='linear', name=None): sig = Signal(ss, wave=cur, name=name or self.name) return sig + def frequency_multiply(self, other, name=None): + ''' + Multiplication in frequency space. + ''' + if isinstance(other, Signal): + other = other.spec + return Signal(self.sampling, spec = self.spec * other, + name=name or self.name) + + def interval_multiply(self, other, name=None): + ''' + Multiplication in interval space. + ''' + if isinstance(other, Signal): + other = other.wave + return Signal(self.sampling, wave = self.wave * other, + name=name or self.name) + def bezout(a, b, eps=1e-6): '''Greated common divisor and Bezout coefficients. @@ -220,9 +260,9 @@ def step(a, b): if a < eps: return (b, 0, 1) else: - print(f'{a=} {b=}') + # print(f'{a=} {b=}') g, x, y = step(b % a, a) - print(f'{g=} {a=} {b=} {x=} {y=}') + # print(f'{g=} {a=} {b=} {x=} {y=}') return (g, y - (b // a) * x, x) return step(a,b) @@ -457,15 +497,18 @@ def interpolate(sig, Tr, time_padding="zero", eps=1e-6, name=None): def convolve(s1, s2, mode='full', name=None): ''' - Return new signal that is the convolution of the two. - - Mode is as taken by numpy.convolve which provides the kernel. + Return new signal that is the linear convolution of the s1 and s2. + Linear convolution is assured by extending both signals to the sum of their + individual size by padding with zeros. Note, this is one more than the + absolute minimum required. ''' if s1.sampling.T != s2.sampling.T: raise ValueError("can not convolve signals if differing sample period") + wave = numpy.convolve(s1.wave, s2.wave, mode) - return Signal(Sampling(T=s1.sampling.T, N=wave.size), wave=wave, name=name) + return Signal(Sampling(T=s1.sampling.T, N=wave.size), wave=wave, + name=name or f'{s1.name} \otimes {s2.name}') # Not exactly lmn, but in DepoTransform we must do a convolution and a @@ -501,16 +544,16 @@ def convolve_downsample(sa, sb, name=None): 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=}') + # 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 Nrf = R*sae.sampling.N Nr = int(Nrf) - print(f'{R=} {Nr=} {Nrf=}') + # print(f'{R=} {Nr=} {Nrf=}') saed = resample(sae, Nr) - print(f'{saed=}') + # print(f'{saed=}') sced = Signal(sbe.sampling, spec=saed.spec * sbe.spec, name=name or f'{sa.name} (x) {sb.name}') return sced