Simulated Ground Telescopes¶
This tutorial focuses on simulating data from a ground-based telescope. We first create a fake telescope with a synthetic focalplane located in Chile. Then we create a synthetic observing schedule and use that to scan the sky. Later notebooks make use of helper functions that create generic focalplanes, but in this example we use low-level functions to customize things a bit more.
# Optionally change logging level
import os
os.environ["TOAST_LOGLEVEL"] = "INFO"
# This is needed before importing toast, and should
# match the value passed to the '-t' option of %toast
os.environ["OMP_NUM_THREADS"] = "4"
# TOAST interactive startup
import toast.interactive
%load_ext toast.interactive
%toast -p 1 -t 4 -a
# Built-in modules
import sys
import os
import re
import datetime
import shutil
# External modules
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.image import imread
import astropy.units as u
from astropy.table import Row, QTable
import healpy as hp
# TOAST
import toast
import toast.schedule_sim_ground
from toast.instrument_sim import plot_focalplane
from toast.tests import helpers
from toast.observation import default_values as defaults
# Display inline plots
%matplotlib inline
from IPython.display import Image
from IPython.display import display
# MPI communicator
world, procs, rank = toast.mpi.get_world()
comm = helpers.create_comm(world, single_group=True)
# Output directory for this tutorial
topdir = "out_sim_ground"
if rank == 0 and os.path.exists(topdir):
shutil.rmtree(topdir)
out_dir = helpers.create_outdir(world, topdir=topdir)
Helper Functions¶
Here are a few functions we will use later in the notebook.
def plot_dets(obs, d_start=0, d_end=None, s_start=0, s_end=None, view=None, signal=defaults.det_data):
"""Plot some detectors in an observation.
Args:
obs (Observation): The observation
d_start (int): The starting local detector index to plot.
d_end (int): The local detector index limit to plot.
s_start (int): The starting sample index to plot.
s_end (int): The sample index limit to plot
view (str): The optional intervals to overplot.
signal (str): The detdata name to plot.
"""
slc = slice(s_start, s_end, 1)
fig = plt.figure(dpi=100, figsize=(18, 12))
ax = fig.add_subplot(2, 1, 1, aspect="auto")
plt.gca().set_prop_cycle(None)
for idet, det in enumerate(obs.select_local_detectors(flagmask=defaults.det_mask_nonscience)):
if idet < d_start:
continue
if d_end is not None and idet >= d_end:
continue
ax.plot(
obs.shared[defaults.times].data[slc],
obs.detdata[signal][det, slc],
'-',
label=det,
)
ax.legend(loc="best")
ax = fig.add_subplot(2, 1, 2, aspect="auto")
if view is not None:
inview = np.zeros_like(obs.shared[defaults.shared_flags].data[slc])
begin = [x.first for x in obs.intervals[view]]
end = [x.last+1 for x in obs.intervals[view]]
for b, e in zip(begin, end):
inview[b:e] = 1
ax.plot(
obs.shared[defaults.times].data[slc],
inview,
'-',
color="red",
label=f"View {view}",
)
ax.plot(
obs.shared[defaults.times].data[slc],
obs.shared[defaults.shared_flags].data[slc],
'-',
color="black",
label="Shared Flags",
)
plt.gca().set_prop_cycle(None)
for idet, det in enumerate(obs.select_local_detectors(flagmask=defaults.det_mask_nonscience)):
if idet < d_start:
continue
if d_end is not None and idet >= d_end:
continue
ax.plot(
obs.shared[defaults.times].data[slc],
obs.detdata[defaults.det_flags][det, slc],
'-',
label=det,
)
ax.legend(loc="best")
plt.show()
plt.close()
def plot_scanning(obs, s_start=0, s_end=None):
slc = slice(s_start, s_end, 1)
times = obs.shared[defaults.times].data[slc]
az = obs.shared[defaults.azimuth].data[slc]
el = obs.shared[defaults.elevation].data[slc]
fig = plt.figure(dpi=100, figsize=(18, 12))
ax = fig.add_subplot(2, 1, 1, aspect="auto")
ax.plot(times, az, label="Azimuth")
ax.set_xlabel("Posix Timestamps")
ax.set_ylabel("Azimuth")
ax = fig.add_subplot(2, 1, 2, aspect="auto")
ax.plot(times, el, label="Elevation")
ax.set_xlabel("Posix Timestamps")
ax.set_ylabel("Elevation")
plt.show()
plt.close()
def plot_noise_model(model, model_fit=None, d_start=0, d_end=None):
fig = plt.figure(dpi=100, figsize=(18, 12))
ax = fig.add_subplot(1, 1, 1)
plt.gca().set_prop_cycle(None)
plot_max = 0
plot_min = 1e100
for idet, det in enumerate(model.detectors):
if idet < d_start:
continue
if d_end is not None and idet >= d_end:
continue
freq = model.freq(det).to_value(u.Hz)
psd = model.psd(det).to_value(u.K**2 * u.s)
plot_min = min(plot_min, np.amin(psd))
plot_max = max(plot_max, np.amax(psd))
ax.loglog(
freq,
psd,
label=det,
)
if model_fit is not None:
# Also plot the fit
plt.gca().set_prop_cycle(None)
for idet, det in enumerate(model.detectors):
if idet < d_start:
continue
if d_end is not None and idet >= d_end:
continue
freq = model_fit.freq(det)
psd = model_fit.psd(det)
ax.loglog(
freq.to_value(u.Hz),
psd.to_value(u.K**2 * u.s),
label=f"{det} Fit",
)
freq = model.freq(model.detectors[0])
ax.set_xlim(freq[0].to_value(u.Hz), freq[-1].to_value(u.Hz))
ax.set_ylim(0.9 * plot_min, 1.1 * plot_max)
ax.set_xlabel("Frequency [Hz]")
ax.set_ylabel("PSD [K$^2$ / Hz]")
ax.legend(loc="best")
plt.show()
plt.close()
Fake Telescope¶
We create just a small number of detectors here since we are running this notebook serially. If you use more processes you can increase the number of detectors on the focalplane. First create a Site for telescope:
site = toast.instrument.GroundSite("atacama", "-22:57:30", "-67:47:10", 5200.0 * u.meter)
Now we will create a focalplane consisting of three rhombus wafers packed into a hexagon.
fp_fwhm = 30.0 * u.arcmin
focalplane = toast.instrument_sim.fake_rhombihex_focalplane(
n_pix_rhombus=16,
width=8.0 * u.degree,
gap=0 * u.radian,
sample_rate=10.0 * u.Hz,
epsilon=0.0,
fwhm=fp_fwhm,
bandcenter=150 * u.GHz,
bandwidth=20 * u.GHz,
psd_net=300.0 * u.uK * np.sqrt(1 * u.second),
psd_fmin=1.0e-5 * u.Hz,
psd_alpha=1.0,
psd_fknee=0.05 * u.Hz,
fwhm_sigma=0.0 * u.arcmin,
bandcenter_sigma=0 * u.GHz,
bandwidth_sigma=0 * u.GHz,
random_seed=123456,
)
fov = focalplane.field_of_view
# Look at the table of detector properties
focalplane.detector_data
# Make a plot of this focalplane layout.
detpolcol = {
x: "red" if re.match(r".*A-.*", x) is not None else "blue" for x in focalplane.detectors
}
if rank == 0:
plot_focalplane(
focalplane=focalplane,
width=1.2 * fov,
height=1.2 * fov,
show_labels=True,
pol_color=detpolcol
)
Atmospheric Monitoring¶
In order to provide a channel to monitor the atmospheric water content, we add a single detector at the boresight whose bandpass is centered on the water line at 183GHz. We make a copy of the previous detector table and construct a new focalplane.
det_props = QTable(focalplane.detector_data)
# Copy the last row into a dictionary
atm_det = {x: det_props[-1][x] for x in det_props.colnames}
print(atm_det)
# Modify the atmosphere detector properties
atm_det["name"] = "ATM0"
atm_det["quat"] = np.array([0.0, 0.0, 0.0, 1.0])
atm_det["bandcenter"] = 183.0 * u.GHz
atm_det["bandwidth"] = 20.0 * u.GHz
det_props.add_row(atm_det)
# Build a new focalplane with the updated table
full_fp = toast.instrument.Focalplane(
detector_data=det_props,
sample_rate=focalplane.sample_rate,
)
detpolcol = {
x: "red" if re.match(r".*A-.*", x) is not None else "blue" for x in full_fp.detectors
}
if rank == 0:
plot_focalplane(
focalplane=full_fp,
width=1.2 * fov,
height=1.2 * fov,
show_labels=True,
pol_color=detpolcol
)
We can check the top-hat bandpasses in this Focalplane. We just look at the first normal detector and the last detector which is the atmosphere monitor.
if rank == 0:
fig = plt.figure(dpi=100, figsize=(12, 6))
ax = fig.add_subplot(1, 1, 1, aspect="auto")
for det in [full_fp.detectors[0], full_fp.detectors[-1]]:
freq = full_fp.bandpass.freqs(det)
bpass = full_fp.bandpass.bandpass(det)
ax.plot(
freq,
bpass,
'-',
label=det,
)
ax.set_xlim(100e9, 250e9)
ax.legend(loc="best")
plt.show()
plt.close()
Finally we build our telescope with this updated focalplane and Site
telescope = toast.instrument.Telescope("telescope", focalplane=full_fp, site=site)
Simulated Observing Schedule¶
Now that we have a telescope, we create an observing schedule.
schedule = None
if rank == 0:
tdir = out_dir
if tdir is None:
tdir = tempfile.mkdtemp()
sch_file = os.path.join(tdir, "ground_schedule.txt")
toast.schedule_sim_ground.run_scheduler(
opts=[
"--site-name",
telescope.site.name,
"--telescope",
telescope.name,
"--site-lon",
"{}".format(telescope.site.earthloc.lon.to_value(u.degree)),
"--site-lat",
"{}".format(telescope.site.earthloc.lat.to_value(u.degree)),
"--site-alt",
"{}".format(telescope.site.earthloc.height.to_value(u.meter)),
"--patch",
"bossn,1,-180,15,-140,2",
"--start",
"2025-02-21 00:00:00",
"--stop",
"2025-02-23 00:00:00",
"--out",
sch_file,
"--equalize-time",
"--patch-coord",
"C",
"--el-min",
"40",
"--el-max",
"70",
"--sun-el-max",
"90",
"--sun-avoidance-angle",
"30",
"--moon-avoidance-angle",
"0",
"--ces-max-time",
"36000",
"--fp-radius",
"0",
"--boresight-angle-step",
"180",
"--boresight-angle-time",
"1440",
"--time-step-s",
"600",
"--lock-az-range",
"--elevations",
"40,50,60,70",
]
)
schedule = toast.schedule.GroundSchedule()
schedule.read(sch_file)
if out_dir is None:
shutil.rmtree(tdir)
if world is not None:
schedule = world.bcast(schedule, root=0)
Simulated Observing¶
Now we use this schedule to create some fake observing with our telescope. This will generate the data containers with boresight pointing, but the detector data is still zero.
# Start with an empty data container
data = toast.Data(comm)
# Populate observations according to the schedule and telescope.
sim_ground = toast.ops.SimGround(
telescope=telescope,
weather="atacama",
detset_key="pixel",
schedule=schedule,
median_weather=True, # no longer random weather, but less chance of an outlier
)
sim_ground.apply(data)
# Print out the result.
data.info()
if rank == 0:
plot_scanning(data.obs[0], s_start=0, s_end=2000)
Simulated Detector Data¶
Now we will simulate several components of our detector data. Before doing that, we set up some operators that compute our detector pointing and response on the sky (the "pointing matrix"). In TOAST, the pointing matrix is split into the pixelization of detector samples and the Stokes weights (response to I/Q/U/V on the sky).
Detector Pointing¶
There are 3 types of operators that define the detector pointing. The first is the geometric offset from the boresight coordinate frame to the detector coordinate frame. In the simplest case this is just a quaternion for each detector (stored in the focalplane table). Once the geometric detector pointing is computed, the pixelization on the sky is specified by a separate operator. Finally, the response of the detector to incoming Stokes parameters is given by another operator.
# Geometric detector pointing from boresight frame to detector frame. We define these
# for both horizontal and equatorial coordinates, since we need the detector pointing
# in horizontal coordinates for atmosphere simulation below.
det_point_azel = toast.ops.PointingDetectorSimple(
boresight=defaults.boresight_azel,
quats="quats_azel"
)
det_point_radec = toast.ops.PointingDetectorSimple(
boresight=defaults.boresight_radec,
quats="quats_radec"
)
# Pixelization. Choose a coarse pixelization for this exercise since there is
# a small patch and only a few detectors.
nside = 256
pixels_radec = toast.ops.PixelsHealpix(
nside=nside,
nest=True,
detector_pointing=det_point_radec,
)
# Stokes weights. This just uses focalplane table properties to treat each detector
# as a linear polarizer with possibly some cross-polar response.
weights_azel = toast.ops.StokesWeights(
mode="IQU",
detector_pointing=det_point_azel,
)
weights_radec = toast.ops.StokesWeights(
mode="IQU",
detector_pointing=det_point_radec,
)
Pixel Distribution¶
When working with sky data for both simulations and mapmaking, each process only stores pixels which are "hit" by the local detectors on that process. Computing this "pixel distribution" requires passing through the pointing. Normally this is done without saving the detector pointing (for memory considerations). In this notebook, we just compute the full detector pointing once at the beginning and save it.
pixels_radec.apply(data)
weights_radec.apply(data)
pix_dist = toast.ops.BuildPixelDistribution(
pixel_dist="pixel_dist",
pixel_pointing=pixels_radec,
)
pix_dist.apply(data)
Synthetic Sky¶
In order to have some kind of sky signal in our data, we generate a fake sky and scan that into timestreams.
input_map_file = os.path.join(out_dir, "fake_sky.fits")
if rank == 0:
c_ell = helpers.fetch_nominal_cmb_cls(
out_file=os.path.join(os.path.dirname(out_dir), "cl_nominal.txt")
)
input_map_ring = hp.synfast(c_ell, nside, fwhm=fp_fwhm.to_value(u.radian))
input_map = 1.0e-6 * hp.reorder(input_map_ring, inp="RING", out="NEST")
hp.write_map(input_map_file, input_map, nest=True)
hp.mollview(input_map[0], nest=True, min=-0.01, max=0.01)
hp.gnomview(
input_map[0], nest=True, min=-0.01, max=0.01, rot=(200.156, 8.466), reso=4.0, xsize=1600
)
hp.mollview(input_map[1], nest=True, min=-0.0002, max=0.0002)
hp.gnomview(
input_map[1], nest=True, min=-0.0002, max=0.0002, rot=(200.156, 8.466), reso=4.0, xsize=1600
)
hp.mollview(input_map[2], nest=True, min=-0.0002, max=0.0002)
hp.gnomview(
input_map[2], nest=True, min=-0.0002, max=0.0002, rot=(200.156, 8.466), reso=4.0, xsize=1600
)
# Scan the map
scan_map = toast.ops.ScanHealpixMap(
file=input_map_file,
pixel_pointing=pixels_radec,
stokes_weights=weights_radec,
)
scan_map.apply(data)
# Plot the last few detectors
if rank == 0:
plot_dets(data.obs[0], d_start=90, d_end=None, s_start=0, s_end=2000, view="scanning")
Instrumental Noise¶
We create a trivial noise model using nominal parameters from the focalplane table and then use this noise model to simulate timestreams.
nominal_noise = toast.ops.DefaultNoiseModel()
nominal_noise.apply(data)
# Plot this nominal noise model for the last few detectors
if rank == 0:
plot_noise_model(
data.obs[0][nominal_noise.noise_model],
model_fit=None,
d_start=50,
d_end=None
)
sim_noise = toast.ops.SimNoise(
noise_model=nominal_noise.noise_model,
)
sim_noise.apply(data)
# Plot the last few detectors
if rank == 0:
plot_dets(data.obs[0], d_start=90, d_end=None, s_start=0, s_end=2000, view="scanning")
We can see that the "atmospheric monitor" detector does not look much different here, since we have only simulated detector noise.
Ground Pickup¶
The ground and environment around the telescope may produce different loading as a function of azimuth as the telescope scans. This kind of signal is one of several possible "scan synchronous signals".
ground_pickup = toast.ops.SimScanSynchronousSignal(
detector_pointing=det_point_azel,
scale=0.001 * u.K,
stokes_weights=weights_azel,
)
ground_pickup.apply(data)
# Plot the last few detectors
if rank == 0:
plot_dets(data.obs[0], d_start=90, d_end=None, s_start=0, s_end=2000, view="scanning")
Simulated Atmosphere¶
Now we will simulate a 3D atmospheric slab moving in front of the telescope and integrate each detector along the line of site and over its bandpass. In order to do this, we have to define an operator which knows how to compute detector pointing in Az/El coordinates. This just uses the boresight pointing and the detector quaternion rotations from the boresight to compute detector pointing.
sim_atm = toast.ops.SimAtmosphere(
detector_pointing=det_point_azel,
add_loading=True,
lmin_center=0.001 * u.m,
lmin_sigma=0.0001 * u.m,
lmax_center=1.0 * u.m,
lmax_sigma=0.1 * u.m,
xstep=20 * u.m,
ystep=20 * u.m,
zstep=20 * u.m,
zmax=200 * u.m,
gain=4e-5,
wind_dist=1000 * u.m,
)
sim_atm.apply(data)
# Plot the last few normal detectors
if rank == 0:
plot_dets(data.obs[0], d_start=90, d_end=96, s_start=0, s_end=2000, view="scanning")
# Plot the atmosphere monitor
if rank == 0:
plot_dets(data.obs[0], d_start=96, d_end=None, s_start=0, s_end=2000, view="scanning")
Here we see that the atmospheric monitor has substantially more power from the 183GHz water line.
Map Making¶
First we want to flag the atmosphere monitor channel so that it is not considered "science" data for mapmaking purposes.
for ob in data.obs:
ob.update_local_detector_flags({"ATM0": defaults.det_mask_processing})
Filtering¶
Since we are using the template regression / destriping mapmaker below, we want to do minimal filtering of the timestreams.
toast.ops.CommonModeFilter(
redistribute=False,
regress=True,
).apply(data)
Noise Estimation¶
The original noise estimate (used above to simulate instrument noise) only captures the nominal readout and detector noise sources. For map making we will treat the timestream as noise-dominated and estimate the noise model directly from the timestream. First, create a raw, binned estimate of the PSD in each detector:
# Estimate noise
estim = toast.ops.NoiseEstim(
out_model="noise_estimate",
lagmax=100,
nbin_psd=32,
nsum=1,
)
estim.apply(data)
This raw estimate can produce undesired effects if we use it directly in mapmaking. Instead, we first fit an analytic 1/f noise model to these.
# Compute a 1/f fit to this
noise_fitter = toast.ops.FitNoiseModel(
noise_model=estim.out_model,
out_model="fit_noise_model",
)
noise_fitter.apply(data)
# Plot this nominal noise model for the last few detectors
if rank == 0:
plot_noise_model(
data.obs[0][estim.out_model],
model_fit=data.obs[0][noise_fitter.out_model],
d_start=90,
d_end=None
)
Binning Operator¶
A central piece of the mapmaking is the "binning" of timestreams into maps. This process accumulates the "noise weighted map" and then multiplies this by the diagonal pixel covariance:
$$ \text{Binned Map} = \left(P^T N^{-1} P\right)^{-1} P^T N^{-1} d $$
You can see from this that in addition to the input timestream data we need the estimated noise model and the pointing matrix.
# Set up binning operator for solving
binner = toast.ops.BinMap(
pixel_dist=pix_dist.pixel_dist,
pixel_pointing=pixels_radec,
stokes_weights=weights_radec,
noise_model=noise_fitter.out_model,
)
Template Matrix¶
The TOAST mapmaker is a "generalized destriper" that solves for "template amplitudes". These templates represent anything in the time ordered data which is not fixed on the sky or simply white noise. For this example, we will use 2 templates. One to model the scan synchronous signal and one to model the 1/f noise (including the atmosphere).
# The Offset template models 1/f noise as a stepwise function, which
# is the same as a "classic" destriper.
tmpl_offset = toast.templates.Offset(
times=defaults.times,
noise_model=noise_fitter.out_model,
step_time=1.0 * u.second,
)
# Build a template matrix with our templates.
tmatrix = toast.ops.TemplateMatrix(
templates=[tmpl_offset],
)
Making the Map¶
Now we are ready to instantiate the mapmaker operator. Note that if we do not specify the template matrix, then this will just produce a binned map.
map_maker = toast.ops.MapMaker(
name="mapmaker",
binning=binner,
template_matrix=tmatrix,
solve_rcond_threshold=1.0e-1,
map_rcond_threshold=1.0e-1,
iter_min=200,
iter_max=300,
write_hits=True,
write_map=True,
write_binmap=True,
write_cov=False,
write_invcov=False,
write_rcond=True,
output_dir=out_dir,
keep_solver_products=True, # We set this to True so we can plot solved template amplitudes later
)
map_maker.apply(data)
Now plot the output maps.
# The output filenames will be based on the name of the mapmaker operator
out_root = os.path.join(out_dir, map_maker.name)
# Helper functions to plot all the maps
def plot_maps(
root,
range_I=(-0.01, 0.01),
range_Q=(-0.0002, 0.0002),
range_U=(-0.0002, 0.0002),
max_hits=1000,
truth=None
):
cmap = "viridis"
gnomres = 8.0
gnomrot = (199.5, 8.3)
xsize = 800
hits_file = f"{root}_hits.fits"
rcond_file = f"{root}_rcond.fits"
binmap_file = f"{root}_binmap.fits"
map_file = f"{root}_map.fits"
# Load hits
hits = hp.read_map(hits_file, field=None, nest=True)
goodhits = hits > 0
badhits = np.logical_not(goodhits)
# Load rcond
rcond = hp.read_map(rcond_file, field=None, nest=True)
rcond[badhits] = hp.UNSEEN
# Maps
maps = hp.read_map(map_file, field=None, nest=True)
binmaps = hp.read_map(binmap_file, field=None, nest=True)
resid = None
resid_bin = None
if truth is not None:
truth_maps = hp.read_map(truth, field=None, nest=True)
resid = list()
resid_bin = list()
for i in range(3):
resid.append(np.array(maps[i]) - truth_maps[i])
resid_bin.append(np.array(binmaps[i]) - truth_maps[i])
for i in range(3):
maps[i][badhits] = hp.UNSEEN
binmaps[i][badhits] = hp.UNSEEN
if truth is not None:
truth_maps[i][badhits] = hp.UNSEEN
resid[i][badhits] = hp.UNSEEN
resid_bin[i][badhits] = hp.UNSEEN
# Plot hits and rcond
fig = plt.figure(dpi=100, figsize=(18, 12))
hp.gnomview(
map=hits,
fig=fig.number,
sub=(1, 2, 1),
rot=gnomrot,
xsize=xsize,
reso=gnomres,
nest=True,
cmap=cmap,
min=0,
max=max_hits,
title="Hits",
)
hp.gnomview(
map=rcond,
fig=fig.number,
sub=(1, 2, 2),
rot=gnomrot,
xsize=xsize,
reso=gnomres,
nest=True,
cmap=cmap,
min=0,
max=0.5,
title="Inverse Condition Number",
)
plt.show()
plt.close()
# Plot maps
plot_cols = 2
if truth is not None:
plot_cols = 4
plot_rows = 3
fig = plt.figure(dpi=100, figsize=(18, 18))
counter = 1
for row, (stokes, rng) in enumerate([("I", range_I), ("Q", range_Q), ("U", range_U)]):
for mps, res, name in [(maps, resid, "Destriped"), (binmaps, resid_bin, "Binned")]:
hp.gnomview(
map=mps[row],
fig=fig.number,
sub=(plot_rows, plot_cols, counter),
rot=gnomrot,
xsize=xsize,
reso=gnomres,
nest=True,
cmap=cmap,
min=rng[0],
max=rng[1],
title=f"{name} Stokes {stokes}",
)
counter += 1
if truth is not None:
hp.gnomview(
map=res[row],
fig=fig.number,
sub=(plot_rows, plot_cols, counter),
rot=gnomrot,
xsize=xsize,
reso=gnomres,
nest=True,
cmap=cmap,
min=rng[0],
max=rng[1],
title=f"{name} Stokes {stokes} Minus Input",
)
counter += 1
plt.show()
plt.close()
plot_maps(out_root, truth=input_map_file)
The solved template amplitudes will usually have degeneracies (for example, the ground pickup across one left-right scan could also be interpreted as just 1/f noise). However, we are just concerned with capturing the degrees of freedom of the relevant non-sky signal content so that it does not contaminate the map. This is something to keep in mind as we plot the solved template amplitudes.
# Write solved offset amplitudes
oamps = data[f"{map_maker.name}_solve_amplitudes"][tmpl_offset.name]
oroot = os.path.join(out_dir, f"{map_maker.name}_offset")
tmpl_offset.write(oamps, oroot)
if rank == 0:
for ob in data.obs:
toast.templates.offset.plot(
f"{oroot}_{ob.name}.h5",
compare={x: ob.detdata[defaults.det_data][x, :] for x in ob.local_detectors},
out=f"{oroot}_{ob.name}",
xlim=(0, 1000),
)