mirror of
https://github.com/Richard-Sti/csiborgtools.git
synced 2024-12-23 01:18:02 +00:00
ee222cd010
* Update nb * Update script * Update script * Rename * Update script * Update script * Remove warning * Ignore minors when extracting MAH * Fix paths bug * Move notebooks * Move files * Rename and delete things * Rename file * Move file * Rename things * Remove old print statement * Add basic MAH plot * Add random MAH path * Output snapshot numbers * Add MAH random extraction * Fix redshift bug * Edit script * Add extracting random MAH * Little updates * Add CB2 redshift * Add some caching * Add diagnostic plots * Add caching * Minor updates * Update nb * Update notebook * Update script * Add Sorce randoms * Add CB2 varysmall * Update nb * Update nb * Update nb * Use catalogue HMF * Move definition of radec2galactic * Update nb * Update import * Update import * Add galatic coords to catalogues * Update nb
55 KiB
55 KiB
Matching of haloes to clusters¶
In [49]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.colors import Normalize
from scipy import stats
from match_clusters import *
import scienceplots
%load_ext autoreload
%autoreload 2
%matplotlib inline
Plot $R - \log M$ as a function of significance¶
In [79]:
cat = csiborgtools.read.CSiBORG2Catalogue(
17417, 99, "main", bounds={"totmass": (1e12, None)})
pos = cat["cartesian_pos"]
logmass = np.log10(cat["totmass"])
model = csiborgtools.match.MatchingProbability(pos, logmass)
In [78]:
xrange = np.linspace(13, 15.2, 250)
log_p_values = np.linspace(-3, -0.25, 1000)
cmap = cm.viridis
norm = Normalize(vmin=log_p_values.min(), vmax=log_p_values.max())
sm = cm.ScalarMappable(cmap=cmap, norm=norm)
sm.set_array([]) # Necessary for the colorbar to use the ScalarMappable
fig, ax = plt.figure(), plt.gca()
for p in log_p_values:
y = [model.inverse_cdf(10**p, x) for x in xrange]
ax.plot(xrange, y, color=cmap(norm(p)))
ls = [":", "--", "-."]
for n in range(1, 4):
p = stats.norm.sf(n) * 2
y = [model.inverse_cdf(p, x) for x in xrange]
ax.plot(xrange, y, color="red", label=fr"${n}\sigma$", linestyle=ls[n-1])
ax.set_xlabel(r"$\log M_{\mathrm{FoF}} ~ [M_\odot / h]$")
ax.set_ylabel(r"$R ~ [\mathrm{Mpc} / h]$")
ax.set_ylim(1, 55)
ax.legend()
cbar = fig.colorbar(sm, ax=ax, label=r"$\log p$", pad=0)
ax.set_xlim(xrange.min(), xrange.max())
fig.tight_layout()
fig.savefig("../../plots/matching_probability.png", dpi=450)
fig.show()
Test matching¶
In [2]:
simname = "csiborg2_main"
bounds = {"totmass": (1e12, None)}
cats = open_cats(simname, bounds)
In [3]:
model = csiborgtools.match.MatchCatalogues(cats)
In [60]:
cluster = csiborgtools.clusters["Leo"]
x0 = cluster.cartesian_pos(cats[0].boxsize)
logmass0 = np.log10(cluster.mass)
print(logmass0)
In [61]:
ps, indxs = model(x0, logmass0, pvalue_threshold=0.05, rmax=30)
print(ps)
In [29]:
# ps_coma = np.array(list(ps.values()))
# ps_centaurus = np.array(list(ps.values()))
# ps_shapley = np.array(list(ps.values()))
ps_virgo = np.array(list(ps.values()))
In [57]:
ps = {"Coma": ps_coma, "Centaurus": ps_centaurus, "Shapley": ps_shapley, "Virgo": ps_virgo}
with plt.style.context('science'):
plt.figure()
for name, p in ps.items():
plt.hist(np.log10(p), bins="auto", label=name, density=True, alpha=0.5)
plt.legend(loc="upper left", ncol=2)
plt.xlabel(r"$\log p$")
plt.ylabel("Counts")
plt.tight_layout()
plt.savefig("../../plots/matching_probability_histogram.png", dpi=450, bbox_inches="tight")
plt.show()
In [6]:
dx = np.asarray([np.linalg.norm(cats[i]["cartesian_pos"][indx] - x0) for i, indx in indxs.items()])
logM = np.asarray([np.log10(cats[i]["totmass"][indx]) for i, indx in indxs.items()])
In [7]:
ps
Out[7]:
In [114]:
plt.figure()
plt.scatter(np.arange(len(dx)), dx)
plt.show()
In [ ]:
In [98]:
ps, dlog, indxs = model.cdf_per_halo(x0, logmass0)
r = np.linalg.norm(pos - x0.reshape(1, -1), axis=1)
In [99]:
p, k = model.match_halo(x0, logmass0, pvalue_threshold=0.1, max_absdlogmass=1)
print(p, k, logmass[k])
In [ ]:
In [111]:
plt.figure()
plt.plot(xrange, y)
plt.plot(xrange, y2)
plt.show()
In [71]:
x, y = model.pdf_marginalized(xrange, 14.0, 42.5, 4)