As we saw in the last blog post, the BOSS dataset contains about 1 million galaxies and their distribution in the Universe. The position of these galaxies is what carries the cosmological information. Imagine that there would be more dark matter in the Universe. The additional matter would have additional gravitational force, which would clump together the material. On the other hand, with less dark matter, the galaxies would be more spread out.
We can now compare these distributions with our data from the BOSS survey, and depending on which of the two distributions looks more like the data, we can determine how much dark matter there is in the Universe.
In practice, we do not want to compare the actual distribution of galaxies, but instead we compare a summary statistic, meaning a compression of these galaxies, which carries all the information we need. There are many choices for such a summary statistic, but here we will use the power spectrum.
To calculate the power spectrum, we first need to download two more catalogues
wget -N https://data.sdss.org/sas/dr12/boss/lss/random0_DR12v5_CMASSLOWZTOT_North.fits.gz -P path/to/folder/ wget -N https://data.sdss.org/sas/dr12/boss/lss/random0_DR12v5_CMASSLOWZTOT_South.fits.gz -P path/to/folder/
To read this file, you can use the read_data() function we defined earlier (the file is significantly larger and it can take a while to read it into memory). Now we can combine the data catalogues and random catalogue, assign the point distributions to a 3D grid and calculate the power spectrum, all using triumvirate.
from astropy import units as u from astropy.coordinates import Distance from astropy.cosmology import Planck15 from triumvirate.parameters import fetch_paramset_template from triumvirate.catalogue import ParticleCatalogue from triumvirate.twopt import compute_powspec def measure_pk(cat_data, cat_random, tag=''): """ Calculate the power spectrum monopole P0(k) using Triumvirate. Assumes cat_data and cat_random were read with read_data() that: - converts columns to native-endian - promotes RA/DEC/Z/NZ/WEIGHT_* to float64 """ # Parameters paramset = fetch_paramset_template('dict') paramset['boxsize']['x'] = 5000.0 paramset['boxsize']['y'] = 10000.0 paramset['boxsize']['z'] = 5000.0 paramset['ngrid']['x'] = 512 paramset['ngrid']['y'] = 1024 paramset['ngrid']['z'] = 512 paramset['tags']['output'] = tag paramset['degrees']['ELL'] = 0 paramset.update({ 'range': [0.005, 0.295], 'num_bins': 29, 'norm_convention': 'particle', 'statistic_type': 'powspec', }) # Data catalogue: spherical -> Cartesian (Mpc) r1_mpc = Planck15.comoving_distance(cat_data['Z']).to_value(u.Mpc) r1_mpc_h = r1_mpc * Planck15.h # now in Mpc/h ra1 = np.deg2rad(cat_data['RA']) dec1 = np.deg2rad(cat_data['DEC']) x1 = r1_mpc_h * np.cos(dec1) * np.cos(ra1) y1 = r1_mpc_h * np.cos(dec1) * np.sin(ra1) z1 = r1_mpc_h * np.sin(dec1) ws1 = cat_data['WEIGHT_SYSTOT'] * (cat_data['WEIGHT_NOZ'] + cat_data['WEIGHT_CP'] - 1.0) wc1 = cat_data['WEIGHT_FKP'] nz1 = cat_data['NZ'] catalogue_data = ParticleCatalogue(x1, y1, z1, nz=nz1, ws=ws1, wc=wc1) # Random catalogue r2_mpc = Planck15.comoving_distance(cat_random['Z']).to_value(u.Mpc) r2_mpc_h = r2_mpc * Planck15.h # now in Mpc/h ra2 = np.deg2rad(cat_random['RA']) dec2 = np.deg2rad(cat_random['DEC']) x2 = r2_mpc_h * np.cos(dec2) * np.cos(ra2) y2 = r2_mpc_h * np.cos(dec2) * np.sin(ra2) z2 = r2_mpc_h * np.sin(dec2) wc2 = cat_random['WEIGHT_FKP'] nz2 = cat_random['NZ'] catalogue_rand = ParticleCatalogue(x2, y2, z2, nz=nz2, wc=wc2) # Compute power spectrum (monopole) results = compute_powspec( catalogue_data, catalogue_rand, degree=0, paramset=paramset, save='.txt' # or a filename/path you prefer ) return results
Here we use a grid with 512 grid points in two dimensions and 1024 in one dimension. There are also two weightings included in this calculation, the 'Weight' and 'WEIGHT_FKP' columns. The first weight refers to a correction of incompleteness in the data catalogue, while the second weight tries to optimise the signal-to-noise.
If we use this code to measure the BOSS power spectrum for the northern and southern part we get
 |
Figure 3: Power spectrum measurements for the northern (orange) and southern (blue) parts of the BOSS dataset. |
Here we see that the power spectra of the northern and southern parts are slightly different. Whether these differences are significant is not clear though, since we don't have uncertainties on these measurements (yet).
To extract cosmological information from this power spectrum we need to get a model. We use the classy package
from classy import Class def get_Pk(Omega_m = 0.3): #Start by specifying the cosmology Omega_b = 0.05 Omega_cdm = Omega_m - Omega_b h = 0.7 #H0/100 A_s = 2.1e-9 n_s = 0.96 #Create a params dictionary #Need to specify the max wavenumber k_max = 10 #UNITS: 1/Mpc params = { 'output':'mPk', 'non linear':'halofit', 'Omega_b':Omega_b, 'Omega_cdm':Omega_cdm, 'h':h, 'A_s':A_s, 'n_s':n_s, 'P_k_max_1/Mpc':k_max, 'z_max_pk':10. #Default value is 10 } #Initialize the cosmology and compute everything cosmo = Class() cosmo.set(params) cosmo.compute() #Specify k and z k = np.logspace(-5, np.log10(k_max), num=1000) #Mpc^-1 z = 1. #Call these for the nonlinear and linear matter power spectra Pnonlin = np.array([cosmo.pk(ki, z) for ki in k]) Plin = np.array([cosmo.pk_lin(ki, z) for ki in k]) #NOTE: You will need to convert these to h/Mpc and (Mpc/h)^3 #to use in the toolkit. To do this you would do: k /= h Plin *= h**3 Pnonlin *= h**3 return k, Plin
 |
Figure 4: Power spectrum models with different amounts of dark matter using the class Boltzmann code. |
However, before we can go ahead and constrain cosmological parameters we have to get an estimate of the uncertainty on the measurements. That will be the subject of the next post.
You can find the code for this project on GitHub.
cheers
Florian
No comments:
Post a Comment