{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cdb35308-1369-4ec0-a468-d104e0ab73cc",
   "metadata": {},
   "source": [
    "## Extract the efficiency of the matching between the shear catalogs\n",
    "\n",
    "Since they all use the same source detection, the only differences should come from differences in the applied noise"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86b99641-5ecc-4a7f-966a-cf2f085803ba",
   "metadata": {},
   "source": [
    "#### Standard imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2361a95-bb8f-4f07-8c48-cab604d54007",
   "metadata": {},
   "outputs": [],
   "source": [
    "import tables_io\n",
    "import numpy as np\n",
    "import hpmcm\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c516209f-845a-40d4-9146-54599e3480b3",
   "metadata": {},
   "source": [
    "#### Set up the configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "157d27ac-8b70-4cd2-b85c-566f8a826b8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "keys = ['_cluster_stats', '_cluster_shear']  # which tables to read\n",
    "st_ = 'wmom'          # which catalog type\n",
    "shear_st_ = \"0p01\"    # Applied shear as a string\n",
    "tract = 10463         # which tract to study\n",
    "\n",
    "dd = tables_io.read(f\"test_data/shear_{st_}_{shear_st_}_match_{tract}.pq\", keys=keys)\n",
    "data = dd['_cluster_stats']\n",
    "data2 = dd['_cluster_shear']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "423949ef-88d7-4a3c-a300-07fbb32e4c9a",
   "metadata": {},
   "source": [
    "#### Merge the two tables we read"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a592e9f0-9eaa-4514-ab09-e3794ba2db2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "data[\"idx\"] = np.arange(len(data))\n",
    "data2[\"idx\"] = np.arange(len(data2))\n",
    "merged = data.merge(data2, on=\"idx\")\n",
    "merged['has_ref_cat'] = merged.n_ns > 0\n",
    "merged['central'] = np.bitwise_and(\n",
    "    np.fabs(merged.x_cent-100) < 75,\n",
    "    np.fabs(merged.y_cent-100) < 75,\n",
    ")\n",
    "central = merged[merged.central]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d0c408f-f50a-4287-aafe-6bb2ef901871",
   "metadata": {},
   "source": [
    "#### Make maskes of different types of matches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d10da659-c31c-49b1-83a4-7cc75ed7b850",
   "metadata": {},
   "outputs": [],
   "source": [
    "good_mask = np.bitwise_and(central.n_src ==5, central.n_unique ==5)\n",
    "missing_md = np.bitwise_and(~good_mask, central.has_ref_cat) \n",
    "missing_ref = np.bitwise_and(~good_mask, ~central.has_ref_cat)\n",
    "extra = central.n_src > central.n_unique"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4b97179-d8b8-44d5-a941-9eca19386efc",
   "metadata": {},
   "source": [
    "#### Make a histogram of the different match types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "274782f6-fe38-4e9f-8e1b-9c89f0a091b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = plt.hist(central.iloc[good_mask.values].snr, bins=np.logspace(0, 5, 101), alpha=0.5, label=\"Good Match\")\n",
    "_ = plt.hist(central.iloc[missing_md.values].snr, bins=np.logspace(0, 5, 101), alpha=0.5, label=\"Has NS\")\n",
    "_ = plt.hist(central.iloc[missing_ref.values].snr, bins=np.logspace(0, 5, 101), alpha=0.5, label=\"Does not have NS\")\n",
    "_ = plt.hist(central.iloc[extra.values].snr, bins=np.logspace(0, 5, 101), alpha=0.5, label=\"Has confusion\")\n",
    "_ = plt.axvline(10, label='SN Cut')\n",
    "_ = plt.xscale('log')\n",
    "_ = plt.yscale('log')\n",
    "_ = plt.legend()\n",
    "_ = plt.xlabel(\"Signal-to-noise [r-band]\")\n",
    "_ = plt.ylabel(\"Objects [per 0.05 dex]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3834b9eb-cbd5-4440-8227-c5d4a16115e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "hist_all =  np.histogram(central.snr, bins=np.logspace(0, 5, 101))[0]\n",
    "hist_missing_md = np.histogram(central.iloc[missing_md.values].snr, bins=np.logspace(0, 5, 101))[0]\n",
    "hist_missing_ref = np.histogram(central.iloc[missing_ref.values].snr, bins=np.logspace(0, 5, 101))[0]\n",
    "hist_extra = np.histogram(central.iloc[extra.values].snr, bins=np.logspace(0, 5, 101))[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "516f6869-03da-4aa1-b251-302cfc2a2d8f",
   "metadata": {},
   "source": [
    "#### Estimate the good match efficiency as a function of SNR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f333a79-85fd-4117-8f55-c3c2be2212ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "ineffic_missing_md = hist_missing_md/hist_all\n",
    "ineffic_missing_ref = hist_missing_ref/hist_all\n",
    "ineffic_extra = hist_extra/hist_all\n",
    "npq_missing_md = np.sqrt(ineffic_missing_md*(1-ineffic_missing_md)/hist_all)\n",
    "npq_missing_ref = np.sqrt(ineffic_missing_ref*(1-ineffic_missing_ref)/hist_all)\n",
    "npq_missing_extra = np.sqrt(ineffic_extra*(1-ineffic_extra)/hist_all)\n",
    "bin_edges = np.logspace(0, 5, 101)\n",
    "bin_centers = np.sqrt(bin_edges[0:-1] * bin_edges[1:])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edb9e0dc-69c1-4b44-b268-dda3dd2326b0",
   "metadata": {},
   "source": [
    "#### Plot the good match efficiency as a function of SNR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "158c7e28-fe3f-43f2-874a-8f66dbda3d07",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = plt.errorbar(bin_centers, ineffic_missing_md, yerr=npq_missing_md, label=\"Has NS\", ls=\"\", marker='.')\n",
    "_ = plt.errorbar(bin_centers, ineffic_missing_ref, yerr=npq_missing_ref, label=\"Does not have NS\", ls=\"\", marker='.')\n",
    "_ = plt.xscale('log')\n",
    "_ = plt.yscale('log')\n",
    "_ = plt.ylim(1e-3,1)\n",
    "_ = plt.axvline(10, label='SN Cut')\n",
    "_ = plt.xlabel(\"Signal-to-noise [r-band]\")\n",
    "_ = plt.ylabel(\"Match Inefficieny\")\n",
    "_ = plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db11f6a6-8372-41b1-b275-edd294f6aa7e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de0fbce0-f0ce-4083-abfb-2f1796f3ccb3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}