<jupyter_start><jupyter_text>Creating a SACC file with the library In this example we will make a Sacc data file using simulated data from CCL, the Core Cosmology Library.<jupyter_code>pylab inline
import sacc
import pyccl as ccl
import datetime
import time<jupyter_output><empty_output><jupyter_text>Getting input dataWe will make some fake window ranges and data vectors from theory predictions using CCL.  We will also need fake N(z) distributions in that case.We will generate 4 bins of galaxy density w(theta) measurements.<jupyter_code>galaxy_density_cl = sacc.standard_types.galaxy_density_cl
nbin = 4<jupyter_output><empty_output><jupyter_text>We will generate 20 ell bins from ell=5 to ell=1000<jupyter_code>ell_edges = np.linspace(5, 1000, 20).astype(int)

# Generate the ranges (ell_min, ell_max)
window_ranges = list(zip(ell_edges[:-1], ell_edges[1:]))
effective_ell = np.array([mean(x) for x in window_ranges])
ell_max = window_ranges[-1][-1]+1
print("Ell window ranges:")
print(window_ranges)<jupyter_output>Ell window ranges:
[(5, 57), (57, 109), (109, 162), (162, 214), (214, 266), (266, 319), (319, 371), (371, 423), (423, 476), (476, 528), (528, 581), (581, 633), (633, 685), (685, 738), (738, 790), (790, 842), (842, 895), (895, 947), (947, 1000)]<jupyter_text>Construct the mock N(z) values and use CCL to calculate the C_ell values in band powers corresponding to it.We use a constant bias within each bin<jupyter_code>cosmo = ccl.Cosmology(Omega_c=0.25, Omega_b=0.05, h=0.72, n_s=0.96, A_s=2.1e-9)

N_z = []
spectra = []
errors = []

# Z sample values
z = np.arange(0., 3.0, 0.01)

for i in range(nbin):
    # Construct n(z)
    mean_z = 0.3*(1+i)
    sigma_z = 0.05
    n_z = np.exp(-0.5*(z - mean_z)**2 / sigma_z**2)
    N_z.append(n_z)
    
    # Use a constant bias in each bin and make a tracer
    bias = np.ones_like(z) * 1.5*(mean_z/0.3)**0.5
    T = ccl.NumberCountsTracer(cosmo, False, (z,N_z[i]), (z, bias))

    # Calculate C_ell over the full ell range (for each integer value)
    ell = np.arange(ell_max+1)
    cl = ccl.angular_cl(cosmo, T, T, ell)

    # Apply window functions to get band powers
    cb = []
    for m,M in window_ranges:
        cb.append(np.mean(cl[m:M]))
    spectra.append(cb)

    # Plot each tomographic bin
    errors.append(np.ones_like(cb)*1e-6)
    line, = loglog(ell, cl)
    errorbar(effective_ell, cb, yerr=errors[-1], fmt='.', color=line.get_color())
    xlim(10,ell_max)
    xlabel("ell")
    ylabel("C_ell")

variances = np.concatenate(errors)**2<jupyter_output><empty_output><jupyter_text>Constructing the data set objectWe now construct the Sacc data set object, using all the different ingredients.We initially start with a blank data set object, and gradually fill it in.<jupyter_code>S = sacc.Sacc()<jupyter_output><empty_output><jupyter_text>Save some metadata.  This is all stored in FITS headers, so can only be simple scalar values<jupyter_code>S.metadata['nbin_lens'] = 4
S.metadata['nbin_source'] = 0 
S.metadata['creator'] = 'Beloved User'
S.metadata['creation'] = datetime.datetime.now().isoformat()
S.metadata['info'] = 'Theory predictions using CCL'<jupyter_output><empty_output><jupyter_text>Now we create some window objects representing the window functions.As an example here we will also create some window functions for mock observations of a pretend new data type.<jupyter_code># In this case we use the same top-hat window for each tomographic bin.
# In reality we would probably use a different window for each
# pair of bins, empirically determined.
windows = [sacc.TopHatWindow(m,M) for m,M in window_ranges]

# Imaginary new statistic's window that I just made up
windows_1pt = []
for i in range(5):
    v = np.arange(i*10., (i+1)*10.)
    w = np.random.uniform(size=10)**2
    windows_1pt.append(sacc.Window(v,w))<jupyter_output><empty_output><jupyter_text>Now we add the tracers and data points<jupyter_code>for i in range(nbin):
    # Add the appropriate tracer
    S.add_tracer('NZ', f'bin_{i}', z, N_z[i])
    # Get the binned spectrum values
    cb = spectra[i]
    # In this case we are just doing LSS, so
    # there is just a single bin, repeated.
    bin_name = f'bin_{i}'

    # Add the values
    S.add_ell_cl(galaxy_density_cl, bin_name, bin_name, effective_ell, cb, window=windows)<jupyter_output><empty_output><jupyter_text>We an specify Sacc covariances in several ways.  If we supply a 1D vector it will be interpreted as the variances of a diagonal covariances<jupyter_code>S.add_covariance(variances)<jupyter_output><empty_output><jupyter_text>We can tag data points with additional information if we want<jupyter_code># Just for fun we add arbitrary tags to a few data points
S.data[4].tags['flag'] = 42
S.data[10].tags['flag'] = 42
# An example cut
cut = S.indices("galaxy_density_cl", ell__gt=900)
print(f"Removing {len(cut)} data points")
S.remove_indices(cut)

# Could also write this to do the cut directly.
# This won't remove anything because we just did the same cut above
S.remove_selection('galaxy_density_cl', ell_gt=900)<jupyter_output>Removing 8 data points<jupyter_text>I/OLet's save to a FITS file now and reload.<jupyter_code>S.save_fits('tmp.fits', overwrite=True)
S2 = sacc.Sacc.load_fits('tmp.fits')<jupyter_output><empty_output><jupyter_text>PlottingLet's plot the original data file and the re-loaded one, to check they agree.<jupyter_code># First the original data set
for b1, b2 in S.get_tracer_combinations('galaxy_density_cl'):
    # We can use this convenience function to pull out the ell values we want
    ell, cl, covmat = S.get_ell_cl('galaxy_density_cl', b1, b2, return_cov=True)
    errorbar(ell, cl, yerr=covmat.diagonal()**0.5, fmt='.', label=f'{b1}-{b2}')
    xscale('log')
    yscale('log')
title("Galaxy C_ell")
legend()

# Now another plot of the second data set that we saved and loaded
for b1, b2 in S2.get_tracer_combinations('galaxy_density_cl'):
    ell, cl, covmat = S2.get_ell_cl('galaxy_density_cl', b1, b2, return_cov=True)
    errorbar(ell, cl, yerr=covmat.diagonal()**0.5, fmt='.', label=f'{b1}-{b2}')
    xscale('log')
    yscale('log')
title("Galaxy C_ell")
legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12358092766785599, 'fraction_numerical': 0.03616607200778463, 'mean_word_length': 4.987378640776699, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 47, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '133341', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2368, 'n_tokens_neox': 2192, 'n_words': 653, 'path': 'examples/Create_Sacc.ipynb', 'repo_name': 'LSSTDESC/sacc', 'size': '6166'}