Skip to content

Commit

Permalink
Merge branch 'master' of github.com:brettviren/wire-cell-python
Browse files Browse the repository at this point in the history
  • Loading branch information
@brettviren
brettviren committed Mar 22, 2024
2 parents 7eff740 + fa0af2a commit 3198a24
Show file tree
Hide file tree
Showing 3 changed files with 205 additions and 92 deletions.
193 changes: 119 additions & 74 deletions wirecell/resp/__main__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -301,8 +301,6 @@ def plot_paths(axes, rfile, pind, n, plane, trange, frange):
printer.savefig()


# lazy.
from wirecell.resp.plots import *

@cli.command("lmn-fr-plots")
@click.option("-i", "--impact", default=6*10+5, # on top of central wire of interest
Expand All @@ -320,20 +318,34 @@ def plot_paths(axes, rfile, pind, n, plane, trange, frange):
@click.option("-d", "--detector-name", default='uboone',
help="The canonical detector name")
@click.option("-r", "--field-file", default=None,
help="Explicit field response file instead of nominal one based on detector name.")
help="Explicit field response file instead of nominal one "
"based on detector name.")
@click.option("-t", "--tick", default='64*ns',
help="Resample the field response to have this sample period with units, eg '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")
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_shift)

slow_tick = unitify(slow_tick)

Tr = unitify(tick)
if vmm:
vmm = unitify_parse(vmm)
Expand All @@ -346,33 +358,41 @@ 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}")

dsigs = lmn.convolve(sigs, ers, name='I_{og} \otimes E_{og}')
dsigr = lmn.convolve(sigr, err, name='I_{rs} \otimes E_{rs}')
print(f'{dsigs.wave.shape=} {dsigr.wave.shape=}')
# 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}')

dsigr2 = lmn.interpolate(dsigs, Tr, name="(I_{og} \otimes E_{og})_{rs}")
dsigr2 = lmn.interpolate(dsigs, Tr, name=r"(I_{og} \otimes E_{og})_{rs}")

dtsigs = multiply_period(dsigs, name='T \cdot I_{og} \otimes E_{og}')
dtsigr = multiply_period(dsigr, name='T \cdot I_{rs} \otimes E_{rs}')
dtsigr2 = multiply_period(dsigr2, name="T \cdot (I_{og} \otimes E_{og})_{rs}")
dtsigs = multiply_period(dsigs, name=r'T \cdot I_{og} \otimes E_{og}')
dtsigr = multiply_period(dsigr, name=r'T \cdot I_{rs} \otimes E_{rs}')
dtsigr2 = multiply_period(dsigr2, name=r"T \cdot (I_{og} \otimes E_{og})_{rs}")

qsigs = multiply_period(sigs, name='Q_{og}')
qsigr = multiply_period(sigr, name='Q_{rs}')
qdsigs = lmn.convolve(qsigs, ers, name='Q_{og} \otimes E_{og}')
qdsigr = lmn.convolve(qsigr, err, name='Q_{rs} \otimes E_{rs}')
qdsigs = lmn.convolve(qsigs, ers, name=r'Q_{og} \otimes E_{og}')
qdsigr = lmn.convolve(qsigr, err, name=r'Q_{rs} \otimes E_{rs}')

# here we set up downsample before FRxER convolution.
sim_tick = 500*units.ns
Expand All @@ -393,13 +413,12 @@ def lmn_fr_plots(impact, plane, period,

## end building data arrays


full_range = dict(tlim=(0, sigs.sampling.T*sigs.sampling.N),
flim=(0*units.MHz, 10*units.MHz))

sigs_dur = sigs.sampling.T*sigs.sampling.N
front_range = dict(tlim=(0, 0.1*sigs_dur),
flim=(0*units.MHz, 0.03*units.MHz))
flim=(0*units.MHz, 0.03*units.MHz))
back_range = dict(tlim=(0.9*sigs_dur, 1.0*sigs_dur),
flim=(0*units.MHz, 0.03*units.MHz))

Expand All @@ -414,10 +433,11 @@ def lmn_fr_plots(impact, plane, period,
flim=(0*units.MHz, 1*units.MHz))

orglines = list()

def orgit(name):
page = 1+len(orglines)
pat = r'#+macro: %s \includegraphics[width=\textwidth, page=%s]{%s}'
text = pat%(name, page, output)
text = pat % (name, page, output)
orglines.append(text)

with pages(output) as printer:
Expand All @@ -430,63 +450,76 @@ def newpage(fig, name, title=''):
printer.savefig()
plt.clf()

frtit = f'current field response: $FR$ ({detector_name} plane {plane} impact {impact})'
frtit = 'current field response: $FR$ ' \
f'({detector_name} plane {plane} impact {impact})'

fig,_ = plot_signals((sigs, sigr), iunits='femtoampere', **full_range)
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), 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)
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')
fig, _ = plot_ends((arrs, arrr), (sigs.name, sigr.name),
iunits='femtoampere')
newpage(fig, 'fig-ends', 'current response ends')

fig,_ = plot_signals((ers, err),
iunits='mV/fC',
tlim=(0*units.us, 20*units.us),
flim=(0*units.MHz, 1*units.MHz))
fig, _ = plot_signals((ers, err),
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((dsigs, dsigr, dsigr2),
iunits='femtoampere*mV/fC',
**conv_range)
newpage(fig, 'fig-fr-er', f'convolution: $FR \otimes ER$ ({detector_name} plane {plane} impact {impact})')


fig,_ = plot_signals((dtsigs, dtsigr, dtsigr2),
iunits='us*femtoampere*mV/fC',
**conv_range)
newpage(fig, 'fig-t-fr-er', f'convolution: $T\cdot FR \otimes ER$ ({detector_name} plane {plane} impact {impact})')


fig,_ = plot_signals((qsigs, qsigr),
iunits='fC',
**zoom_kwds)
newpage(fig, 'fig-q', f'charge field response $Q = T\cdot FR$ ({detector_name} plane {plane} impact {impact})')

fig,_ = plot_signals((qdsigs, qdsigr), iunits='mV', **conv_range)
newpage(fig, 'fig-v', f'voltage response ({detector_name} plane {plane} impact {impact})')

fig,axes = plot_signals((qdsigs, qdsigr), iunits='mV', **front_range)
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((dsigs, dsigr, dsigr2),
iunits='femtoampere*mV/fC',
**conv_range)
newpage(fig, 'fig-fr-er', r'convolution: $FR \otimes ER$ '
f'({detector_name} plane {plane} impact {impact})')

fig, _ = plot_signals((dtsigs, dtsigr, dtsigr2),
iunits='us*femtoampere*mV/fC',
**conv_range)
newpage(fig, 'fig-t-fr-er', r'convolution: $T\cdot FR \otimes ER$'
f'({detector_name} plane {plane} impact {impact})')

fig, _ = plot_signals((qsigs, qsigr), iunits='fC', **zoom_kwds)
newpage(fig, 'fig-q', r'charge field response $Q = T\cdot FR$'
f'({detector_name} plane {plane} impact {impact})')

fig, _ = plot_signals((qdsigs, qdsigr), iunits='mV', **conv_range)
newpage(fig, 'fig-v', f'voltage response ({detector_name} '
f'plane {plane} impact {impact})')

fig, axes = plot_signals((qdsigs, qdsigr), iunits='mV', **front_range)
axes[0].set_ylim(-0.5e-6, 0)
newpage(fig, 'fig-v-front', f'voltage response ({detector_name} plane {plane} impact {impact}, zoom front)')
newpage(fig, 'fig-v-front', f'voltage response ({detector_name} plane'
f' {plane} impact {impact}, zoom front)')

fig,axes = plot_signals((qdsigs, qdsigr), iunits='mV', **back_range)
fig, axes = plot_signals((qdsigs, qdsigr), iunits='mV', **back_range)
tiny = 1e-11
axes[0].set_ylim(-tiny, tiny)
newpage(fig, 'fig-v-back', f'voltage response ({detector_name} plane {plane} impact {impact}, zoom back)')
print('q=', 100*numpy.sum(qdsigs.wave) / numpy.sum(numpy.abs(qdsigs.wave)), '%')
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)), '%')

# ER in fast and slow binning
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)
newpage(fig, 'fig-fr-ds', 'FR(I) (downsample vs decimate) $\leftrightarrow$ convolve')
Expand All @@ -502,35 +535,46 @@ def newpage(fig, name, title=''):

# yet more checks

colors=["black","red","blue","green","yellow"]
colors = ["black", "red", "blue", "green", "yellow"]

# Check shift
tshift_ns = 1

# current
ss = lmn.Sampling(T=sigr.sampling.T, N=sigr.sampling.N, t0=tshift_ns*units.ns)
sigr_shift = lmn.Signal(ss, wave=sigr.wave, name='I_{rs} + %dns'%tshift_ns)
ss = lmn.Sampling(T=sigr.sampling.T, N=sigr.sampling.N,
t0=tshift_ns*units.ns)
sigr_shift = lmn.Signal(ss, wave=sigr.wave,
name='I_{rs} + %dns' % tshift_ns)
diff_sigs = [sigs, sigr, sigr_shift]
ndiff_sigs = len(diff_sigs)
fig, axes = plot_wave_diffs(diff_sigs, xlim=zoom_kwds['tlim'], marker='.',
per = [dict(markersize=ndiff_sigs-ind,
color=colors[ind]) for ind in range(ndiff_sigs)])
newpage(fig, 'fig-fr-diff', f'FR differences ({detector_name} plane {plane} impact {impact})')

fig, axes = plot_wave_diffs(
diff_sigs, xlim=zoom_kwds['tlim'], marker='.',
per=[dict(markersize=ndiff_sigs-ind,
color=colors[ind]) for ind in range(ndiff_sigs)])
newpage(fig, 'fig-fr-diff', 'FR differences '
f'({detector_name} plane {plane} impact {impact})')

# coldelec
ss = lmn.Sampling(T=err.sampling.T, N=err.sampling.N, t0=tshift_ns*units.ns)
err_shift = lmn.Signal(ss, wave=err.wave, name='I_{rs} + %dns'%tshift_ns)
ss = lmn.Sampling(T=err.sampling.T, N=err.sampling.N,
t0=tshift_ns*units.ns)
err_shift = lmn.Signal(ss, wave=err.wave,
name='I_{rs} + %dns'%tshift_ns)
diff_sigs = [ers, err, err_shift]
ndiff_sigs = len(diff_sigs)
fig, axes = plot_wave_diffs(diff_sigs, xlim=zoom_kwds['tlim'], marker='.',
per = [dict(markersize=ndiff_sigs-ind,
color=colors[ind]) for ind in range(ndiff_sigs)])
newpage(fig, 'fig-fr-diff', f'FR differences ({detector_name} plane {plane} impact {impact})')
fig, axes = plot_wave_diffs(
diff_sigs, xlim=zoom_kwds['tlim'], marker='.',
per=[dict(markersize=ndiff_sigs-ind,
color=colors[ind]) for ind in range(ndiff_sigs)])
newpage(fig, 'fig-fr-diff', 'FR differences '
f'({detector_name} plane {plane} impact {impact})')

# voltage
ss = lmn.Sampling(T=qdsigr.sampling.T, N=qdsigr.sampling.N, t0=tshift_ns*units.ns)
vr_shift = lmn.Signal(ss, wave=qdsigr.wave, name='Q_{rs} \otimes E_{rs} + %dns'%tshift_ns)

ss = lmn.Sampling(T=qdsigr.sampling.T, N=qdsigr.sampling.N,
t0=tshift_ns*units.ns)
vr_shift = lmn.Signal(
ss, wave=qdsigr.wave,
name=r'Q_{rs} \otimes E_{rs} + %dns' % tshift_ns)

diff_sigs = [qdsigs, qdsigr, vr_shift]
ndiff_sigs = len(diff_sigs)
fig, axes = plot_wave_diffs(diff_sigs, xlim=zoom_kwds['tlim'], marker='.',
Expand All @@ -545,9 +589,10 @@ def newpage(fig, name, title=''):
oo.write('\n'.join(orglines))
oo.write('\n')


def main():
cli(obj=dict())


if '__main__' == __name__:
main()

21 changes: 14 additions & 7 deletions wirecell/util/__main__.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ def cli(ctx):
'''
pass


@cli.command("lmn")
@click.option("-n", "--sampling-number", default=None, type=int,
help="Original number of samples.")
Expand All @@ -28,7 +29,8 @@ def cli(ctx):
@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")
help="Precision by which integer and "
"rationality conditions are judged")
def cmd_lmn(sampling_number, sampling_period, resampling_period, error):
'''Print various LMN parameters for a given resampling.
Expand All @@ -44,8 +46,12 @@ def cmd_lmn(sampling_number, sampling_period, resampling_period, error):
from wirecell.util import lmn
print(f'initial sampling: {Ns=} {Ts=}')

gcd = lmn.egcd(Tr, Ts-Tr)
print(f'egcd({Ts=}, {Tr=}): {gcd}')
dn = lmn.rational_deltan(Ts, Tr)
print(f'minimum delta-n: {dn}')
nrat = lmn.rational_size(Ts, Tr, error)
print(f'rationality factor: {nrat}')
print(f'minimum size: {nrat}')

nrag = Ns % nrat
if nrag:
Expand All @@ -57,19 +63,20 @@ def cmd_lmn(sampling_number, sampling_period, resampling_period, error):
Ns_rational = Ns

Nr_target = Ns_rational*Ts/Tr
assert abs(Nr_target - round(Nr_target)) < error # assured by rational_size()
assert abs(Nr_target - round(Nr_target)) < error
Nr_target = round(Nr_target)
print(f'resampling target: {Nr_target=} {Tr=}')

ndiff = Nr_target - Ns_rational
print(f'final resampling: {Ns_rational=}, {Ts=} -> {Nr_target=}, {Tr=} diff of {ndiff}')

print(f'final resampling: {Ns_rational=}, {Ts=} '
f'-> {Nr_target=}, {Tr=} diff of {ndiff}')

Nr_wanted = Ns*Ts/Tr
Nr = math.ceil(Nr_wanted)
if abs(Nr_wanted - Nr) > error:
print(f'--> warning: noninteger {Nr_wanted=} for {Tr=}. Duration will change from {Ns*Ts} to {Nr*Tr} due to rounding.')
print(f'requested resampling target: {Nr=} {Tr=}')
print('--> note: noninteger nominal target size:'
f' {Nr_wanted} for {Ns=} {Ts=} {Tr=}.')
print(f'resampling target size: {Nr=} {Tr=}')

ntrunc = Nr - Nr_target
if ntrunc < 0:
Expand Down
Loading

0 comments on commit 3198a24

Please sign in to comment.