Load and scale 1/f and white noise simulations
import logging as log
thinfp_table
Table length=11
telescopetubethinfp
str3str5int64
LATULFPL1
LATLFL4
LATLFPL4
LATMFL16
LATHFL16
LATMFPL16
LATHFPL16
SATLFS1
SATMFLS4
SATMFHS4
SATHFS8

get_thinfp[source]

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
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

get_tube_years[source]

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'
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
Warning: Python module not loaded, you already have Python loaded via conda init
# we can also specify defaults for sites (notice capitalization) or telescopes
# observing_efficiency.Pole.SAT.default = 0.2
# observing_efficiency.Chile.default = 0.2

[telescopes]

    # Telescope names are not used within the software

    [telescopes.SAT]
    # Tubes for each SAT should be 3

    SAT1 = { LFS=0, MFLS= 3, MFHS=0, HFS=0, site="Pole", years=7 }
    SAT3 = { LFS=0, MFLS= 0, MFHS=3, HFS=0, site="Pole", years=7 }
    SAT5 = { LFS=1, MFLS= 0, MFHS=0, HFS=2, site="Pole", years=7 }

    SAT0 = { LFS=0, MFLS= 3, MFHS=0, HFS=0, site="Pole", years=7 }
    SAT2 = { LFS=0, MFLS= 0, MFHS=3, HFS=0, site="Pole", years=7 }
    SAT4 = { LFS=1, MFLS= 0, MFHS=0, HFS=2, site="Pole", years=7 }

    [telescopes.LAT]
    # Tubes for each LAT should be 19

    # Legacy survey LAT deployed in Chile

    LAT0 = { LFL=8, MFL=54, HFL=23, site="Chile", years=7 }
    LAT1 = { LFL=8, MFL=54, HFL=23, site="Chile", years=7 }

    # Delensing survey LAT deployed at Pole

    LAT2 = { ULFPL=4, LFPL=9, MFPL=54, HFPL=18, site="Pole", years=7 }
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"            
Pole ULFPL1 28
Pole LFL1 0
Pole LFPL1 63
Pole LFL2 0
Pole LFPL2 63
Pole MFPL1 378
Pole MFL1 0
Pole MFL2 0
Pole MFPL2 378
Pole HFL1 0
Pole HFPL1 126
Pole HFL2 0
Pole HFPL2 126
Pole LFS1 14
Pole LFS2 14
Pole MFLS1 42
Pole MFHS1 42
Pole MFLS2 42
Pole MFHS2 42
Pole HFS1 28
Pole HFS2 28
Chile ULFPL1 0
Chile LFL1 112
Chile LFPL1 0
Chile LFL2 112
Chile LFPL2 0
Chile MFPL1 0
Chile MFL1 756
Chile MFL2 756
Chile MFPL2 0
Chile HFL1 322
Chile HFPL1 0
Chile HFL2 322
Chile HFPL2 0
Chile LFS1 0
Chile LFS2 0
Chile MFLS1 0
Chile MFHS1 0
Chile MFLS2 0
Chile MFHS2 0
Chile HFS1 0
Chile HFS2 0

load_noise[source]

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
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

Available input maps

filenames = !ls $base_folder/00000000/*noise*bmap*
import os.path
for f in map(os.path.basename, filenames):
    print(f)
Warning: Python module not loaded, you already have Python loaded via conda init
chile_noise_LAT_HFL1_filtered_telescope_all_time_all_bmap.fits
chile_noise_LAT_HFL2_filtered_telescope_all_time_all_bmap.fits
chile_noise_LAT_LFL1_filtered_telescope_all_time_all_bmap.fits
chile_noise_LAT_LFL2_filtered_telescope_all_time_all_bmap.fits
chile_noise_LAT_MFL1_filtered_telescope_all_time_all_bmap.fits
chile_noise_LAT_MFL2_filtered_telescope_all_time_all_bmap.fits
pole_noise_LAT_HFL1_filtered_telescope_all_time_all_bmap.fits
pole_noise_LAT_HFL2_filtered_telescope_all_time_all_bmap.fits
pole_noise_LAT_LFL1_filtered_telescope_all_time_all_bmap.fits
pole_noise_LAT_LFL1_test_filtered_telescope_all_time_all_bmap.fits
pole_noise_LAT_LFL2_filtered_telescope_all_time_all_bmap.fits
pole_noise_LAT_MFL1_filtered_telescope_all_time_all_bmap.fits
pole_noise_LAT_MFL2_filtered_telescope_all_time_all_bmap.fits
pole_noise_LAT_ULFL1_filtered_telescope_all_time_all_bmap.fits
pole_raster_1h_noise_LAT_MFL1_filtered_telescope_all_time_all_bmap.fits
pole_raster_noise_LAT_MFL1_filtered_telescope_all_time_all_bmap.fits
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))
/global/homes/z/zonca/condajupynersc/lib/python3.7/site-packages/healpy/fitsfunc.py:369: UserWarning: If you are not specifying the input dtype and using the default np.float64 dtype of read_map(), please consider that it will change in a future version to None as to keep the same dtype of the input file: please explicitly set the dtype if it is important to you.
  "If you are not specifying the input dtype and using the default "
/global/homes/z/zonca/condajupynersc/lib/python3.7/site-packages/healpy/fitsfunc.py:391: UserWarning: NSIDE = 512
  warnings.warn("NSIDE = {0:d}".format(nside))
/global/homes/z/zonca/condajupynersc/lib/python3.7/site-packages/healpy/fitsfunc.py:400: UserWarning: ORDERING = NESTED in fits file
  warnings.warn("ORDERING = {0:s} in fits file".format(ordering))
/global/homes/z/zonca/condajupynersc/lib/python3.7/site-packages/healpy/fitsfunc.py:428: UserWarning: INDXSCHM = IMPLICIT
  warnings.warn("INDXSCHM = {0:s}".format(schm))
/global/homes/z/zonca/condajupynersc/lib/python3.7/site-packages/healpy/fitsfunc.py:486: UserWarning: Ordering converted to RING
  warnings.warn("Ordering converted to RING")
input_map[input_map == 0] = hp.UNSEEN
output_map = load_noise(config, site, channel, realization=0)
INFO:root:Base folder: /global/cscratch1/sd/keskital/s4sim/reference_tool_round_2/out
INFO:root:Reading 00000000/pole_noise_SAT_HFS1_telescope_all_time_all_filtered.fits.gz
output_map[np.isnan(output_map)] = hp.UNSEEN
np.testing.assert_allclose(input_map == hp.UNSEEN, output_map == hp.UNSEEN)
simulations_observing_efficiency
{'pole': {'ULFPL1': 0.7094065625,
  'LFPL1': 0.7094065625,
  'LFPL2': 0.7094065625,
  'MFPL1': 0.7094065625,
  'MFPL2': 0.7094065625,
  'HFPL1': 0.7094065625,
  'HFPL2': 0.7094065625,
  'LFS1': 0.5410630787037037,
  'LFS2': 0.5410630787037037,
  'MFLS1': 0.5410630787037037,
  'MFHS1': 0.5410630787037037,
  'MFLS2': 0.5410630787037037,
  'MFHS2': 0.5410630787037037,
  'HFS1': 0.5410630787037037,
  'HFS2': 0.5410630787037037},
 'chile': {'LFL1': 1.0367161458333334,
  'LFL2': 1.0367161458333334,
  'MFL1': 1.0367161458333334,
  'MFL2': 1.0367161458333334,
  'HFL1': 1.0350467156468064,
  'HFL2': 1.0350467156468064}}
observing_efficiency = get_observing_efficiency(
                config["experiment"]["observing_efficiency"], site, telescope, channel
            ) * get_observing_efficiency(
                config["experiment"]["sensitivity_factor"], site, telescope, channel
            )
observing_efficiency
0.08839315576000001
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")
/global/homes/z/zonca/condajupynersc/lib/python3.7/site-packages/healpy/projaxes.py:920: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = copy.copy(mpl.cm.get_cmap("viridis"))
  newcm.set_over(newcm(1.0))
/global/homes/z/zonca/condajupynersc/lib/python3.7/site-packages/healpy/projaxes.py:921: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = copy.copy(mpl.cm.get_cmap("viridis"))
  newcm.set_under(bgcolor)
/global/homes/z/zonca/condajupynersc/lib/python3.7/site-packages/healpy/projaxes.py:922: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = copy.copy(mpl.cm.get_cmap("viridis"))
  newcm.set_bad(badcolor)
/global/homes/z/zonca/condajupynersc/lib/python3.7/site-packages/healpy/projaxes.py:211: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.
  **kwds
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")