In [1]:
Copied!
# Run this notebook anywhere: install MeasureIA if it is not already available (e.g. on
# Google Colab). Everything else the notebook needs, including the mock catalogues, comes
# with the package - no simulation data and no clone of the repository are required.
import importlib.util
import subprocess
import sys
if importlib.util.find_spec("measureia") is None:
subprocess.run([sys.executable, "-m", "pip", "install", "-q", "measureia"], check=True)
# Run this notebook anywhere: install MeasureIA if it is not already available (e.g. on
# Google Colab). Everything else the notebook needs, including the mock catalogues, comes
# with the package - no simulation data and no clone of the repository are required.
import importlib.util
import subprocess
import sys
if importlib.util.find_spec("measureia") is None:
subprocess.run([sys.executable, "-m", "pip", "install", "-q", "measureia"], check=True)
In [2]:
Copied!
import matplotlib.pyplot as plt
from measureia import MeasureIABox
from measureia import ReadData
from measureia.mocks import radial_alignment_box_mock
import matplotlib.pyplot as plt
from measureia import MeasureIABox
from measureia import ReadData
from measureia.mocks import radial_alignment_box_mock
MeasureIA usage¶
MeasureIA can be used to measure intrinsic alignment correlation functions and clustering in hydrodynamic simulations. Both the w_gg, w_g+ and multipole moment estimators can be measured.
This notebook runs on a synthetic mock catalogue with a known radial-alignment signal, so it can be executed as-is without any simulation data. Replace the data dictionary below with your own catalogue to measure your own data.
Parameters for MeasureIA object¶
These are the parameters we need to enter to create the MeasureIA object. They govern measurement parameters, e.g. r(p) and pi bins (bounds, number of bins...).
In [3]:
Copied!
simulation = None # Indicator of simulation.
# Choose from [TNG100, TNG100_2, TNG300, EAGLE, HorizonAGN, FLAMINGO_L1,
# FLAMINGO_L2p8, COLIBRE_L400, COLIBRE_L200] for now, which sets the
# boxsize automatically. If your simulation is not included (as for this mock): input None and make sure to add
# the boxsize input parameter in the object creation.
# If this is used, the boxsize is in Mpc/h so the positions and r/rp limits must also be in these units.
boxsize = 205.0 # size of the periodic box, in the same units as the positions. Only used if simulation is None.
snapshot = None # Number of snapshot - if given, the output is saved in group 'Snapshot_[snapshot]'
outfile = './example_IA_box.hdf5' # file path to output data and datafile name (must be hdf5)
num_nodes = 1 # number of CPU cores available for calculation
r_lims = [0.3, 8.0] # r or rp bin edges. Must be in same units as position and boxsize.
num_r = 8 # number of r or rp bins
num_pi = 8 # number of pi bins. (trivial for wg+, not too little for multipoles)
pimax = None # pi max value, if None, half the boxsize is used.
periodicity = True # take periodic boundary conditions into account
simulation = None # Indicator of simulation.
# Choose from [TNG100, TNG100_2, TNG300, EAGLE, HorizonAGN, FLAMINGO_L1,
# FLAMINGO_L2p8, COLIBRE_L400, COLIBRE_L200] for now, which sets the
# boxsize automatically. If your simulation is not included (as for this mock): input None and make sure to add
# the boxsize input parameter in the object creation.
# If this is used, the boxsize is in Mpc/h so the positions and r/rp limits must also be in these units.
boxsize = 205.0 # size of the periodic box, in the same units as the positions. Only used if simulation is None.
snapshot = None # Number of snapshot - if given, the output is saved in group 'Snapshot_[snapshot]'
outfile = './example_IA_box.hdf5' # file path to output data and datafile name (must be hdf5)
num_nodes = 1 # number of CPU cores available for calculation
r_lims = [0.3, 8.0] # r or rp bin edges. Must be in same units as position and boxsize.
num_r = 8 # number of r or rp bins
num_pi = 8 # number of pi bins. (trivial for wg+, not too little for multipoles)
pimax = None # pi max value, if None, half the boxsize is used.
periodicity = True # take periodic boundary conditions into account
Parameters for methods¶
These are the parameters that are needed to configure the methods that are called. They govern what is measured and which methods are used. Check the method docstrings for extensive explanations.
In [4]:
Copied!
tree_path = f"./" # path where tree can be temporarily stored. For large samples this file can grow large.
# If tree_path=None, no trees will be used (slower calculation)
num_jk = 27 # number of jackknife regions. Must be x^3 with x an int. Use 0 for no covariance.
corr_type = "both" # type of correlation to be calculated, choose g+, gg or both
masks = None # optional directory in form of data_dir containing masks to be placed over data in data_dir
rp_cut = None # optional minimum cut on rp for multipoles calculation
tree_path = f"./" # path where tree can be temporarily stored. For large samples this file can grow large.
# If tree_path=None, no trees will be used (slower calculation)
num_jk = 27 # number of jackknife regions. Must be x^3 with x an int. Use 0 for no covariance.
corr_type = "both" # type of correlation to be calculated, choose g+, gg or both
masks = None # optional directory in form of data_dir containing masks to be placed over data in data_dir
rp_cut = None # optional minimum cut on rp for multipoles calculation
Create data dictionary¶
The data input is managed via a dictionary. These need to have specific keys (see below) and contain the position coordinates of both samples, eigen vectors of the chosen axis (shape) and the ratio of the projected axis lengths.
Here the arrays come from the seeded mock catalogue in measureia.mocks: centrals placed uniformly in a periodic box with satellites scattered around them, whose projected major axes point at their own central. That gives a strong, reproducible w_g+ and w_gg signal. Read your own arrays from file instead to measure your own data.
In [5]:
Copied!
mock = radial_alignment_box_mock(n_centrals=600, n_sat=8, boxsize=boxsize)
COM = mock["Position"] # positions of the density sample (centrals + satellites)
COM_shapes = mock["Position_shape_sample"] # positions of the shape sample (satellites)
Semimajor_Axis_Direction = mock["Axis_Direction"]
q = mock["q"]
data_dir = {
"Position": COM, # positions of the position (D) sample
"Position_shape_sample": COM_shapes, # positions of the shape (S) sample
"Axis_Direction": Semimajor_Axis_Direction,
"LOS": 2, # column index of the line of sight parameter (2 assumes the shapes are projected over the z-axis)
"q": q,
}
print(f"{len(COM)} galaxies in the position sample, {len(COM_shapes)} in the shape sample")
mock = radial_alignment_box_mock(n_centrals=600, n_sat=8, boxsize=boxsize)
COM = mock["Position"] # positions of the density sample (centrals + satellites)
COM_shapes = mock["Position_shape_sample"] # positions of the shape sample (satellites)
Semimajor_Axis_Direction = mock["Axis_Direction"]
q = mock["q"]
data_dir = {
"Position": COM, # positions of the position (D) sample
"Position_shape_sample": COM_shapes, # positions of the shape (S) sample
"Axis_Direction": Semimajor_Axis_Direction,
"LOS": 2, # column index of the line of sight parameter (2 assumes the shapes are projected over the z-axis)
"q": q,
}
print(f"{len(COM)} galaxies in the position sample, {len(COM_shapes)} in the shape sample")
5400 galaxies in the position sample, 4800 in the shape sample
Run¶
Providing a dataset name and all parameters defined above, the runs for wgg, wg+ and the multipoles are called bellow.
In [6]:
Copied!
dataset_name = "mock"
# if __name__ == "__main__": # when using multiprocessing, this statement is needed.
MeasureIA_test = MeasureIABox(data_dir, simulation=simulation, snapshot=snapshot, boxsize=boxsize,
separation_limits=r_lims, pi_max=pimax,
num_bins_r=num_r, num_nodes=num_nodes,
num_bins_pi=num_pi, output_file_name=outfile, periodicity=periodicity)
# measure wgg, wg+
MeasureIA_test.measure_xi_w(dataset_name, corr_type, num_jk, temp_file_path=tree_path, masks=masks, chunk_size=100)
# measure multipoles
MeasureIA_test.measure_xi_multipoles(dataset_name, corr_type, num_jk, temp_file_path=tree_path, masks=masks,
rp_cut=rp_cut, chunk_size=100)
dataset_name = "mock"
# if __name__ == "__main__": # when using multiprocessing, this statement is needed.
MeasureIA_test = MeasureIABox(data_dir, simulation=simulation, snapshot=snapshot, boxsize=boxsize,
separation_limits=r_lims, pi_max=pimax,
num_bins_r=num_r, num_nodes=num_nodes,
num_bins_pi=num_pi, output_file_name=outfile, periodicity=periodicity)
# measure wgg, wg+
MeasureIA_test.measure_xi_w(dataset_name, corr_type, num_jk, temp_file_path=tree_path, masks=masks, chunk_size=100)
# measure multipoles
MeasureIA_test.measure_xi_multipoles(dataset_name, corr_type, num_jk, temp_file_path=tree_path, masks=masks,
rp_cut=rp_cut, chunk_size=100)
MeasureIA object initialised with: simulation None that has a periodic boxsize of 205.0 cMpc/h. There are 4800 galaxies in the shape sample and 5400 galaxies in the position sample. The separation bin edges are given by [0.3 0.45223831 0.68173162 1.02768384 1.54919334 2.33534857 3.52044695 5.30693655 8. ] cMpc/h. There are 8 r or r_p bins and 8 pi bins. The maximum pi used for binning is 102.5. The data will be written to ./example_IA_box.hdf5 There are 4800 galaxies in the shape sample and 5400 galaxies in the position sample.
There are 4800 galaxies in the shape sample and 5400 galaxies in the position sample.
Read data¶
There is an internal read class that reads the output file. It provides two ways of reading the data. One, is using the read_MeasureIA_output method, and the second is using the read_cat method. Both are shown below. You can also read the output file manually, using the described output file structure in the README
In [7]:
Copied!
# version 1
IA_data = ReadData(simulation, "example_IA_box", snapshot, data_path="./")
IA_data.read_MeasureIA_output(dataset_name, num_jk)
# version 2
IA_data_wgg = ReadData(simulation, "example_IA_box", snapshot, sub_group="w_gg/", data_path="./")
wgg = IA_data_wgg.read_cat(dataset_name)
rp = IA_data_wgg.read_cat(f"{dataset_name}_rp")
wgg_error = IA_data_wgg.read_cat(f"{dataset_name}_jackknife_{num_jk}")
IA_data_wgp = ReadData(simulation, "example_IA_box", snapshot, sub_group="w_g_plus/", data_path="./")
wgp = IA_data_wgp.read_cat(dataset_name)
rp = IA_data_wgp.read_cat(f"{dataset_name}_rp")
wgp_error = IA_data_wgp.read_cat(f"{dataset_name}_jackknife_{num_jk}")
# version 1
IA_data = ReadData(simulation, "example_IA_box", snapshot, data_path="./")
IA_data.read_MeasureIA_output(dataset_name, num_jk)
# version 2
IA_data_wgg = ReadData(simulation, "example_IA_box", snapshot, sub_group="w_gg/", data_path="./")
wgg = IA_data_wgg.read_cat(dataset_name)
rp = IA_data_wgg.read_cat(f"{dataset_name}_rp")
wgg_error = IA_data_wgg.read_cat(f"{dataset_name}_jackknife_{num_jk}")
IA_data_wgp = ReadData(simulation, "example_IA_box", snapshot, sub_group="w_g_plus/", data_path="./")
wgp = IA_data_wgp.read_cat(dataset_name)
rp = IA_data_wgp.read_cat(f"{dataset_name}_rp")
wgp_error = IA_data_wgp.read_cat(f"{dataset_name}_jackknife_{num_jk}")
Plot w_gg¶
In [8]:
Copied!
plt.figure()
plt.errorbar(IA_data.rp, IA_data.w_gg, yerr=IA_data.errors_w_gg)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r"$w_{gg}$")
plt.xlabel(r"$r_p [\mathrm{cMpc}/h]$")
plt.show()
plt.figure()
plt.errorbar(IA_data.rp, IA_data.w_gg, yerr=IA_data.errors_w_gg)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r"$w_{gg}$")
plt.xlabel(r"$r_p [\mathrm{cMpc}/h]$")
plt.show()
Plot w_g+¶
In [9]:
Copied!
plt.figure()
plt.errorbar(IA_data.rp, IA_data.w_gp, yerr=IA_data.errors_w_gp)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r"$w_{g+}$")
plt.xlabel(r"$r_p [\mathrm{cMpc}/h]$")
plt.show()
plt.figure()
plt.errorbar(IA_data.rp, IA_data.w_gp, yerr=IA_data.errors_w_gp)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r"$w_{g+}$")
plt.xlabel(r"$r_p [\mathrm{cMpc}/h]$")
plt.show()
Plot multipoles clustering¶
In [10]:
Copied!
plt.figure()
plt.errorbar(IA_data.r, IA_data.multipoles_gg, yerr=IA_data.errors_multipoles_gg)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r"$\tilde{\xi}_{gg,0}$")
plt.xlabel(r"$r [\mathrm{cMpc}/h]$")
plt.show()
plt.figure()
plt.errorbar(IA_data.r, IA_data.multipoles_gg, yerr=IA_data.errors_multipoles_gg)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r"$\tilde{\xi}_{gg,0}$")
plt.xlabel(r"$r [\mathrm{cMpc}/h]$")
plt.show()
Plot multipoles g+¶
In [11]:
Copied!
plt.figure()
plt.errorbar(IA_data.r, IA_data.multipoles_gp, yerr=IA_data.errors_multipoles_gp)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r"$\tilde{\xi}_{g+,2}$")
plt.xlabel(r"$r [\mathrm{cMpc}/h]$")
plt.show()
plt.figure()
plt.errorbar(IA_data.r, IA_data.multipoles_gp, yerr=IA_data.errors_multipoles_gp)
plt.xscale('log')
plt.yscale('log')
plt.ylabel(r"$\tilde{\xi}_{g+,2}$")
plt.xlabel(r"$r [\mathrm{cMpc}/h]$")
plt.show()
