Load and scale 1/f and white noise simulations
import logging as log
thinfp_table
def get_thinfp(channel):
"""Get the focal plane thinning factor for noise simulations
Parameters
----------
channel : str
CMB-S4 channel tag e.g. HFL2
Returns
-------
thinfp : int
thinning factor
"""
return (thinfp_table[thinfp_table["tube"] == channel[:-1]])["thinfp"][0]
assert get_thinfp("ULFPL1") == 1
assert get_thinfp("HFS2") == 8
assert get_thinfp("MFLS1") == 4
def get_tube_years(config, site, channel):
"""Compute the number of tube/years in the CMB-S4 configuration
config_telescopes : dict
CMB-S4 telescopes configuration,
generally loaded from a TOML file
site : str
'Pole' or 'Chile', case doesn't matter
channel : str
Channel tag, e.g. 'MFHS1'
"""
tube_years = 0
for telescope_name, telescope_config in config["telescopes"][
get_telescope(channel)
].items():
if telescope_config["site"].lower() == site.lower():
num_tubes = telescope_config.get(channel[:-1], 0)
tube_years += num_tubes * telescope_config.get(
"years", config["experiment"]["total_experiment_length_years"]
)
return tube_years
s4 = read_instrument_model()
!tail -30 s4_design.toml
import toml
config = toml.load("s4_design.toml")
log.basicConfig(level=log.INFO)
for site in ["Pole", "Chile"]:
for row in s4:
channel = row["band"]
tube_years = get_tube_years(config, site, channel)
print(site, channel, tube_years)
telescope = get_telescope(channel)
if site == "Chile":
if telescope == "SAT":
assert tube_years == 0, "All SAT at Pole"
elif channel.startswith("ULFL"):
assert tube_years == 0, "No ULFL in Chile"
else:
if channel.startswith("MFL"):
assert tube_years == 54*2*7, "54+54 MFL channels"
if site == "Pole":
if telescope == "SAT":
if channel.startswith("LFS"):
assert tube_years == 2*7, "2 LFS Tubes"
elif channel.startswith("MF"):
assert tube_years == 6*7, "2 LFS Tubes"
else:
if channel.startswith("HFPL"):
assert tube_years == 18*7, "18 HFL tubes"
def load_noise(config, site, channel, realization=0):
"""Load noise maps for a channel
Parameters
----------
config : dict
CMB-S4 configuration,
generally loaded from a TOML file
site : str
'Pole' or 'Chile', case doesn't matter
channel : str
Channel tag, e.g. 'MFHS1'
realization : int
Choose one of the available 8 realizations
Returns
-------
output_map : numpy array
Output map with all emissions combined, uses nan for missing pixels
"""
telescope = get_telescope(channel)
channel_noP = channel.replace("P", "")
map_filename = (
Path(f"{realization:08d}")
/ f"{site.lower()}_noise_{telescope}_{channel_noP}_{mapmaking_naming[telescope]}"
)
log.info(f"Base folder: {base_folder}")
log.info(f"Reading {map_filename}")
output_map = hp.read_map(
Path(base_folder) / map_filename, (0, 1, 2), dtype=None, verbose=False
)
output_map[output_map == hp.UNSEEN] = np.nan
output_map[output_map == 0] = np.nan
# input map is 10 days
output_map *= np.sqrt(
10
* simulations_observing_efficiency[site.lower()][channel]
/ (
365.25
* get_observing_efficiency(
config["experiment"]["observing_efficiency"], site, telescope, channel
)
* get_observing_efficiency(
config["experiment"]["sensitivity_factor"], site, telescope, channel
)
)
)
output_map /= np.sqrt(
get_tube_years(config, site, channel) / simulated_tubes[channel[:-1]]
)
# focal plane thinning factor of TOD simulations
output_map /= np.sqrt(get_thinfp(channel))
return output_map
filenames = !ls $base_folder/00000000/*noise*bmap*
import os.path
for f in map(os.path.basename, filenames):
print(f)
channel = "HFS1"
site = "Pole"
telescope = get_telescope(channel)
input_map = hp.read_map(
base_folder + "/00000000/" + \
f"{site.lower()}_noise_{telescope}_{channel}_telescope_all_time_all_filtered.fits.gz"
, (0,1,2))
input_map[input_map == 0] = hp.UNSEEN
output_map = load_noise(config, site, channel, realization=0)
output_map[np.isnan(output_map)] = hp.UNSEEN
np.testing.assert_allclose(input_map == hp.UNSEEN, output_map == hp.UNSEEN)
simulations_observing_efficiency
observing_efficiency = get_observing_efficiency(
config["experiment"]["observing_efficiency"], site, telescope, channel
) * get_observing_efficiency(
config["experiment"]["sensitivity_factor"], site, telescope, channel
)
observing_efficiency
np.testing.assert_allclose(input_map[input_map != hp.UNSEEN] * \
np.sqrt(10 * simulations_observing_efficiency[site.lower()][channel] / 8 / (7 * 365.25 * observing_efficiency) / 4),
output_map[output_map != hp.UNSEEN], rtol=1e-6)
%matplotlib inline
import matplotlib.pyplot as plt
hp.mollview(output_map[0], min=-1e-6, max=1e-6, unit="K", title="Noise I")
hp.mollview(output_map[1], min=-1e-5, max=1e-5, unit="K", title="Noise Q")
hp.mollview(output_map[2], min=-1e-5, max=1e-5, unit="K", title="Noise U")
