csiborgtools/notebooks/MAH/mah.ipynb

900 lines
390 KiB
Plaintext
Raw Normal View History

{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Constraining MAH in $\\texttt{BORG}$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.cm import ScalarMappable\n",
"from matplotlib.colors import Normalize\n",
"from cache_to_disk import delete_disk_caches_for_function\n",
"from scipy.stats import norm\n",
"from scipy.signal import savgol_filter\n",
"\n",
"from mah import *\n",
"\n",
"%load_ext autoreload\n",
"%autoreload 2\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## $\\texttt{MDPL2}$ to $\\texttt{CB2}$ MAH comparison"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHVCAYAAAB8NLYkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOydd3wkdf3/n5+Z7T39Uq/3CgdH7yAcIEhHURFRFMGGP1Es2L8WUBE8EQVBQZAmAtKk93pwvff0nu1lyuf3xyR7ud42l+Runo9HHsnOzs58dpPsvuZdXm8hpZTY2NjY2NjY2NgMe5TBXoCNjY2NjY2NjU1hsIWdjY2NjY2Njc0Bgi3sbGxsbGxsbGwOEGxhZ2NjY2NjY2NzgGALOxsbGxsbGxubAwRb2NnY2NjY2NjYHCDYws7GxsbGxsbG5gDBFnY2NjY2NjY2NgcItrCzsbGxsbGxsTlAsIWdjY2NjY2Njc0Bgi3sbGxsbGxsbGwOEA46Yfff//6XiRMnMn78eO68887BXo6NjY2NjY2NTcEQUko52IvYX+i6zpQpU3j55ZcJh8PMnj2bt956i5KSksFemo2NjY2NjY3NPnNQRezee+89pk6dSnV1NYFAgLlz5/K///1vsJdlY2NjY2NjY1MQHIO9gD3htdde46abbmL+/Pk0Nzfz2GOP8YlPfGKLfebNm8dNN91ES0sLM2fO5LbbbmPOnDkANDU1UV1dnd+3urqaxsbG3T6/aZo0NTURDAYRQhTkOdnY2NjY2NjY7AwpJfF4nKqqKhRl5zG5YSXskskkM2fO5POf/zznn3/+Nvc/+OCDXHfddfz5z3/miCOO4JZbbuH0009n5cqVlJeX7/H5stks2Ww2f7uxsZEpU6bs03OwsbGxsbGxsdkb6uvrqamp2ek+w0rYzZ07l7lz5+7w/t/97nd88Ytf5IorrgDgz3/+M0899RR/+9vf+O53v0tVVdUWEbrGxsZ8NG97/PKXv+QnP/nJNtvr6+sJhUL78ExsbGxsdoKUYJpIRSGrm2QXL4MXX8BIpzESKYxMFtIpZCYLmQyN51xI5yFHYkiT8KIPmfK7n6PmsijZLGoug5rLoeayqFqOBV+7gXXnfRopJaULP+CEb35uh8tYftU3WfvpqwAIrVzK8V+8cIf7rv70l1h51TcA8NVv4OTLdvxevfbiz7H82u8A4G5r5tirLtnhvo1zP8GKL10HgDPWw/GfPaf/C5V/uSSw/thTee3Kb6PpJmY6xZVfPb/3fklfMXlfVfniKUfw90u/hYlESvjdjy7e8sT9qs+XjpvFnz75nfzt3//qClxaZrvrXVM3md9ffmP+9q9/+2UCqWj+tooJgIHCpsox/PoLv2Bu5r9UGk1Mums+/mgMAAWJR2YwUcgINy1l1fzsyzdzfvphRutrKb6vgUhHp/W6oBM04xhCJSpCdIdL+cHXbuO7sZ8zU1/IhkfLCLV0957fIGzGMIRCVIRJeoNc///+wnXxmzhce491T1QQqe/MrzVsRjER9CgRNIeLb9zwd74S/wPHaW+w9ulKita3AyCQFJk9SKBLFAHwyevu5ir9fo4336X1WR+Va1qQEgTgIYOJIIsbgPO+cBtjXI3U0sbRL73PnBWL878INzkAMr37XvrZm6jztzNebGLaqys4aumC3jWAkxwg0HAgEVz+qV+ghxxESHDW26/w8YWvomDmvwQSiSCLi6su+jGiBM5S3uaE99+j+v1GRO9zs/a1iOLnK+f9gFhFiMuVp7lkwf/grc1BoP70EOCbZ3+HtppSPqs8wyVLn4NXt7+vgYJ6lhvG9Eqz5Rq8kMHMQUKDWiAYDG73sf0Zts0TQogtUrG5XA6fz8cjjzyyRXr28ssvp6enh8cffxxd15k8eTKvvPLKbjVPbB2xi8Vi1NbWEo1GbWFnY2OzbzQ2wuuvQ2MjZkMDRkMjNDQgmppQWppZ99vb2XjKWWQ0g5L/PcWR375qh4da+sNf0XTp51AElLz3BrM+t2MBtuG7P6b5C9cC4F+ygMlXXIzp9mC63Ui3B9Nj/Wy6vbSffwkd51jHcrY2U/XXP265j8eL6XJjuj2kJk0hNXkaACKbxbd6Rf641v7W46TTBf1SSVJKDFOSM0xyuolm9P/Z+p4zrJ+zmklaM0jnDJI5naxmoJvW4w1TopuSLT/RJIoQqIr1JQQoovc71vf+24QQKPR+7729t0S0VuoyK1nnnU7CYQmdafE3+ELjD9jgmcwto27P7/v1jdcwOr2UO6t/zpLgsQCMS37EtfXfpNk1il+PuSe/79WbrmNi6kPurfwe88MfA6Ams5L/t+FLtLlq+b8x9+b3vbzxJ4xPzeffFV/nw9Ap1rpyzVza/Bu6lRL+UvodNMPEMCVzkq9SqdfzvvMwVqnjMExJQO/hzOzTpKWTe9VPoBvW63yosYgaWljEBFbJOiSgSo0prCMpPaymdovXzvqdSIQQllASlgijd5++Pfsekn9k7/75bVvf3+8EAVIUE81/lcgoz4pjiIsAAFeY/+bL8kEcvcJ6az6n/JzFYiIAnzL/y7fk37e7H8A1yg94R8wE4DzzBX4g78BAkMFNCi+dhOkmTLcIcb84i+ViLBWyg1PNt5jJSirpoJpW1lDHD5WvkcFFynTwpPwaZSKKKUERcDfn8IE+m1/99ofITIIw7Jb+GFYRu53R0dGBYRhUVFRssb2iooIVK1YA4HA4+O1vf8tJJ52EaZpcf/31O+2IdbvduN3uAV23jY3N8ENKuYUAyeqGJUB0k1QmB/UNeNauxrNuNZ61q3E31ONsbab+G9+l+9QzkBKK/vcCk6/5HGB1sW1dNZNat5H0cQZuh4oydRrtH7/AEkZ5YeXJCyd9zpEE3NbbuTZ1BsvverDfPu4tfjYCm6/4k9Nm8cH7q3brOWsVlWz8wS927/Vxu4lNmUFaM8hqBmnNIJM1ySSyZLQU6ZxBPKuRyBpkcgaaaYkL05QYcrNQQwDSElh94qBPpDl6vztVBa9z8zaxj4JsZ8zueZYirZXlgaNo9E4AIKB3c2zXo2jCzYtln8nv+/nGG6nLrOQf1T9iYfhkANLOCHE1QlYNoCqb1/hB5AxWBebQ5h2T3x5zV/B86WdIquEt9n2z5AIWhk+m3j8VISCnm6wSo7iy9iky0oEezaAZJroh+Z7jOnI+k1y3idbRaW03Hfxb3oAEZKIrH+1cKGeBmNV7lqT1e8TFEj7RG7XS8sLrXTGN95iGEODIv94uVolJCCCwD69/UCaooYUy2UU5XQRJ4pMZPGTwyix/Uy6gSVif8xeYz3GF+W+KiOJC3+ZYq5VRLBKTAEgJHw5poqPSSgkdFJEWbjJYX0kRQOld9wJlCrfIz5LBTZrN+2R6999IVX7f/yon8RQnouPYrEz7kJJz5Et8w7yXw+TSbdZXJVs5XX+FEhHl1+LzfNPxPdK+as4QbxP31rLYPYsL7/0dwUyCjkmTYcXy3XoNDxhht7ucc845nHPOObve0cbG5qAnoxkksjqprEHOMEjnTBJZjVTOQI8ncK9fRyJcRLKkAgGMePsVjr7+KhzZ7afptFWraZtzEgC50irKDj2CTHkluYoR6JVVGJVV6JVV5Coq0coqKHG5ADAmTGDN7+/YrTXrRcX0nHBKQZ7/jpBSktFMkjmdRFYnmdVJZg2SOY1ExhJt6ZyBYVoCWDdNTLMvlGZ9+DnUzeKsT6CpToGqgCr6Imz7sUlNyvwHc3GumVM67kWVOv+q/l5+l5mxV5iaeJuEoygv7HxGlNM67iWphrcQdht9U1EwkWKzZN/gm86PJz6+zanfKdr8mdR30bBRVrDG/1nr4qErlb+IWKBPI6MZZNtNNKMdKSWm3DLdvPVT6ouSKb0RMKVvm6L0u3/gXuuIjFJFGxEZJ0yCCDHCMk6YOBFi3KJcTqsoA6xo2RfkIzs81pOcRBOWsFMwqaAzf18CL1290bIuESaHM3/fs+I4XlaPoIswplB3ut4VYiwrxNjdem6G2FJGKdLYfHwhOMd4iZmsBOBDOZHFjGW1rGOtUofPCXfp38NAYVndZ8n4j0AIwTJZQ0s
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"min_age = 2.8\n",
"pmin, pmax = 5, 95\n",
"bounds = [14.5, 14.7]\n",
"mass_norm = 10**np.mean(bounds)\n",
"alphas = {3: 0.25, 2: 0.5, 1: 0.75}\n",
"xrange_cb2, random_mah_cb2 = extract_mah(\"csiborg2_random\", bounds, \"MainProgenitorMass\", min_age=2.8)\n",
"xrange_mdpl2, random_mah_mdpl2 = extract_mah_mdpl2(bounds, min_age=2.8)\n",
"\n",
"random_mah_cb2 /= random_mah_cb2[:, -1].reshape(-1, 1)\n",
"random_mah_mdpl2 /= random_mah_mdpl2[:, -1].reshape(-1, 1)\n",
"# random_mah_cb2 /= mass_norm\n",
"# random_mah_mdpl2 /= mass_norm\n",
"\n",
"\n",
"\n",
"plt.figure()\n",
"cols = plt.rcParams[\"axes.prop_cycle\"].by_key()[\"color\"]\n",
"for n in [3, 2, 1]:\n",
" pmin, pmax = norm.cdf(-n) * 100 , norm.cdf(n) * 100\n",
" ylow, yhigh = np.percentile(random_mah_mdpl2, [pmin, pmax], axis=0)\n",
" plt.fill_between(xrange_mdpl2, ylow, yhigh, alpha=alphas[n], color=cols[0],\n",
" facecolor=cols[0], label=f\"MDPL2 Random ({len(random_mah_mdpl2)})\" if n == 1 else None)\n",
"\n",
"\n",
"pmin, pmax = norm.cdf(-1) * 100 , norm.cdf(1) * 100\n",
"ylow, yhigh = np.percentile(random_mah_cb2, [pmin, pmax], axis=0)\n",
"plt.plot(xrange_cb2, ylow, c=cols[1], ls=\"--\", label=rf\"$1\\sigma$ CB2 Random ({len(random_mah_cb2)})\", zorder=1)\n",
"plt.plot(xrange_cb2, yhigh, c=cols[1], ls=\"--\", zorder=1)\n",
"\n",
"pmin, pmax = norm.cdf(-2) * 100 , norm.cdf(2) * 100\n",
"ylow, yhigh = np.percentile(random_mah_cb2, [pmin, pmax], axis=0)\n",
"plt.plot(xrange_cb2, ylow, c=cols[1], ls=\"dotted\", label=rf\"$2\\sigma$ CB2 Random ({len(random_mah_cb2)})\", zorder=1)\n",
"plt.plot(xrange_cb2, yhigh, c=cols[1], ls=\"dotted\", zorder=1)\n",
"\n",
"# plt.fill_between(xrange_mdpl2,\n",
"# np.ones_like(xrange_mdpl2) * 10**bounds[0] / mass_norm,\n",
"# np.ones_like(xrange_mdpl2) * 10**bounds[1] / mass_norm, alpha=0.25,\n",
"# label=r\"$z = 0$ selection\", color=\"gray\", zorder=0,\n",
"# hatch=\"///\")\n",
"\n",
"# They appear to be plotting the min-max\n",
"plt.plot(RANDOM_MAH_Sorce_Virgo_UPPER[:, 0],\n",
" RANDOM_MAH_Sorce_Virgo_UPPER[:, 1], c=\"red\", ls=\"--\", label=\"Sorce+2016 Virgo randoms\")\n",
"plt.plot(RANDOM_MAH_SORCE_Virgo_LOWER[:, 0],\n",
" RANDOM_MAH_SORCE_Virgo_LOWER[:, 1], c=\"red\", ls=\"--\")\n",
"\n",
"\n",
"plt.yscale(\"log\")\n",
"plt.legend(loc=\"lower right\")\n",
"plt.xlabel(r\"$\\mathrm{Age} ~ [\\mathrm{Gyr}]$\")\n",
"plt.ylabel(r\"$M_{\\rm main}(t) / M_{\\rm main}(z = 0)$\")\n",
"plt.xlim(xrange_mdpl2.min(), xrange_mdpl2.max())\n",
"plt.ylim(1e-2, 1.01)\n",
"\n",
"plt.tight_layout()\n",
"plt.savefig(\"../../plots/MAH_MDPL2_CB2_comparison_VIRGO.png\", dpi=450)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## MAH of specific clusters in $\\texttt{CB}$\n",
"\n",
"### Summary plot"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"nsim0 = 17417\n",
"simname = \"csiborg2_main\"\n",
"min_logmass = 12.25\n",
"delete_cache = False\n",
"\n",
"if delete_cache:\n",
" delete_disk_caches_for_function(\"load_data\")\n",
"cat0, catxs, merger_trees, overlaps = load_data(nsim0, simname, min_logmass)\n",
"nsimxs = [cat.nsim for cat in catxs]"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"catx = cat0"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"pos = catx[\"spherical_pos\"]\n",
"\n",
"ra, dec = pos[:, 1], pos[:, 2]\n",
"l, b = csiborgtools.flow.radec_to_galactic(pos[:, 1], pos[:, 2])"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"70.92198581560284"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"50 / 0.705"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"79.66499999999999"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"113 * 0.705"
]
},
{
"cell_type": "code",
"execution_count": 61,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGdCAYAAAA8F1jjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAu1klEQVR4nO3de3xU5Z3H8e/kMgMIMwmQawkQkIIoggaM8baupAQLrqysl0oVlBesNrpFKEpAQVwlFazFO+pWoLtYqVsvLRWQgmBdIyiKikoAxYYSJuCFGQTJ9dk/pkwykECiTM4z4fN+vc7r9cw5T8bfPJ6Z+XKec864jDFGAAAAFopzugAAAICmEFQAAIC1CCoAAMBaBBUAAGAtggoAALAWQQUAAFiLoAIAAKxFUAEAANZKcLqA76uurk7l5eXq1KmTXC6X0+UAAIBmMMZo//79yszMVFxc08dNYj6olJeXKysry+kyAADAd7Bz505169atye0xH1Q6deokKfRCvV6vw9UAaLGqKulXvwq1p0yR3G5n6wHQKoLBoLKyssLf402J+aByeLrH6/USVIBYVFUleTyhttdLUAFOMsc7bYOTaQEAgLVi/ogKgBgXFycNGlTfBoAGCCoAnJWQII0a5XQVACzFP18AAIC1OKICwFnGSNXVoXZiosT9kAA0wBEVAM6qrpbmzAkthwMLAPwDQQUAAFiLoAIAAKxFUAEAANYiqAAAAGsRVAAAgLUIKgAAwFrcRwWAs+LipP7969sA0EBUPxVqa2t11113KTs7W+3bt1fv3r31n//5nzLGhPsYYzRz5kxlZGSoffv2ys/P17Zt26JZFgCbJCRIV10VWhL4txOASFENKvfff7+eeOIJPfroo/rkk090//33a+7cuXrkkUfCfebOnauHH35YCxYs0Pr163XKKaeooKBAhw4dimZpAAAgBrhMw8MbJ9jIkSOVlpam3/zmN+F1o0ePVvv27fU///M/MsYoMzNTU6ZM0S9+8QtJUiAQUFpamhYtWqRrrrnmuP+NYDAon8+nQCAgr9cbrZcCAABOoOZ+f0f1iMp5552n1atXa+vWrZKk999/X2+88YYuvfRSSdKOHTvk9/uVn58f/hufz6fc3FyVlJQ0+pyVlZUKBoMRC4AYVlUl3X13aKmqcroaAJaJ6oTwtGnTFAwG1a9fP8XHx6u2tlb33XefxowZI0ny+/2SpLS0tIi/S0tLC287UnFxsWbPnh3NsgEAgCWiekTl97//vZYsWaJnn31W7777rhYvXqwHHnhAixcv/s7PWVRUpEAgEF527tx5AisGAAA2ieoRlalTp2ratGnhc00GDBigv/3tbyouLtbYsWOVnp4uSaqoqFBGRkb47yoqKjRo0KBGn9Pj8cjj8USzbAAAYImoHlE5ePCg4o64L0J8fLzq6uokSdnZ2UpPT9fq1avD24PBoNavX6+8vLxolgYAAGJAVI+oXHbZZbrvvvvUvXt3nX766Xrvvff04IMP6sYbb5QkuVwuTZo0Sffee6/69Omj7Oxs3XXXXcrMzNSoUaOiWRoAAIgBUQ0qjzzyiO666y797Gc/0549e5SZmal///d/18yZM8N9br/9dh04cEATJ07Uvn37dMEFF2jFihVq165dNEsDAAAxIKr3UWkN3EcFiHE1NdLSpaH21Vdzd1rgJNHc728+EQA4KyFB+sctCwDgSPwCGAAAsBZBBQAAWIupHwDOqqqS5s0LtadOldxuZ+sBYBWCCgDnVVc7XQEASzH1AwAArEVQAQAA1iKoAAAAaxFUAACAtQgqAADAWlz1A8BZLpfUs2d9GwAaIKgAcFZiojRunNNVALAUUz8AAMBaBBUAAGAtpn4AOKuqSpo/P9SeNIlb6AOIQFAB4LyDB52uAIClmPoBAADWIqgAAABrEVQAAIC1CCoAAMBaBBUAAGAtrvoB4CyXS8rMrG8DQAMEFQDOSkyUJk50ugoAlmLqBwAAWIugAgAArMXUDwBnVVdLjz0WahcWhqaCAOAfCCoAnGWMtG9ffRsAGmDqBwAAWIugAgAArEVQAQAA1op6UNm1a5d++tOfqkuXLmrfvr0GDBigd955J7zdGKOZM2cqIyND7du3V35+vrZt2xbtsgAAQAyIalD5+uuvdf755ysxMVHLly/Xxx9/rF/96ldKTk4O95k7d64efvhhLViwQOvXr9cpp5yigoICHTp0KJqlAQCAGBDVq37uv/9+ZWVlaeHCheF12dnZ4bYxRvPnz9edd96pyy+/XJL029/+VmlpaXrppZd0zTXXRLM8ADZwuaSUlPo2ADQQ1SMqf/zjHzV48GBdeeWVSk1N1VlnnaWnn346vH3Hjh3y+/3Kz88Pr/P5fMrNzVVJSUmjz1lZWalgMBixAIhhiYmh+6dwDxUAjYhqUPnss8/0xBNPqE+fPlq5cqVuvvlm/cd//IcWL14sSfL7/ZKktLS0iL9LS0sLbztScXGxfD5feMnKyormSwAAAA6KalCpq6vT2WefrTlz5uiss87SxIkTNWHCBC1YsOA7P2dRUZECgUB42blz5wmsGAAA2CSqQSUjI0P9+/ePWHfaaaeprKxMkpSeni5JqqioiOhTUVER3nYkj8cjr9cbsQCIYYdvof/YY6E2ADQQ1aBy/vnnq7S0NGLd1q1b1aNHD0mhE2vT09O1evXq8PZgMKj169crLy8vmqUBsIUx0t69oYVb6AM4QlSv+rntttt03nnnac6cObrqqqu0YcMGPfXUU3rqqackSS6XS5MmTdK9996rPn36KDs7W3fddZcyMzM1atSoaJYGAABiQFSDypAhQ/Tiiy+qqKhI99xzj7KzszV//nyNGTMm3Of222/XgQMHNHHiRO3bt08XXHCBVqxYoXbt2kWzNAAAEANcxsT2sdZgMCifz6dAIMD5KkAsqqqS5swJtadPl9xuZ+sB0Cqa+/3Nb/0AAABrEVQAAIC1onqOCgAcl8slJSXVtwGgAYIKAGclJkqTJjldBQBLMfUDAACsRVABAADWYuoHgLOqq6WFC0PtG27gF5QBRCCoAHCWMVJ5eX0bABpg6gcAAFiLoAIAAKxFUAEAANYiqAAAAGsRVAAAgLW46geA8zp0cLoCAJYiqABwltst3X6701UAsBRTPwAAwFoEFQAAYC2mfgA4q7paWrIk1B4zhlvoA4hAUAHgLGOkzz+vbwNAA0z9AAAAaxFUAACAtQgqAADAWgQVAABgLYIKAACwFlf9AHAelyQDaAJBBYCz3G5pxgynqwBgKaZ+AACAtQgqAADAWkz9AHBWTY20dGmoffXVUgIfSwDq8YkAwFl1ddK2bfVtAGig1aZ+fvnLX8rlcmnSpEnhdYcOHVJhYaG6dOmijh07avTo0aqoqGitkgAAgOVaJai8/fbbevLJJ3XmmWdGrL/tttv0pz/9Sc8//7zWrVun8vJyXXHFFa1REgAAiAFRDyrffPONxowZo6efflrJycnh9YFAQL/5zW/04IMP6pJLLlFOTo4WLlyoN998U2+99Va0ywIAADEg6kGlsLBQI0aMUH5+fsT6jRs3qrq6OmJ9v3791L17d5WUlDT5fJWVlQoGgxELAABom6J6Mu1zzz2nd999V2+//fZR2/x+v9xut5KSkiLWp6Wlye/3N/mcxcXFmj179okuFQAAWChqR1R27typn//851qyZInatWt3wp63qKhIgUAgvOzcufOEPTcAALBL1I6obNy4UXv27NHZZ58dXldbW6vXX39djz76qFauXKmqqirt27cv4qhKRUWF0tPTm3xej8cjj8cTrbIBtDa3W7r7bqerAGCpqAWVoUOH6sMPP4xYd8MNN6hfv3664447lJWVpcTERK1evVqjR4+WJJWWlqqsrEx5eXnRKgsAAMSQqAWVTp066YwzzohYd8opp6hLly7h9ePHj9fkyZPVuXNneb1e3XrrrcrLy9O5554brbIAAEAMcfTOtL/+9a8VFxen0aNHq7KyUgUFBXr88cedLAlAa6upkV54IdS+4gpuoQ8ggssYY5wu4vsIBoPy+XwKBALyer1OlwOgpaqqpDlzQu3p00P
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"mask = (pos[:, 0] > 83.5) & (catx[\"totmass\"] > 1e14) & (pos[:, 0] < 84) & (catx[\"totmass\"] < 4e14)\n",
"\n",
"plt.figure()\n",
"plt.scatter(ra[mask], dec[mask], s=20, c=np.log10(catx[\"totmass\"][mask]))\n",
"plt.axvline(csiborgtools.hms_to_degrees(11, 44), zorder=0, ls=\"--\", color=\"red\", alpha=0.5)\n",
"plt.axhline(csiborgtools.dms_to_degrees(19, 45), zorder=0, ls=\"--\", color=\"red\", alpha=0.5)\n",
"# plt.axvline(285, zorder=0)\n",
"# plt.axhline(74, zorder=0)\n",
"plt.xlim(0, 360,)\n",
"plt.ylim(-90, 90)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([1.1764427e+14], dtype=float32),\n",
" array([83.93529855]),\n",
" array([184.62410722]),\n",
" array([315], dtype=int32))"
]
},
"execution_count": 62,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cat0[\"totmass\"][mask], cat0[\"dist\"][mask], pos[:, 1][mask], cat0[\"index\"][mask]"
]
},
{
"cell_type": "code",
"execution_count": 120,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"24"
]
},
"execution_count": 120,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mask = (cat0[\"dist\"] < 100) & (cat0[\"totmass\"] > 2e14)\n",
"mask.sum()"
]
},
{
"cell_type": "code",
"execution_count": 121,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([1.3149970e+15, 4.9394832e+14, 3.9707187e+14, 3.6969589e+14,\n",
" 3.5838989e+14, 3.4054645e+14, 3.2880589e+14, 3.2858516e+14,\n",
" 3.1561675e+14, 3.0209864e+14, 2.9957615e+14, 2.9370144e+14,\n",
" 2.8257113e+14, 2.7131389e+14, 2.5883651e+14, 2.5817466e+14,\n",
" 2.3839617e+14, 2.3532438e+14, 2.2714606e+14, 2.2489048e+14,\n",
" 2.1780584e+14, 2.1020473e+14, 2.1004685e+14, 2.0899980e+14],\n",
" dtype=float32),\n",
" array([ 2, 23, 31, 38, 39, 44, 49, 50, 56, 62, 65, 66, 69,\n",
" 75, 82, 83, 97, 100, 108, 109, 114, 122, 124, 125], dtype=int32))"
]
},
"execution_count": 121,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cat0[\"totmass\"][mask], cat0[\"index\"][mask]"
]
},
{
"cell_type": "code",
"execution_count": 133,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.25819606, 0.43737528, 0.51124454, 0.46384627, 0.26452401,\n",
" 0.32776746, 0.3836844 , 0.56561893, 0.22747411, 0.57970226,\n",
" 0.19699863, 0.55168611, 0.38103878, 0.19168761, 0.34359506,\n",
" 0.39914343, 0.74848348, 0.65658271, 0.61172086])"
]
},
"execution_count": 133,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"overlaps.max_overlap(0, True)[108]"
]
},
{
"cell_type": "code",
"execution_count": 135,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Cross main progenitors: 0%| | 0/19 [00:00<?, ?it/s]"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"Cross main progenitors: 100%|██████████| 19/19 [00:06<00:00, 2.98it/s]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending main progenitor for 17417.\n"
]
}
],
"source": [
"data = extract_main_progenitor_maxoverlap(108, overlaps, merger_trees)"
]
},
{
"cell_type": "code",
"execution_count": 136,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"mean = 14.2311, std = 0.1513\n",
"Found 2629 haloes in CB2 random.\n",
"Found 18572 haloes in MDPL2.\n"
]
}
],
"source": [
"logmp = np.log10(np.array([data[nsim][\"MainProgenitorMass\"][0] for nsim in data.keys()]))\n",
"\n",
"mean, std = np.mean(logmp), np.std(logmp)\n",
"print(f\"mean = {mean:.4f}, std = {std:.4f}\")\n",
"\n",
"bounds = [mean - 1 * std, mean + 1 * std]\n",
"\n",
"xrange_cb2, random_mah_cb2 = extract_mah(\"csiborg2_random\", bounds, \"MainProgenitorMass\", min_age=2.8)\n",
"xrange_mdpl2, random_mah_mdpl2 = extract_mah_mdpl2(bounds, min_age=2.8)\n",
"\n",
"print(f\"Found {len(random_mah_cb2)} haloes in CB2 random.\")\n",
"print(f\"Found {len(random_mah_mdpl2)} haloes in MDPL2.\")\n",
"\n",
"age, mah = summarize_extracted_mah(simname, data, nsim0, nsimxs, \"MainProgenitorMass\", min_age=2.8)\n",
"\n",
"\n",
"random_mah_cb2 /= random_mah_cb2[:, -1].reshape(-1, 1)\n",
"random_mah_mdpl2 /= random_mah_mdpl2[:, -1].reshape(-1, 1)\n",
"mah /= mah[:, -1].reshape(-1, 1)"
]
},
{
"cell_type": "code",
"execution_count": 137,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA+GklEQVR4nO3de3wU5d3///cmgSySZDVBTALhYCyhISQSaihYtZggyU1jUB8euCUiaFtT0VIULd5fjWltEW0RWynaVsWWHtT7LiBqoUjloAUDCalEFCGNyGEhQkpO3gm6O78/+GVvlmSTbNjj5PV8PPahM3vNtZ8dJuHNzDXXWAzDMAQAAICwFxHsAgAAAOAbBDsAAACTINgBAACYBMEOAADAJAh2AAAAJkGwAwAAMAmCHQAAgEkQ7AAAAEwiKtgFhAOn06kjR44oNjZWFosl2OUAAIA+xDAMNTU1KTk5WRERXZ+TI9j1wJEjR5SSkhLsMgAAQB928OBBDR06tMs2BLseiI2NlXR6h8bFxQW5GgAA0Jc0NjYqJSXFlUe6QrDrgfbLr3FxcQQ7AAAQFD0ZDsbNEwAAACZBsAMAADAJgh0AAIBJEOwAAABMgmAHAABgEgQ7AAAAkyDYAQAAmATBDgAAwCQIdgAAACZBsAMAADAJHikGAPArh9NQeW296ppaNTjWqpyR8YqM6P7RSAC8R7ADAPjNumq7ytbukb2h1bUuyWZVaWG68jOSglgZYE5cigUA+MW6artKVla6hTpJOtrQqpKVlVpXbQ9SZYB5EewAAD7ncBoqW7tHRifvta8rW7tHDmdnLQD0FsEOAOBz5bX1Hc7UncmQZG9oVXltfeCKAvoAgh0AwOfqmjyHut60A9AzBDsAgM8NjrX6tB2AnuGuWACAz+WMjFeSzaqjDa2djrOzSEq0nZ76BH0HU9/4H8EOAOBzkREWlRamq2RlpSySW7hr/2u8tDCdv9T7EKa+CQwuxQIA/CI/I0nLZ2Yr0eZ+uTXRZtXymdn8Zd6HMPVN4HDGDgDgN/kZSZqSnsjltz6su6lvLDo99c2U9ESOCx8g2AEA/CoywqKJqQnBLgNB4s3UNxwn545LsQAAwG+Y+iawCHYAAMBvmPomsAh2AADAb9qnvvE0es6i03fHMvWNbxDsAACA37RPfSOpQ7hj6hvfC6lgt2XLFhUWFio5OVkWi0WrV6/22Pauu+6SxWLR0qVLe9z/448/LovFonnz5p1zrQAAoGeY+iZwQuqu2JaWFmVlZWnOnDm6/vrrPbZbtWqVtm/fruTk5B73vWPHDj333HPKzMz0RakAAMALTH0TGCEV7AoKClRQUNBlm8OHD+uee+7R+vXrNW3atB7129zcrFtvvVW/+c1v9Nhjj/miVAAA4CWmvvG/kLoU2x2n06ni4mItWLBAY8aM6fF2d999t6ZNm6a8vLwetW9ra1NjY6PbCwAAINSF1Bm77ixevFhRUVG69957e7zNn//8Z1VWVmrHjh093mbRokUqKyvrTYkAAABBEzbBrqKiQk8//bQqKytlsfTsevzBgwf1/e9/Xxs2bJDV2vP5cRYuXKj58+e7lhsbG5WSkuJ1zQAAnAuH02BMGrwSNsFu69atqqur07Bhw1zrHA6H7rvvPi1dulSffPJJh20qKipUV1en7Oxst222bNmiZ555Rm1tbYqMjOywXXR0tKKjo/3yPQAA6Il11XaVrd3j9jiuJJtVpYXp3EUKj8Im2BUXF3cYIzd16lQVFxdr9uzZnW6Tm5ur3bt3u62bPXu2Ro8erQcffLDTUAcAQLCtq7arZGWljLPWH21oVcnKSqYIgUchFeyam5u1f/9+13Jtba2qqqoUHx+vYcOGKSHB/U6afv36KTExUWlpaa51ubm5uu666zR37lzFxsYqIyPDbZuBAwcqISGhw3oAAEKBw2mobO2eDqFOkgydntS3bO0eTUlP5LIsOgipu2J37typcePGady4cZKk+fPna9y4cXrkkUd63EdNTY2OHz/urxIBAPCr8tp6t8uvZzMk2RtaVV5bH7iiEDZC6ozdN7/5TRlGZ/9G6Vxn4+o6W3emTZs2eVcUAJgUA/NDU12T51DXm3boW0Iq2AEAAoOB+aFrcGzPZnHoaTv0LSF1KRYA4H/tA/PPvtzXPjB/XbU9SJVBknJGxivJZpWnc6cWnQ7hOSPjA1kWwgTBDgD6kO4G5kunB+Y7nD0fFgPfioywqLQwXZI6hLv25dLCdC6bo1MEOwDoQxiYHx7yM5K0fGa2Em3ul1sTbVamOkGXGGMHAH0IA/PDR35GkqakJ3KDC7xCsAOAPoSB+eElMsKiiakJ3TcE/n9cigWAPoSB+YC5EewAoA9hYD5gbgQ7AOhjGJgPmBdj7ACgD2JgPmBOBDsA6KMYmA+YD5diAQAATIIzdgAAwFQcTqPPDjMg2AEAANNYV21X2do9bk9YSbJZVVqY3iduDOJSLAAAMIV11XaVrKzs8Ni8ow2tKllZqXXV9iBVFjgEOwAAEPYcTkNla/fI6OS99nVla/fI4eyshXkQ7AAAQNgrr63vcKbuTIYke0OrymvrA1dUEBDsAABA2Ktr8hzqetMuXBHsAABA2Bsca+2+kRftwhXBDgAAhL2ckfFKslk7PAO5nUWn747NGRkfyLICjmAHADAFh9PQtpoTWlN1WNtqTph+kDzcRUZYVFqYLkkdwl37cmlhuunns2MeOwBA2Ovrc5fhtPyMJC2fmd3hWEjsQ8eCxTAM/knTjcbGRtlsNjU0NCguLi7Y5QAAztA+d9nZf5m1n5dZPjO7T/yFjv9jtidPeJNDOGMHAAhb3c1dZtHpucumpCeG9V/s8E5khEUTUxOCXUZQMMYOABC2mLsMcEewAwCELeYuA9wR7AAAYYu5ywB3BDsAQNhi7jLAHcEOABC2mLsMcEewAwCEtfa5yxJt7pdbE21WpjpBn8N0JwCAsJefkaQp6YmmmrsM6A2CHQDAFPry3GVAOy7FAgAAmATBDgAAwCS4FAvAb8z2vEYACHUEOwB+sa7arrK1e9we95Rks6q0MJ27FAHAT7gUC8Dn1lXbVbKyssMzPI82tKpkZaXWVduDVBkAmBvBDoBPOZyGytbukdHJe+3rytbukcPZWQsAwLkg2AHwqfLa+g5n6s5kSLI3tKq8tj5wRQFAH0GwA+BTdU2eQ11v2gEAeo5gB8CnBsdau2/kRTsAQM8R7AD4VM7IeCXZrB0eyN7OotN3x+aMjA9kWQDQJxDsAPhUZIRFpYXpktQh3LUvlxamM58dAPgBwQ6Az+VnJGn5zGwl2twvtybarFo+M5t57ADAT5igGIBf5GckaUp6Ik+eAIAAItgB8JvICIsmpiYEuwwA8ItQfGwiwQ4AAMBLofrYRMbYAQAAeCGUH5tIsAMAAOihUH9sIsEOAACgh0L9sYkEOwAAgB4K9ccmEuwAAAB6KNQfm0iwAwAA6KFQf2wiwQ4AAKCHQv2xiQQ7AAAAL4TyYxOZoBgAAMBLofrYRIIdAABAL4TiYxND6lLsli1bVFhYqOTkZFksFq1evdpj27vuuksWi0VLly7tss/ly5crMzNTcXFxiouL08SJE/XXv/7Vt4UDAACEgJAKdi0tLcrKytKyZcu6bLdq1Spt375dycnJ3fY5dOhQPf7446qoqNDOnTt19dVXq6ioSB988IGvygYAAAgJIXUptqCgQAUFBV22OXz4sO655x6tX79e06ZN67bPwsJCt+Wf/OQnWr58ubZv364xY8acU70AAAChJKSCXXecTqeKi4u1YMGCXoUyh8OhV199VS0tLZo4caIfKgQAAAiesAp2ixcvVlRUlO69916vttu9e7cmTpyo1tZWxcTEaNWqVUpPT/fYvq2tTW1tba7lxsbGXtcMAAAQKCE1xq4rFRUVevrpp7VixQpZLN7dSpyWlqaqqiq99957Kikp0axZs7Rnzx6P7RctWiSbzeZ6paSknGv5AAAAfmcxDMMIdhGdsVgsWrVqlaZPny5JWrp0qebPn6+IiP/Log6HQxEREUpJSdEnn3zS477z8vKUmpqq5557rtP3Ozt
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure()\n",
"# plt.hist(logmp, bins=\"auto\")\n",
"plt.scatter(np.arange(len(logmp)), logmp)\n",
"plt.xlabel(\"Simulation index\")\n",
"plt.ylabel(r\"$\\log M ~ [M_\\odot / h]$\")\n",
"# plt.axhline(mean, c=\"red\", ls=\"--\")\n",
"plt.tight_layout()\n",
"plt.savefig(\"../../plots/COMA_MASS.png\", dpi=450)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 138,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAGGCAYAAACNCg6xAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAADimklEQVR4nOz9aaxkeX7XDX7OvsS+3X3Nfa2lq7p6p922x6ZhbDEjBsQbDEh4kMwbLITAQvgF4h0YBCpkDY/AIzHPwACP/WAb48Ztt912t129Vld3V2VWbjfvHvt29m1enHuzsrq2XO4S9+b5SEcRcW/kjX9GnDjne37L9yckSZKQkZGRkZGRkZHxkYjHvYCMjIyMjIyMjJNCJpwyMjIyMjIyMh6RTDhlZGRkZGRkZDwimXDKyMjIyMjIyHhEMuGUkZGRkZGRkfGIZMIpIyMjIyMjI+MRyYRTRkZGRkZGRsYjkgmnjIyMjIyMjIxHJBNOGRkZGRkZGRmPSCacMjIyMjIyMjIekUw4ZWRkZGRkZGQ8Is+ccPqt3/otLl68yPnz5/nf/rf/7biXk5GRkZGRkXGCEJ6lIb9hGHLlyhX+4A/+gFKpxEsvvcTXvvY1arXacS8tIyMjIyMj4wTwTEWcXnvtNa5evcr8/Dz5fJ4vfvGLfOlLXzruZWVkZGRkZGScEE6UcPqjP/ojfuZnfoa5uTkEQeA3fuM33vOcV199lZWVFXRd5xOf+ASvvfbag99tbW0xPz//4PH8/Dybm5tHsfSMjIyMjIyMU8CJEk6WZfH888/z6quvvu/v//N//s/84i/+Ir/8y7/Mt7/9bZ5//nl++qd/mmazecQrzcjIyMjIyDiNyMe9gMfhi1/8Il/84hc/8Pe/8iu/wt/+23+bv/k3/yYAv/qrv8pv//Zv8+///b/nH/7Df8jc3Ny7Ikybm5u88sorH/j3PM/D87wHj+M4ptvtUqvVEAThAP5HGRkZGRkZGcdNkiSMRiPm5uYQxY+IKSUnFCD59V//9QePPc9LJEl618+SJEn++l//68nP/uzPJkmSJEEQJOfOnUs2NjaS0WiUXLhwIWm32x/4Gr/8y7+cANmWbdmWbdmWbdn2DGzr6+sfqT9OVMTpw2i320RRxPT09Lt+Pj09zVtvvQWALMv8i3/xL/jCF75AHMf8g3/wDz60o+4f/aN/xC/+4i8+eDwYDFhaWmJ9fZ1isXg4/5GMjGecMATPSzfff2dz3QTLibHdBM9NCJKEKE6IooQwSYjjBEGMQYwRpATEBFFMkCUQJZAkEKUkvRVBkvd/BpKY3oonKJCcJOAH4DnguQLO3u3+5jrgumDZYO9tnpcQRRAlCXECcZwQJQlJDHGSkAgJohwjSDFI6XuJlCAIMaKUpO+RlLyzie9+LIggiO/cigIIwkM/F0AU3/1YEAAhgYfe+zSgn+z9bu9X+889BpIYQl/Ed0R8TyJ0RXxXJNi/9UQCVyLwREJPwPdFAk9If+4LhJ5AGIhAQpLsnaETSEjeeY2EdP+VEoS9z+DhjYfui0r6mQhSsvd+P3w/3RBjEJMf+fnea+x9Lux9NgjJ3uP0c9l/vP87KYlQ4gA5CpHjACUKkZMQJQlRogA18tEjHy32UOIAUUiQiJGECDGJUKMALfaR4wg5CZHjECUOkeIYKU4fq1GAnERIRIhCgkC6ZpF0LbEskkgCyAJIkAiAKJCIIokoEIuk9wWBWBKJBZEgVAh9lcBXCX2N2JMp2EPyYxsfnUFxluHSLEoR1n4o8G+//A8oFAofuT+cGuH0qPzsz/4sP/uzP/tIz9U0DU3T3vPzYrGYCaeMjCckDNMTuuOkm+uCZSWMnZixneD5MUGcEEYxQbR3UpDTE4qigFIQ0CuQlwVkBWQ5QZZBVo7vxHqQhAHYtoBjCTi2gG3z4L41FhiNE4YjGI8hCBPCeE9AxglRvHdyVaIHm6zGyHqMUorRlRhJSZDlBGnvvvTQ/Y/KULwb4UduJ58kgcAV8WwJz5ZSIbR3mz5Of+7aIp4t4jnvCKMkSd4lePZvRSVC1CJEJQQlQpBDRCVGzIdIlQhJjlGUCFGN0+fKUfp7Ndp7HD8Qq8IRVR2LcYQSBqmQifaETBSm96MIJUyFkB55aJGHQpgKIaK92xAkAURQxAA5DpGEMBU+SYQcRWiRjxp5aFGAKMSIQoK89++REgQVUAQSWSCSJSJZJpI1IlnceywRSRKxJJJIH/3GxJ5M4CiEtkLsqESugpAIiALkZZsZuY0Z2gSNAoOXL2EvzoJZwl/X4bZFKZEAHqkM59QIp3q9jiRJ7O7uvuvnu7u7zMzMHNOqMjKeTYLgnSjHvjgajmJGVozlxgRRjB8mhMneSUSOUdUERQOlmJBTQVFTofR4J/PJJY7BscG2UgFkWwL2WMAai4xGCYMBDIfguKkYCuOEMIqIkr2TrBoiqBGKHqEYMUolJqfFKHubrMXIanxq3q9HIYnBcyQ8K93cvVvPlnDH6X3HEnEsCc9OI0Zxsid+Eh6IIUGOELUwfY81F1ENkcwIsRyhaiG6FiGpIZK6J5L276vRRIp1KQrRQh8lDNLbKEALfIzQJh84aLGH+i4xFJGIgARIAokoEKkikSQTSSKRICIHEbKf/h3N9ynYYwq+hRr7KISoQgiKADIkukCkpOInVCQCRX1HDCkSyVPupHEoEjoqkb0vlBSIRURBQDNjzFKENu9SjodUxh0YwShfY3fuAr3ZWQJJpbulM3wzxxQhL8712X5uBX7/0V7/1AgnVVV56aWX+PKXv8xf+kt/CUiLub/85S/zd//u3z3exWVknELiOBVFDwuk4SimN4yxvRg/iPGiGEGKEOT0pK5qoNUT8ipoeoKsHPf/4uDwPRiPhHQbioyGAtZIoD/cF0VJKhajhCCOieJ47yQcIGkhihGjNmIMI0LRY1Q9RtEjZDWZyJPzYZEk4Dsi7ljGHUs4o/T2weOxhD0SU5Fki8TxfuorTT8ixkh6gKinIkjSAqRqiDQbktdDJG1/i5D0MBVMUvJRyzo2hDhOIzpxtJfaipDiCDGJkffuK2GAHriYoYsZOahJgEyIQoQkRCQSCIqQihhDJpRkHNkgkkRCSSYWRRAEhDhGDXxM18F0HErjIWV7SNEboSU+GgGCBmgQFGQCVSHQNBwtT6Aqh5LrThIhFUgPRZNiX0IUBRQloVCMMKYCjEKEnouQpAijM8TYHRG4Kp3CPO0rywwbDWJE+jsavY0ceVnjlbMeV/VNwqWr3PPOPPKaTpRwGo/H3Lp168Hju3fv8t3vfpdqtcrS0hK/+Iu/yM/93M/x8ssv88orr/Cv/tW/wrKsB112GRkZj08cg2Wl4siyYDxO6A1i+uMIP0rwg5iINO0gqTG6AVohoWQkaPrpiRjtC6Nh/x1RNBpCtwf9AYythGAvvRgmfnpC1kJkPUQzItSliJIZo+oRqplGiZ41QRR4Is5QxhlK2CMZdyjjjmWckYQ9krCGaaQojNLnx0maEhOUCMnwEfQASbeQiiHyVEDeCJH0AFkPkfbui8rkv69iHD2UGkvTZEoUPPiZFnnooZdGhuLgoRRZgkiEQFpDlIgCggiJDIkoEJipKApkGVfWCSWZSJLek8MWo4icY1OwXQzXRXcdit6YgjdGS9IIkixFYAgEBRm/oWEbefq6+khps6chDiRCSyWwVCJHIXYVRARECcx8hDEToed9jEKIosUPMsWiH2I2B6hNCzfOsVG7QOfiAnalQpLAsKXSvZ9DR+P58yKXzjlUm+vY1fPYl67CG96HL+whTpRw+uY3v8kXvvCFB4/3C7d/7ud+jl/7tV/jr/7Vv0qr1eKf/JN/ws7ODi+88AL/83/+z/cUjGdkZLyXJEnF0XicCiTLSugOYvqjCD+M8YKYRErrZjQ9QS8nlA3QT0nkKEnAHgsMhwLDvsBoKDAaiPR6CZ0uDEfJXooxImIvqqEFKGaEVokwFyI0M91U/ejqVSaBJAHXknCHMvZQxhml4sgZylhDmfFAwh5KhL6QFqXvRYckLUA0AkTDRTJ85NkQw/SRjRDZ9JHMANkIEOXJjwgp0Tt1QlI
"text/plain": [
"<Figure size 600x400 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"pmin, pmax = 0.15, 100 - 0.15\n",
"# norm_mass = 10**np.nanmean(np.log10(mah[:, -1]))\n",
"norm_mass = 1\n",
"alphas = {3: 0.25, 2: 0.5, 1: 0.75}\n",
"\n",
"\n",
"plt.figure(figsize=(6, 4))\n",
"cols = plt.rcParams[\"axes.prop_cycle\"].by_key()[\"color\"]\n",
"lw = plt.rcParams[\"lines.linewidth\"]\n",
"\n",
"for n in [3, 2, 1]:\n",
" pmin, pmax = norm.cdf(-n) * 100 , norm.cdf(n) * 100\n",
" ylow, yhigh = np.percentile(random_mah_mdpl2 / norm_mass, [pmin, pmax], axis=0)\n",
"\n",
" ylow = savgol_filter(ylow, 5, 1, mode=\"interp\")\n",
" yhigh = savgol_filter(yhigh, 5, 1, mode=\"interp\")\n",
" plt.fill_between(xrange_mdpl2, ylow, yhigh, alpha=alphas[n], color=cols[0],\n",
" facecolor=cols[0], zorder=0, edgecolor=\"blue\",\n",
" label=f\"MDPL2 Random ({len(random_mah_mdpl2)})\" if n == 1 else None)\n",
" \n",
" # pmin, pmax = norm.cdf(-n) * 100 , norm.cdf(n) * 100\n",
" # ylow, yhigh = np.percentile(random_mah_cb2 / norm_mass, [pmin, pmax], axis=0)\n",
" # ylow = savgol_filter(ylow, 5, 1, mode=\"interp\")\n",
" # yhigh = savgol_filter(yhigh, 5, 1, mode=\"interp\")\n",
" # plt.fill_between(xrange_cb2, ylow, yhigh, alpha=alphas[n], color=cols[0],\n",
" # facecolor=cols[0], zorder=0, edgecolor=\"blue\",\n",
" # label=f\"CB2 Random ({len(random_mah_cb2)})\" if n == 1 else None)\n",
"\n",
" if n == 3:\n",
" continue\n",
"\n",
" ylow, yhigh = np.percentile(mah / norm_mass, [pmin, pmax], axis=0)\n",
" ylow = savgol_filter(ylow, 5, 1, mode=\"interp\")\n",
" yhigh = savgol_filter(yhigh, 5, 1, mode=\"interp\")\n",
" plt.fill_between(age, ylow, yhigh, alpha=alphas[n] / 2,\n",
" label=\"Coma\" if n == 1 else None, zorder=1,\n",
" color=cols[1], edgecolor=\"red\")\n",
"\n",
"# pmin, pmax = norm.cdf(-1) * 100 , norm.cdf(1) * 100\n",
"# ylow, yhigh = np.percentile(random_mah_cb2 / norm_mass, [pmin, pmax], axis=0)\n",
"# plt.plot(xrange_cb2, ylow, c=cols[2], ls=\"--\", lw=1.5 * lw,\n",
"# label=rf\"$1\\sigma$ CB2 Random ({len(random_mah_cb2)})\")\n",
"# plt.plot(xrange_cb2, yhigh, c=cols[2], ls=\"--\", lw=1.5 * lw)\n",
"\n",
"# for i in range(len(mah)):\n",
" # plt.plot(age, mah[i] / norm_mass, c=\"black\", lw=lw/4, alpha=0.5, zorder=4)\n",
"\n",
"# plt.plot(RANDOM_MAH_Sorce_Virgo_UPPER[:, 0],\n",
" # RANDOM_MAH_Sorce_Virgo_UPPER[:, 1], c=\"red\", ls=\"--\", label=\"Sorce+2016 Virgo randoms\")\n",
"# plt.plot(RANDOM_MAH_SORCE_Virgo_LOWER[:, 0],\n",
" # RANDOM_MAH_SORCE_Virgo_LOWER[:, 1], c=\"red\", ls=\"--\")\n",
"\n",
"\n",
"\n",
"plt.yscale(\"log\")\n",
"plt.legend(loc=\"lower right\")\n",
"\n",
"plt.xlabel(r\"$\\mathrm{Age} ~ [\\mathrm{Gyr}]$\")\n",
"plt.ylabel(r\"$M_{\\rm main}(t) / M_{\\rm main}(z = 0)$\")\n",
"plt.xlim(xrange_mdpl2.min(), xrange_mdpl2.max())\n",
"plt.ylim(0.01, 1)\n",
"plt.tight_layout()\n",
"# plt.savefig(\"../../plots/coma_MAIN.png\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 79,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/tmp/ipykernel_3803914/1215582327.py:31: RuntimeWarning: All-NaN slice encountered\n",
" mu = np.nanmedian(random_mah, axis=0)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnAAAAHWCAYAAAD3vrTNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd5hc5Xn4/e8p0/vM7s72rrq76hKIJqopNjY2tjF2CI5xS4xL+NlxL7ET4/LGIQ4EEvcaG/cYMBhEVUW9rrb3Nr3XU94/FhaEChIsRXA+1zXXrGbOeU6Z0Tn3PO0WdF3XMRgMBoPBYDCcMcRXegcMBoPBYDAYDKfHCOAMBoPBYDAYzjBGAGcwGAwGg8FwhjECOIPBYDAYDIYzjBHAGQwGg8FgMJxhjADOYDAYDAaD4QxjBHAGg8FgMBgMZxgjgDMYDAaDwWA4w8iv9A6cCTRNY3JyEpfLhSAIr/TuGAwGg8FwynRdJ51OU1tbiyi+PPU2hUKBUqk0b+WZzWasVuu8lfdaYARwp2BycpKGhoZXejcMBoPBYHjBxsbGqK+vf8m3UygUaGlyMh1S563M6upqhoaGjCDuWYwA7hS4XC5g9svvdrtf4b0xGAwGg+HUpVIpGhoa5u5lL7VSqcR0SGVkVzNu14uv8UulNZpWD1MqlYwA7lmMAO4UPN1s6na7jQDOYDAYDGekl7sLkNMl4HS9+G1qGF2XjscI4AwGg8FgMMw7VddQ9fkpx3AsYxSqwWAwGAwGwxnGqIEzGAwGg8Ew7zR0NF58Fdx8lPFaZARw80RVVcrl8iu9G4bXILPZ/LIN/TcYDIb5oqExH42f81PKa48RwL1Iuq4zPT1NIpF4pXfF8BoliiItLS2YzeZXelcMBoPB8CphBHAv0tPBW1VVFXa73Zjo1zCvnp5EempqisbGRuP7ZTAYzhiqrqPqL775cz7KeC0yAriTuOOOO7jjjjtQ1eNPRqiq6lzwFggEXua9M7xeVFZWMjk5iaIomEymV3p3DAaD4ZQYfeBeWkbHmpP4yEc+wuHDh9mxY8dx33+6z5vdbn85d8vwOvN00+mJfkgYDAaD4fXHqIGbB0azluGlZHy/DAbDmUhDRzVq4F4yRg2cwWAwGAwGwxnGCOAML8rmzZvp6urCZDJxzTXXvNK7M+++8pWvsGLFihdVxvDwMIIgsHfv3nnZJ4PBYDgTPN0Hbj4ehmMZAdzr1Hvf+14EQUAQBEwmEy0tLfzTP/0ThULhtMq55ZZbWLFiBUNDQ/z4xz9+aXbWYDAYDGecp0ehzsfDcCyjD9zr2BVXXMGPfvQjyuUyu3bt4sYbb0QQBL75zW+echkDAwN8+MMfpr6+/gXvR6lUMuY4MxgMBoPhNBg1cK9jFouF6upqGhoauOaaa7j00kt58MEH597XNI1bb72VlpYWbDYby5cv57e//S3wTLNgNBrlfe97H4IgzNXAHTx4kCuvvBKn00kwGOSGG24gEonMlXvhhRdy880384lPfIKKigouv/zyU17vYx/7GP/0T/+E3++nurqar3zlK0cdUyKR4EMf+hDBYBCr1UpnZyf33HPP3PubNm3i/PPPx2az0dDQwMc+9jGy2ezznquf/exnNDc34/F4eNe73kU6nZ577/777+e8887D6/USCAR405vexMDAwEnLe+yxx1i3bh0Wi4Wamho+85nPoCjK8+6HwWAwnCm0eXwYjmUEcAZgNnjasmXLUTVht956Kz/96U+56667OHToEP/4j//I3/zN3/DYY4/R0NDA1NQUbreb2267jampKa677joSiQQXX3wxK1euZOfOndx///3MzMzwzne+86jt/eQnP8FsNrN582buuuuu01rP4XCwfft2vvWtb/HVr351LujUNI0rr7ySzZs38/Of/5zDhw/zjW98A0mSgNnawiuuuIJrr72W/fv38+tf/5pNmzZx8803n/TcDAwM8Mc//pF77rmHe+65h8cee4xvfOMbc+9ns1luueUWdu7cycaNGxFFkbe+9a1o2vEvOxMTE1x11VWsXbuWffv2ceedd/KDH/yAf/mXfzn1D8xgMBhe5dSnRqHOx8NwLKMJ9aVw001w6NDLu82ODvjBD05rlXvuuQen04miKBSLRURR5PbbbwegWCzy9a9/nYceeoj169cD0NrayqZNm/jv//5vNmzYQHV1NYIg4PF4qK6uBuDf/u3fWLlyJV//+tfntvPDH/6QhoYGent7WbhwIQALFizgW9/61twy//Iv/3JK6y1btowvf/nLc2XcfvvtbNy4kcsuu4yHHnqIJ598ku7u7rnlW1tb58q79dZbec973sMnPvGJufW/+93vsmHDBu68806sVutxz5Omafz4xz/G5XIBcMMNN7Bx40b+9V//FYBrr732qOV/+MMfUllZyeHDh+ns7DymvP/6r/+ioaGB22+/HUEQWLx4MZOTk3z605/mS1/6kpH31GB4ndB1nXw+TyqVIhaLEY1GiUQiTExM4PP5eM973vNK76LhVcwI4F7HLrroIu68806y2Sz//u//jizLc8FIf38/uVyOyy677Kh1SqUSK1euPGGZ+/bt45FHHsHpdB7z3sDAwFxgtXr16he03rJly456r6amhlAoBMDevXupr6+fW/Z4+7Z//35+8YtfzL2m6zqapjE0NMSSJUuOu15zc/Nc8PbcbQL09fXxpS99ie3btxOJROZq3kZHR48bwHV3d7N+/fqj5nc799xzyWQyjI+P09jYeNz9MBgMZ7ZkMkl3dzdHjhxhamqKSCRCNpulWCySz+fnulEUCgU6OjrO+ABO1Wcf81GO4VhGAPdSOM2asFeKw+Ggvb0dmK01Wr58OT/4wQ+46aabyGQyANx7773U1dUdtZ7FYjlhmZlMhquvvvq4AyFqamqO2vYLWe+5qaQEQZgLmGw22wn36+ltfOhDH+JjH/vYMe+dLGg62TYBrr76apqamvje975HbW0tmqbR2dlJqVQ66f4YDIbXNl3XCYVC9Pf3z/2AjEaj6LqOzWbDYrFgsVjw+XzU1NTMdWEZGxt7hfd8fsxX/zWjD9zxGQGcAQBRFPnc5z7HLbfcwrvf/W6WLl2KxWJhdHSUDRs2nHI5q1at4ne/+x3Nzc3I8ql/vV7oes+2bNkyxsfHj2pyfe42Dh8+PBe0zodoNEpPTw/f+973OP/884HZgRIns2TJEn73u9+h6/pcLdzmzZtxuVwvajSvwWB4Zem6zuTkJENDQwwPD3Pw4EFCoRCpVApRFAkEAixatOgFX+MMhmczOtsY5rzjHe9AkiTuuOMOXC4Xn/zkJ/nHf/xHfvKTnzAwMMDu3bv5z//8T37yk5+csIyPfOQjxGIxrr/+enbs2MHAwAAPPPAAf/d3f3fSXJ4vdL1n27BhAxdccAHXXnstDz74IENDQ/zlL3/h/vvvB+DTn/40W7Zs4eabb2bv3r309fXxpz/96XkHMZyMz+cjEAjwP//zP/T39/Pwww9zyy23nHSdf/iHf2BsbIyPfvSjHDlyhD/96U98+ctf5pZbbjH6vxkMZxhVVRkYGODee+/lq1/9Kl/4whe47bbb+OMf/8j09DQul4slS5awdOlSgsHg6yp40xBQ5+GhYaQTPJ7XzzfJ8LxkWebmm2/mW9/6Fn//93/P1772NSorK7n11lsZHBzE6/WyatUqPve5z52wjNraWjZv3synP/1p3vCGN1AsFmlqauKKK644aXDyQtd7rt/97nd88pOf5PrrryebzdLe3j43YnTZsmU89thjfP7zn+f8889H13Xa2tq47rrrTv0kPYcoivzqV7/iYx/7GJ2dnSxatIjvfve7XHjhhSdcp66ujvvuu49PfepTLF++HL/fz0033cQXvvCFF7wfBoPh5aPrOiMjI+zfv5+tW7cyPj5ONpvFbrdTUVFBY2PjSXMY64AiQFkQKIvPPKsIaAKE/U4Uu+mE6xsMAIKuG1McP59UKoXH4yGZTOJ2u+deLxQKDA0N0dLScsIRjAbDi2V8zwyGV97T/dkOHjzI1m3b6B7oJ1EsYPG6sQX8iA47RUmgKAqzAZkAijj7XBI
"text/plain": [
"<Figure size 640x480 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"\n",
"cmap = plt.cm.viridis\n",
"\n",
"x = np.asarray([data[nsimx][\"Overlap\"] for nsimx in nsimxs if nsimx in data])\n",
"w = x - x.min()\n",
"w /= w.max()\n",
"\n",
"norm = Normalize(vmin=x.min(), vmax=x.max())\n",
"sm = ScalarMappable(cmap=cmap, norm=norm)\n",
"sm.set_array(x)\n",
"\n",
"\n",
"fig, ax = plt.subplots()\n",
"lw = plt.rcParams[\"lines.linewidth\"] / 2\n",
"d = data[nsim0]\n",
"ax.plot(d[\"Age\"], d[\"MainProgenitorMass\"], color=\"red\", lw=1.5*lw,\n",
" label=\"Reference halo\", zorder=1)\n",
"\n",
"i = 0\n",
"for nsimx in nsimxs:\n",
" try:\n",
" d = data[nsimx]\n",
" except KeyError:\n",
" continue\n",
"\n",
" ax.plot(d[\"Age\"], d[\"MainProgenitorMass\"], color=cmap(norm(x[i])),\n",
" zorder=0, lw=lw * (x[i] / x.max()))\n",
"\n",
" i += 1\n",
"\n",
"\n",
"mu = np.nanmedian(random_mah, axis=0)\n",
"ymin = np.nanpercentile(random_mah, 16, axis=0)\n",
"ymax = np.nanpercentile(random_mah, 84, axis=0)\n",
"ax.plot(xrange, mu, color=\"black\", label=\"Random\")\n",
"ax.fill_between(xrange, ymin, ymax, color=\"black\", alpha=0.5, zorder=-1)\n",
"\n",
"cbar = fig.colorbar(sm, ax=ax, orientation='vertical')\n",
"cbar.set_label('Overlap')\n",
"\n",
"ax.legend()\n",
"ax.set_yscale(\"log\")\n",
"ax.set_xlabel(r\"$t ~ [\\mathrm{Gyr}]$\")\n",
"ax.set_ylabel(r\"$M_{\\rm main} ~ [M_{\\odot}]$\")\n",
"ax.set_xlim(0)\n",
"plt.tight_layout()\n",
"# plt.savefig(\"../../plots/example_mah.png\")\n",
"fig.show()"
]
},
{
"cell_type": "code",
"execution_count": 136,
"metadata": {},
"outputs": [],
"source": [
"nsim0 = 1\n",
"simname = \"csiborg2_random\"\n",
"kind = simname.split(\"_\")[-1]\n",
"min_logmass = 12.25\n",
"\n",
"# NOTE: These can possibly be pickled to avoid doing this long process every\n",
"# single time.\n",
"\n",
"cat = csiborgtools.read.CSiBORG2Catalogue(\n",
" nsim0, 99, kind, bounds={\"totmass\": (1e13, None), \"dist\": (0, 135)})\n",
"# merger_reader = csiborgtools.read.CSiBORG2MergerTreeReader(nsim0, kind)"
]
},
{
"cell_type": "code",
"execution_count": 149,
"metadata": {},
"outputs": [],
"source": [
"a = np.array([0.01428571, 0.02424242, 0.03419913, 0.04415584, 0.05411255,\n",
" 0.06406926, 0.07402597, 0.08398268, 0.09393939, 0.1038961 ,\n",
" 0.11385281, 0.12380952, 0.13376623, 0.14372294, 0.15367965,\n",
" 0.16363636, 0.17359307, 0.18354978, 0.19350649, 0.2034632 ,\n",
" 0.21341991, 0.22337662, 0.23333333, 0.24329004, 0.25324675,\n",
" 0.26320346, 0.27316017, 0.28311688, 0.29307359, 0.3030303 ,\n",
" 0.31298701, 0.32294372, 0.33290043, 0.34285714, 0.35281385,\n",
" 0.36277056, 0.37272727, 0.38268398, 0.39264069, 0.4025974 ,\n",
" 0.41255411, 0.42251082, 0.43246753, 0.44242424, 0.45238095,\n",
" 0.46233766, 0.47229437, 0.48225108, 0.49220779, 0.5021645 ,\n",
" 0.51212121, 0.52207792, 0.53203463, 0.54199134, 0.55194805,\n",
" 0.56190476, 0.57186147, 0.58181818, 0.59177489, 0.6017316 ,\n",
" 0.61168831, 0.62164502, 0.63160173, 0.64155844, 0.65151515,\n",
" 0.66147186, 0.67142857, 0.68138528, 0.69134199, 0.7012987 ,\n",
" 0.71125541, 0.72121212, 0.73116883, 0.74112554, 0.75108225,\n",
" 0.76103896, 0.77099567, 0.78095238, 0.79090909, 0.8008658 ,\n",
" 0.81082251, 0.82077922, 0.83073593, 0.84069264, 0.85064935,\n",
" 0.86060606, 0.87056277, 0.88051948, 0.89047619, 0.9004329 ,\n",
" 0.91038961, 0.92034632, 0.93030303, 0.94025974, 0.95021645,\n",
" 0.96017316, 0.97012987, 0.98008658, 0.99004329, 1. ])"
]
},
{
"cell_type": "code",
"execution_count": 155,
"metadata": {},
"outputs": [],
"source": [
"from h5py import File"
]
},
{
"cell_type": "code",
"execution_count": 156,
"metadata": {},
"outputs": [],
"source": [
"f = File(\"/mnt/extraspace/rstiskalek/csiborg_postprocessing/random_mah/random_mah_csiborg2_random_1.hdf5\")"
]
},
{
"cell_type": "code",
"execution_count": 168,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ nan, nan, nan, nan,\n",
" nan, nan, nan, nan,\n",
" nan, nan, nan, nan,\n",
" nan, 9.6507642e+10, 3.0197550e+11, 3.2376760e+11,\n",
" 7.0668498e+11, 7.7206113e+11, 9.3394495e+11, 1.5814801e+12,\n",
" 1.9799633e+12, 2.8672110e+12, 3.4369173e+12, 4.2152047e+12,\n",
" 4.0595470e+12, 6.7960056e+12, 8.0506048e+12, 9.8562317e+12,\n",
" 1.5811687e+13, 2.1929025e+13, 2.7361473e+13, 3.6311776e+13,\n",
" 4.0548774e+13, 4.6049712e+13, 5.1784131e+13, 5.8592594e+13,\n",
" 6.2322144e+13, 6.8669858e+13, 7.1316036e+13, 1.2149067e+14,\n",
" 1.4343837e+14, 1.5714247e+14, 1.6746878e+14, 1.8382215e+14,\n",
" 1.8707539e+14, 1.8938536e+14, 1.9296859e+14, 1.9942214e+14,\n",
" 2.0502892e+14, 2.0806736e+14, 2.1257209e+14, 2.3088363e+14,\n",
" 2.4707200e+14, 2.5878680e+14, 2.6422236e+14, 2.7754664e+14,\n",
" 2.8651874e+14, 2.9785370e+14, 3.0968989e+14, 3.1373079e+14,\n",
" 3.1944651e+14, 3.2246628e+14, 3.2166309e+14, 3.2575374e+14,\n",
" 3.2939614e+14, 3.3505896e+14, 3.4118875e+14, 3.4710996e+14,\n",
" 3.5778809e+14, 3.6322364e+14, 3.7472670e+14, 3.8465764e+14,\n",
" 4.0075266e+14, 4.1667327e+14, 4.3123039e+14, 4.3992541e+14,\n",
" 4.6098277e+14, 4.7190678e+14, 4.9967297e+14, 5.0868244e+14,\n",
" 5.2763220e+14, 5.5073174e+14, 5.9031232e+14, 5.9930933e+14,\n",
" 6.1641921e+14, 6.2424880e+14, 6.3106035e+14, 6.4602838e+14,\n",
" 6.5189973e+14, 6.5782102e+14, 6.6686783e+14, 6.7132272e+14,\n",
" 6.7323109e+14, 6.7524221e+14, 6.7571224e+14, 6.7606094e+14,\n",
" 6.8144662e+14, 6.8422983e+14, 6.8603237e+14, 6.8886215e+14],\n",
" dtype=float32)"
]
},
"execution_count": 168,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f[\"MainProgenitorMass\"][0]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 154,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[69.0000210000063, 40.250007218751264, 28.24050991940438, 21.6470609550175, 17.480001404480106, 14.608109099433955, 12.508772664512199, 10.90721705951751, 9.64516173673259, 8.625000360937513, 7.7832702592057235, 7.0769233254437935, 6.475728365821477, 5.95783150553419, 5.50704240932355, 5.111111246913583, 4.760598622974984, 4.448113312911626, 4.1677853285437605, 3.914893700679041, 3.685598452365574, 3.476744253718227, 3.285714346938776, 3.1103203402819117, 2.9487179993425383, 2.7993421515051513, 2.6608558268213116, 2.5321101306287352, 2.4121122957547967, 2.3000000330000008, 2.1950207773798662, 2.096514773533915, 2.003901196522936, 1.9166666909722223, 1.8343558508261513, 1.7565632668759008, 1.6829268488994646, 1.613122190273029, 1.5468577900064306, 1.4838709837669097, 1.4239244641145379, 1.366803292753544, 1.3123123255056859, 1.2602739849878026, 1.210526327423823, 1.162921359250726, 1.117323566656109, 1.0736086272735772, 1.0316622782422846, 0.9913793189283591, 0.9526627299814432, 0.9154228931957131, 0.8795768989699038, 0.8450479301016136, 0.8117647122768166, 0.7796610229819017, 0.7486752517178681, 0.7187500053710938, 0.6898317534223188, 0.6618705083794834, 0.6348195374209455, 0.6086351017498701, 0.5832762206018658, 0.5587044572276223, 0.5348837244997295, 0.5117801080759505, 0.48936170529651424, 0.46759847820604516, 0.4464621192761633, 0.42592592856652933, 0.4059647012034677, 0.3865546241790834, 0.3676731815824261, 0.34929906746973005, 0.3314121056648591, 0.31399317585528075, 0.2970241454144613, 0.28048780643961924, 0.2643678175452504, 0.2486486499985392, 0.23331553782343795, 0.21835443153641232, 0.20375195520916023, 0.18949536658248856, 0.17557251998135315, 0.1619718318042056, 0.14868224838055033, 0.13569321600925854, 0.122994653006949, 0.11057692361085425, 0.09843081359419292, 0.08654750746436402, 0.0749185671253807, 0.06353591189600438, 0.05239179978414388, 0.04147880992632613, 0.03078982610853953, 0.020318021291547472, 0.010056843069963017, 0.0]\n"
]
}
],
"source": [
"print(list(1 / a - 1))"
]
},
{
"cell_type": "code",
"execution_count": 142,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"17"
]
},
"execution_count": 142,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.sum(cat[\"totmass\"] > 5e14)"
]
},
{
"cell_type": "code",
"execution_count": 130,
"metadata": {},
"outputs": [],
"source": [
"d1 = merger_reader.main_progenitor(10000)\n"
]
},
{
"cell_type": "code",
"execution_count": 132,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([99, 98, 97, 96, 95, 94, 93, 92, 91, 90, 89, 88, 87, 86, 85, 84, 83,\n",
" 82, 81, 80, 79, 78, 77, 76, 75, 74, 73, 72, 71, 70, 69, 68, 67, 66,\n",
" 65, 64, 63, 62, 61, 60, 59, 58, 57, 56, 55, 54, 53, 52, 51, 50, 49,\n",
" 48, 47, 46, 45, 44, 43, 42, 41, 40, 39, 38, 37, 36, 35, 34, 33, 32,\n",
" 31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17],\n",
" dtype=int32)"
]
},
"execution_count": 132,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"d1[\"SnapNum\"]"
]
},
{
"cell_type": "code",
"execution_count": 127,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"8156"
]
},
"execution_count": 127,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.sum(cat[\"totmass\"] > 1e13)"
]
},
{
"cell_type": "code",
"execution_count": 128,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"33.333333333333336"
]
},
"execution_count": 128,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"250e-3 * 8000 / 60"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 116,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAi0AAAGiCAYAAAAr5/biAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABmw0lEQVR4nO3dd3xV9f3H8de92XvvQcIeQgIyRRE0FrFqxapYrVL81S6ctK72V2x/bbXWalGbih3OaovaOqqoFZShsveeAUL2vsnNurn3/P44JBjDCtzk5ibv5+NxH/fec86998MVkzffaTEMw0BERESkh7N6ugARERGRM6HQIiIiIl5BoUVERES8gkKLiIiIeAWFFhEREfEKCi0iIiLiFRRaRERExCsotIiIiIhXUGgRERERr6DQIiIiIl5BoUVERES8gkdCy8yZM4mKiuK6665rdzwjI4NRo0aRnZ3NtGnTPFGaiIiI9FAWT2yYuGzZMmpra3nppZd48803245nZGSwfft2QkNDu7skERER6eF8PfGhU6dOZdmyZW59T5fLRWFhIWFhYVgsFre+t4iIiHQNwzCora0lOTkZq/U0HUBGJy1fvty48sorjaSkJAMw3nrrrQ7X/PGPfzT69etnBAQEGOPHjzfWrFnT4ZpPP/3U+OY3v9nuWEZGhjFmzBhj7Nixxt///vdO1ZWfn28Auummm2666aabF97y8/NP+7u+0y0tdrudrKwsbrvtNq699toO5xctWsS8efNYuHAhEyZMYMGCBUyfPp09e/YQHx9/yvf+7LPPSElJoaioiJycHEaOHMmoUaNOeG1TUxNNTU1tz41jvVz5+fmEh4d39o8lIiIiHmCz2UhLSyMsLOy013Y6tMyYMYMZM2ac9PyTTz7J7bffzpw5cwBYuHAh77//Ps8//zwPPvjgKd87JSUFgKSkJK644go2btx40tDy6KOP8stf/rLD8fDwcIUWERERL3MmQzvcOnuoubmZDRs2kJOTc/wDrFZycnJYtWrVKV9rt9upra0FoK6ujk8++YQRI0ac9PqHHnqImpqatlt+fr57/hAiIiLSI7l1IG55eTlOp5OEhIR2xxMSEti9e3fb85ycHLZs2YLdbic1NZU33niDhIQEZs6cCYDT6eT2229n3LhxJ/2sgIAAAgIC3Fm+iIiI9GAemT20ZMmSEx7fsmVLN1ciIiIi3sKt3UOxsbH4+PhQUlLS7nhJSQmJiYnu/Kg2ubm5DB8+/JStMiIiIuL93Bpa/P39Of/881m6dGnbMZfLxdKlS5k0aZI7P6rN3Llz2blzJ+vWreuS9xcREZGeodPdQ3V1dezfv7/teV5eHps3byY6Opr09HTmzZvH7NmzGTt2LOPHj2fBggXY7fa22UQiIiIiZ6PToWX9+vXt9gWaN28eALNnz+bFF19k1qxZlJWVMX/+fIqLi8nOzubDDz/sMDhXREREpDM8svdQV7DZbERERFBTU6N1WkRERLxEZ35/e2SXZ3fSQFwREZG+QS0tIiIi4jF9qqVFRERE+gaFFhEREfEKHlkRV0RERLpHc4uL6oZmbA0tOJwuWpwGDpcLp8toe97icuFwms8dTheOFoMmp4smh5NGh5MGh5OGZhdp0UHMmZzpsT+L14eW3NxccnNzcTqdni5FRETkhBxOF7YGBzVfutkaW9qO2Rod1DW20NTiovnYzeF00eIycH7p1uJy4TTA6XLhdIFhGLgMA8PAvAcwwMAMKzUNDuqaWtz25xifGe3R0KKBuCIiIp3kdBmU1zVRamuivK711kx5XRMVX3rcGkrszZ79h7XFAmEBvvj7WvG1WvH1seDnY8XXasHXx4qfjwVfq3nM39eKv495TZCfD4HHbkH+PvSLDubG8elura0zv7+9vqVFRETEnQzDoLreQbGtkaKaBgqrGymsbjBvNebj4ppGWlyd/zd/WIAv4UF+hAf5ERHkS0SQH+GB5vOwQF8CfH3M0OBrxd/Hgo/VDBZWqxkqrBbz3ufYMasFrBYLFgtYMJ9bjj33tVqIDPYnKtj8DKvV0gXfVvdSaBERkT6lNZQcqrBzuKKe/Mp6jlY1cLS6noKqBopqGmlqcZ32fawWiA0NMG9hAcSG+Jv3of7EhAQQE+pPVLA/EUF+RBwLJb4+mv9yLhRaRESkV6qpd3CgvI5D5XYOVdRzqNzO4Qo7eeV2bI2nH+cRHeJPYnggyZFBpEQGkhQZdPxxRBDxYQEKId1MoUVERHqFIxX1vL25gM/2lXOgrI4Ke/Mpr08MD6RfTDDp0cGkRgWTGhVEalSQGUjCAwj08+mmyuVMeX1o0ewhEZG+q6Kuife2FvH25gI2HanucD4xPJDM2BAyYkPIiAmmX0wImbEhpEcHE+SvUOJtNHtIRES8Sn1zC//dUcLbmwtYua8c57EBsVYLTB4Yy9dHJjEiOYLMuBBCA7z+3+a9nmYPiYhIr5NfWc+LXxxi0br8dmuPjEqN4BvZKVw1Kon48EAPVti71TvqqXPUER8c77EaFFpERKRHcrkM8irs7Cy08cH2Ij7cXkzrLON+McF8IzuFb2QnMyAu1LOF9jKNLY3k1+Zz2HaYfVX72Fu1l71Ve8mvzWda2jSeuuQpj9Wm0CIiIh7X1OJkb3EdO4tq2FFoY0ehjV1FNuq/sijbRYNiue3CTC4eFNcr1h3pbs3OZsoayiirL6OsoYzS+lLKG8oprS+lpL6EI7YjFNuLOba2bgflDeXdXHF7Ci0iItKtDMNgd3Etn+8vZ2ehjZ1FNvaX1p1wsbZAPytDE8PJTovkxvFpDE3UmMXTsTvs7KvaR15NHkdqj5Bfm09+bT6FdYVUN1Wf0XuE+YXRL7wfA6MGMihyEIOjBzMochAxQTFdW/xpKLSIiEiXa2pxsuZgJUt2lbB0VykF1Q0drokM9mNEcjgjkiMYnhTOiORwMmNDtBbKl7gMF3WOOmqba6lqrKLYXkyxvZgiexH5tfnsq9rH0bqjp3wPP6sf8cHxxAbFEh8cT1xQHHHBccQFxZEenk6/8H5EBURhsfS8liyvDy2a8iwi0jOV1Taxcl8ZS3aVsHxPWbv9dwL9rFwwIJbstEiGJ4UzPDmcpIjAHvmLsqs5XU4K7YXk2/IpbSiltN68VTRUYGu2Udtci63Zhq3ZRl1z3Um7br4sLiiOgZEDSQ9PJy0sjdSwVFJDU0kITiAiIMJrv2dNeRYREbcor2tizcFKVh0sZ/XBSvaX1rU7Hx8WwKXD4skZlsDkgbG9fvG2mqYayhvKqWmqwdZso6apxrw112BrslHRWMEh2yEO1xym2XXqhfC+KtAnkPCAcBJDEkkMTiQxJJHk0GQGRQ5iUNQgogKjuuhP5X6a8iwiIl2u0t7MmoMVrDpYweqDFewtaR9SLBYYlhhOzrB4Lh2WwMiUiD4xeLagroCnNz7N4rzFZ/waP6sfaWFpJIYkEh8c39Z9E+EfQXhAOGH+YYT7H7/39/Hvwj8B4HJBbRFUHYLqw1B12LyP7g8X39+1n30KCi0iInJahmFwtKqBjUeq2HSkmtUHK9hdXNvhuqGJYUwaEMPE/jFMyIwmMriLf7l2I8MwOGQ7xOqi1awqXMWGkg2E+YcxMWkiE5MnMiJ6BIv2LOK13a/hcDkAiAiIMIOHfzgRAWYAaX0cGRBJv/B+ZIZnkhyajI/VQy1PTXVw+HM4/AVU7IeKA1CVBy2NHa9NHafQIiIiPYfLZXC4sv7YzB5zCvL2ghrK6zp2YQxJCGNi/2gmDYhhQmYMUSG9J6QAVDRUsKZoDauKVrG6aDXF9uJ2523NNv6171/8a9+/2h2fkDSBeefPY3jM8O4s98y4nFC4CQ58Cgc/hfy1cCxktWP1hYhUiOwHUf3M+/hh3V/vlyi0iIj0YS6XwcFyO5vzq9l6tJqdx9ZHsTd3nNzg52NheHIEo9MiGZ8ZzfjMaGJDAzxQtftVNVaxsXQjOyt2UlhXaN7shR1Cip/Vj9Hxo5mUPInxieOpbqpmVaEZaPZX72dQ1CDmnT+PycmTe85gV5cLyveYLSkHP4W8FdBY0/6ayH7Q/2JIOA+iB0B0JkS
"text/plain": [
"<Figure size 640x480 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure()\n",
"plt.plot(d1[\"Age\"], d1[\"MainProgenitorMass\"])\n",
"plt.plot(d2[\"Age\"], d2[\"MainProgenitorMass\"])\n",
"plt.plot(d3[\"Age\"], d3[\"MainProgenitorMass\"])\n",
"plt.yscale(\"log\")\n",
"\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 110,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1.4623428e+15, 7.4588315e+14, 7.1264091e+14, ..., 9.9620782e+10,\n",
" 9.9620782e+10, 9.9620782e+10], dtype=float32)"
]
},
"execution_count": 110,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cat[\"totmass\"]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "venv_csiborg",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}