Infusion_postprocessing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f2482bb0-6aef-409f-9491-70e47881efa7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from astropy.table import Table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c1cc066-b7aa-49ad-9764-011b6f343194",
   "metadata": {},
   "source": [
    "We use the constants for different IA models as inputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "628e5cc6-2685-4151-b4bc-bcdec02f1d05",
   "metadata": {},
   "outputs": [],
   "source": [
    "A_IA = 1.0\n",
    "b_TA = 0.0\n",
    "C_2  = 0.0\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d298a125-6c33-404d-89ec-58598ab919ab",
   "metadata": {},
   "source": [
    "The intensity of the NLA, delta-NLA and TT models can be chosen by tweaking the numbers above, this can probably be changed for some calculation involving weights given to each model but that's up for discussion, for now, we setup for the observed epsilons with and without noise, using the provided parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f0f13178-07aa-402f-a862-cf2f7a8e253a",
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-5-4ff29c0d13b0>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     13\u001b[0m     \u001b[0;31m#catalog = Table.read('../../GalCat/StageIV_nz/V1/GalCat_tomo{}_All_0.6GpAM_RA_Dec_g1_g2_w_z_den_Sij_linear_bias_smooth.asc'.format(i),\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m     \u001b[0;31m#catalog = Table.read('/global/homes/j/jharno/IA-infusion/SkySim5000/GalCat/StageIV_nz/GalCat_tomo{}_All_0.06GpAM_RA_Dec_g1_g2_w_den_Sij_smooth.dat'.format(i),\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 15\u001b[0;31m                      format='ascii')  # descriptions=['ra_deg', 'dec_deg', 'shear1', 'shear2', 'weight', 'z_true', 'delta', 's11', 's22', 's12'])\n\u001b[0m\u001b[1;32m     16\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     17\u001b[0m     \u001b[0;31m#V0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/lsst/software/stack/python/miniconda3-4.7.10/envs/lsst-scipipe-4d7b902/lib/python3.7/site-packages/astropy/table/connect.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m     50\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     51\u001b[0m         \u001b[0mcls\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cls\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 52\u001b[0;31m         \u001b[0mout\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mregistry\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcls\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     53\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     54\u001b[0m         \u001b[0;31m# For some readers (e.g., ascii.ecsv), the returned `out` class is not\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/lsst/software/stack/python/miniconda3-4.7.10/envs/lsst-scipipe-4d7b902/lib/python3.7/site-packages/astropy/io/registry.py\u001b[0m in \u001b[0;36mread\u001b[0;34m(cls, format, *args, **kwargs)\u001b[0m\n\u001b[1;32m    521\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    522\u001b[0m         \u001b[0mreader\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_reader\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcls\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 523\u001b[0;31m         \u001b[0mdata\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mreader\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    524\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    525\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcls\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/lsst/software/stack/python/miniconda3-4.7.10/envs/lsst-scipipe-4d7b902/lib/python3.7/site-packages/astropy/io/ascii/connect.py\u001b[0m in \u001b[0;36mio_read\u001b[0;34m(format, filename, **kwargs)\u001b[0m\n\u001b[1;32m     17\u001b[0m         \u001b[0mformat\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mre\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msub\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mr'^ascii\\.'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m''\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mformat\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     18\u001b[0m         \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'format'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mformat\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 19\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mread\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     20\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     21\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/lsst/software/stack/python/miniconda3-4.7.10/envs/lsst-scipipe-4d7b902/lib/python3.7/site-packages/astropy/io/ascii/ui.py\u001b[0m in \u001b[0;36mread\u001b[0;34m(table, guess, **kwargs)\u001b[0m\n\u001b[1;32m    285\u001b[0m         \u001b[0;31m# through below to the non-guess way so that any problems result in a\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    286\u001b[0m         \u001b[0;31m# more useful traceback.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 287\u001b[0;31m         \u001b[0mdat\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_guess\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtable\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnew_kwargs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mformat\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfast_reader\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    288\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mdat\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    289\u001b[0m             \u001b[0mguess\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/lsst/software/stack/python/miniconda3-4.7.10/envs/lsst-scipipe-4d7b902/lib/python3.7/site-packages/astropy/io/ascii/ui.py\u001b[0m in \u001b[0;36m_guess\u001b[0;34m(table, read_kwargs, format, fast_reader)\u001b[0m\n\u001b[1;32m    444\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    445\u001b[0m             \u001b[0mreader\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mguessing\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 446\u001b[0;31m             \u001b[0mdat\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mreader\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtable\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    447\u001b[0m             _read_trace.append({'kwargs': copy.deepcopy(guess_kwargs),\n\u001b[1;32m    448\u001b[0m                                 \u001b[0;34m'Reader'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mreader\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__class__\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/lsst/software/stack/python/miniconda3-4.7.10/envs/lsst-scipipe-4d7b902/lib/python3.7/site-packages/astropy/io/ascii/core.py\u001b[0m in \u001b[0;36mread\u001b[0;34m(self, table)\u001b[0m\n\u001b[1;32m   1168\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1169\u001b[0m         \u001b[0;31m# Extract table meta values (e.g. keywords, comments, etc).  Updates self.meta.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1170\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mheader\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mupdate_meta\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlines\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmeta\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1171\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1172\u001b[0m         \u001b[0;31m# Get the table column definitions\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/lsst/software/stack/python/miniconda3-4.7.10/envs/lsst-scipipe-4d7b902/lib/python3.7/site-packages/astropy/io/ascii/core.py\u001b[0m in \u001b[0;36mupdate_meta\u001b[0;34m(self, lines, meta)\u001b[0m\n\u001b[1;32m    547\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcomment\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    548\u001b[0m             \u001b[0mre_comment\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mre\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcomment\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 549\u001b[0;31m             \u001b[0mcomment_lines\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mx\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mx\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mlines\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mre_comment\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmatch\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    550\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    551\u001b[0m             \u001b[0mcomment_lines\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/lsst/software/stack/python/miniconda3-4.7.10/envs/lsst-scipipe-4d7b902/lib/python3.7/site-packages/astropy/io/ascii/core.py\u001b[0m in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m    547\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcomment\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    548\u001b[0m             \u001b[0mre_comment\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mre\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcomment\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 549\u001b[0;31m             \u001b[0mcomment_lines\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mx\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mx\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mlines\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mre_comment\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmatch\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    550\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    551\u001b[0m             \u001b[0mcomment_lines\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "for i in range(1,6):\n",
    "#for i in np.arange(5)+1:\n",
    "\n",
    "    #V0\n",
    "    #catalog = Table.read('../../GalCat/StageIV_nz/V0/GalCat_tomo{}_All_0.06GpAM_RA_Dec_g1_g2_w_den_Sij.dat'.format(i),\n",
    "    #catalog = Table.read('../../GalCat/StageIV_nz/V0/GalCat_tomo{}_All_0.06GpAM_RA_Dec_g1_g2_w_den_Sij_smooth.dat'.format(i),\n",
    "    catalog = Table.read('../../GalCat/SRD-Y1/random_pos/GalCat_tomo{}_All_0.6GpAM_RA_Dec_g1_g2_w_z_den_Sij_smooth0.5Mpc.asc'.format(i),\n",
    "\n",
    "    #V1\n",
    "    #catalog = Table.read('../../GalCat/StageIV_nz/V1/GalCat_tomo{}_All_0.06GpAM_RA_Dec_g1_g2_w_z_den_Sij_linear_bias.asc'.format(i),\n",
    "    #catalog = Table.read('../../GalCat/StageIV_nz/V1/GalCat_tomo{}_All_0.06GpAM_RA_Dec_g1_g2_w_z_den_Sij_linear_bias_smooth.asc'.format(i),\n",
    "    #catalog = Table.read('../../GalCat/StageIV_nz/V1/GalCat_tomo{}_All_0.06GpAM_RA_Dec_g1_g2_w_z_den_Sij_linear_bias_smooth_1Mpc.asc'.format(i),\n",
    "    #catalog = Table.read('../../GalCat/StageIV_nz/V1/GalCat_tomo{}_All_0.6GpAM_RA_Dec_g1_g2_w_z_den_Sij_linear_bias_smooth.asc'.format(i),\n",
    "    #catalog = Table.read('/global/homes/j/jharno/IA-infusion/SkySim5000/GalCat/StageIV_nz/GalCat_tomo{}_All_0.06GpAM_RA_Dec_g1_g2_w_den_Sij_smooth.dat'.format(i),\n",
    "                     format='ascii')  # descriptions=['ra_deg', 'dec_deg', 'shear1', 'shear2', 'weight', 'z_true', 'delta', 's11', 's22', 's12'])\n",
    "\n",
    "    #V0\n",
    "    #epsilons = Table.read('../../GalCat/StageIV_nz/V0/epsilon_{}_0.06GpAM.asc'.format(i), format='ascii')  #descriptions=['e1','e2','e1_noise','e2_noise']\n",
    "    #epsilons = Table.read('../../GalCat/StageIV_nz/V0/epsilon_{}_0.06GpAM_V2.asc'.format(i), format='ascii')  #descriptions=['e1','e2','e1_noise','e2_noise']\n",
    "    #epsilons = Table.read('../../GalCat/StageIV_nz/V0/epsilon_{}_0.06GpAM_smooth_V2.asc'.format(i), format='ascii')  #descriptions=['e1','e2','e1_noise','e2_noise']\n",
    "    #epsilons = Table.read('../../GalCat/SRD-Y1/random_pos/epsilon_{}_0.6GpAM_smooth0.5Mpc.asc'.format(i), format='ascii')  #descriptions=['e1','e2','e1_noise','e2_noise']\n",
    "    epsilons = Table.read('../../GalCat/SRD-Y1/random_pos/epsilon_{}_0.6GpAM_smooth0.5Mpc_nofactor2.asc'.format(i), format='ascii')  #descriptions=['e1','e2','e1_noise','e2_noise']\n",
    "\n",
    "    #V1\n",
    "    #epsilons = Table.read('../../GalCat/StageIV_nz/V1/epsilon_{}_0.06GpAM_linear_bias.asc'.format(i), format='ascii')  #descriptions=['e1','e2','e1_noise','e2_noise']\n",
    "    #epsilons = Table.read('../../GalCat/StageIV_nz/V1/epsilon_{}_0.06GpAM_linear_bias_smooth.asc'.format(i), format='ascii')  #descriptions=['e1','e2','e1_noise','e2_noise']\n",
    "    #epsilons = Table.read('../../GalCat/StageIV_nz/V1/epsilon_{}_0.6GpAM_linear_bias_smooth.asc'.format(i), format='ascii')  #descriptions=['e1','e2','e1_noise','e2_noise']\n",
    "    #epsilons = Table.read('../../GalCat/StageIV_nz/V1/epsilon_{}_0.06GpAM_linear_bias_smooth_1Mpc.asc'.format(i), format='ascii')  #descriptions=['e1','e2','e1_noise','e2_noise']\n",
    "\n",
    "    \n",
    "    ra=np.asarray(catalog['col1'])\n",
    "    dec=np.asarray(catalog['col2'])\n",
    "    weight=np.asarray(catalog['col5'])\n",
    "    z=np.asarray(catalog['col6'])\n",
    "    gamma1 = np.asarray(catalog['col3'])\n",
    "    gamma2 = np.asarray(catalog['col4'])\n",
    "    e1_NLA = np.asarray(epsilons['e1_NLA'])\n",
    "    e2_NLA = np.asarray(epsilons['e2_NLA'])\n",
    "    e1_TT = np.asarray(epsilons['e1_TT'])\n",
    "    e2_TT = np.asarray(epsilons['e2_TT'])\n",
    "    e1_noise = np.asarray(epsilons['e1_noise'])\n",
    "    e2_noise = np.asarray(epsilons['e2_noise'])\n",
    "    delta = np.asarray(catalog['col7'])\n",
    "    epsilon_NLA = np.empty(e1_NLA.shape, dtype = complex)\n",
    "    epsilon_NLA.real = A_IA * e1_NLA  \n",
    "    epsilon_NLA.imag = A_IA * e2_NLA  \n",
    "    epsilon_TT = np.empty(e1_TT.shape, dtype = complex)\n",
    "    epsilon_TT.real = C_2 * e1_TT  \n",
    "    epsilon_TT.imag = C_2 * e2_TT\n",
    "    epsilon_random = np.empty(e1_TT.shape, dtype = complex)\n",
    "    epsilon_random.real = e1_noise\n",
    "    epsilon_random.imag = e2_noise\n",
    "    g = np.empty(gamma1.shape, dtype=complex)\n",
    "    g.real = gamma1\n",
    "    g.imag = gamma2\n",
    "    \n",
    "    \n",
    "    #calibrate all IA:\n",
    "    # Needed after finding the factor of 2.0 bug in the e_IA2\n",
    "    epsilon_NLA*=1.5\n",
    "    epsilon_TT*=1.5\n",
    "    \n",
    "    #calibrate delta:\n",
    "    #delta /=2.5\n",
    "    #delta*=0.0\n",
    "    \n",
    "    #calibrate the low-redshift e_TT:\n",
    "    epsilon_TT[z<0.5]/=2.5\n",
    "\n",
    "    #V1 mocks further require an overall calibration (not sure why)\n",
    "    #epsilon_TT /=20.0\n",
    "\n",
    "    \n",
    "    #And finally, we calculate the observed epsilons using the summation as in equation 11 of the paper\n",
    "    # The V1 (linear bias) mocks already have the delta-weighted term with b_ta=1.0, so we keep it or we undo it. \n",
    "    # In the later case, we must use the same b_TA parameter as that used in the construction of the mocks\n",
    "    \n",
    "    # 1-Keep delta-weights:\n",
    "    epsilon_IA_TATT = epsilon_NLA + epsilon_TT\n",
    "    # 2-Undo delta-weights:\n",
    "    #epsilon_IA_TATT = epsilon_NLA / (1 + b_TA * delta) + epsilon_TT\n",
    "    #epsilon_IA_TATT /= (1 + b_TA * delta)\n",
    "    # 3-Change from b-TA = 1 to b_TA = 0 or 2.0 \n",
    "    #epsilon_IA_TATT *= (1 + b_TA * delta)/(1 + 1.0 * delta)\n",
    "    \n",
    "    epsilon_int = (epsilon_IA_TATT + epsilon_random)/(1 + np.conjugate(epsilon_IA_TATT) * epsilon_random)\n",
    "    epsilon_obs = (epsilon_int + g)/(1 + epsilon_int * np.conjugate(g))\n",
    "    epsilon_obs_no_noise = (epsilon_IA_TATT + g)/(1 + epsilon_IA_TATT * np.conjugate(g))\n",
    "   \n",
    "    fin_epsilon = Table([ra,dec,z,gamma1, gamma2, weight,delta,epsilon_obs.real,epsilon_obs.imag,epsilon_obs_no_noise.real,epsilon_obs_no_noise.imag,epsilon_IA_TATT.real,epsilon_IA_TATT.imag], names = ('ra','dec','z','g1', 'g2', 'weight','delta','e_obs_1','e_obs_2', 'e_obs_no_noise_1','e_obs_no_noise_2','e_IA_TATT_1', 'e_IA_TATT_2' ))\n",
    "    #fin_epsilon.write('/global/cscratch1/sd/jharno/fin_epsilon_IA_{}.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('/global/cscratch1/sd/jharno/fin_epsilon_IA_{}_linear_bias.dat'.format(i), format='ascii', overwrite=True)    \n",
    "\n",
    "    #V0:\n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V0/fin_epsilon_IA_{}_0.06GpAM.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/random_pos/fin_epsilon_IA_{}_0.6GpAM_smooth0.1Mpc.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/random_pos/fin_epsilon_IA_{}_0.6GpAM_smooth0.5Mpc_nofactor2.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/NLA/galcat_skysim5000_noIA_tomo{}.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    fin_epsilon.write('../../GalCat/SRD-Y1/NLA/galcat_skysim5000_IAIp1_tomo{}_nofactor2.dat'.format(i), format='ascii', overwrite=True)    \n",
    "\n",
    "    #V0-TT                     \n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V0/fin_epsilon_IA_{}_0.06GpAM_TT.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V0/fin_epsilon_IA_{}_0.06GpAM_TT_V2.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V0/fin_epsilon_IA_{}_0.06GpAM_TT_V2corr.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V0/fin_epsilon_IA_{}_0.06GpAM_TT_smooth_V2corr.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V0/fin_epsilon_IA_{}_0.06GpAM_TT_smooth_V3corr.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/random_pos/fin_epsilon_IA_{}_0.6GpAM_TT_smooth0.1Mpc.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/random_pos/fin_epsilon_IA_{}_0.6GpAM_TT_smooth0.5Mpc_m1.dat'.format(i), format='ascii', overwrite=True)    \n",
    "\n",
    "    #V1:\n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V1/fin_epsilon_IA_{}_0.06GpAM_linear_bias.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V1/fin_epsilon_IA_{}_0.06GpAM_linear_bias_smooth.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V1/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V1/fin_epsilon_IA_{}_0.06GpAM_linear_bias_smooth_1Mpc.dat'.format(i), format='ascii', overwrite=True)    \n",
    "\n",
    "    #V1-NLA:\n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V1/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth_NLA.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    \n",
    "    print('Done tomo',i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "561c04d0-8c88-40ce-8564-621080c88d32",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Done tomo 1\n",
      "Done tomo 2\n",
      "Done tomo 3\n",
      "Done tomo 4\n",
      "Done tomo 5\n"
     ]
    }
   ],
   "source": [
    "for i in range(1,6):\n",
    "    ra,dec,z,gamma1,gamma2 = np.loadtxt('../../GalCat/HOD/FromCharlie/GalCat_tomo'+str(i)+'_06GpAM_new.dat',usecols = (3,6,1,8,9), unpack=True)\n",
    "    epsilons = Table.read('../../GalCat/HOD/FromCharlie/epsilon_{}_0.6GpAM.asc'.format(i), format='ascii')\n",
    "\n",
    "    #ra=np.asarray(catalog[:,4])\n",
    "    #dec=np.asarray(catalog[:,7])\n",
    "    weight=np.asarray(np.ones(np.shape(ra)[0]))\n",
    "    #z=np.asarray(catalog[:,2])\n",
    "    #gamma1 = np.asarray(catalog[:,8])\n",
    "    #gamma2 = np.asarray(catalog[:,9])\n",
    "    e1_NLA = np.asarray(epsilons['e1_NLA'])\n",
    "    e2_NLA = np.asarray(epsilons['e2_NLA'])\n",
    "    e1_TT = np.asarray(epsilons['e1_TT'])\n",
    "    e2_TT = np.asarray(epsilons['e2_TT'])\n",
    "    e1_noise = np.asarray(epsilons['e1_noise'])\n",
    "    e2_noise = np.asarray(epsilons['e2_noise'])\n",
    "    #delta = np.asarray(catalog[:,6])\n",
    "    delta = np.asarray(np.ones(np.shape(ra)[0]))\n",
    "    epsilon_NLA = np.empty(e1_NLA.shape, dtype = complex)\n",
    "    epsilon_NLA.real = A_IA * e1_NLA *(-1.0)\n",
    "    epsilon_NLA.imag = A_IA * e2_NLA *(-1.0)\n",
    "    epsilon_TT = np.empty(e1_TT.shape, dtype = complex)\n",
    "    epsilon_TT.real = C_2 * e1_TT *(-1.0)\n",
    "    epsilon_TT.imag = C_2 * e2_TT *(-1.0)\n",
    "    epsilon_random = np.empty(e1_TT.shape, dtype = complex)\n",
    "    epsilon_random.real = e1_noise\n",
    "    epsilon_random.imag = e2_noise\n",
    "    g = np.empty(gamma1.shape, dtype=complex)\n",
    "    g.real = gamma1\n",
    "    g.imag = gamma2\n",
    "    \n",
    "    \n",
    "    #calibrate all IA:\n",
    "    # Needed after finding the factor of 2.0 bug in the e_IA2\n",
    "    epsilon_NLA*=1.5\n",
    "    epsilon_TT*=1.5\n",
    "    \n",
    "    \n",
    "    #------------------\n",
    "    #calibrate delta:\n",
    "    #delta*=0.0\n",
    "\n",
    "    #calibrate HOD_NLA:\n",
    "    epsilon_NLA/=np.sqrt(2.0) \n",
    "\n",
    "    #calibrate TT:\n",
    "    # for TT: \n",
    "    #epsilon_TT[z<0.5]/=2.5 \n",
    "   \n",
    "    \n",
    "    #calibrate HOD_TT:\n",
    "    # for TT: \n",
    "    #epsilon_TT[z<0.75]/=4.0 \n",
    "    #epsilon_TT/=7\n",
    "    \n",
    "    # for and delta-TT\n",
    "    #epsilon_TT[z<0.5]/=5.0 \n",
    "    #epsilon_TT/=7 \n",
    "    #------------------\n",
    "    \n",
    "    #And finally, we calculate the observed epsilons using the summation as in equation 11 of the paper\n",
    "    # The V1 mocks already sample the galaxy by linearly tracing the matter field, so they already have the \n",
    "    # delta-weighted term. \n",
    "    # 1-Keep as is:\n",
    "    epsilon_IA_TATT = epsilon_NLA + epsilon_TT\n",
    "    # 2-Undo delta-weight:\n",
    "    #epsilon_IA_TATT = epsilon_NLA / (1 + b_TA * delta) + epsilon_TT\n",
    "    #epsilon_IA_TATT /= (1 + b_TA * delta) \n",
    "    # 3-Change from b_TA = 1 to b_TA = 0 or 2.0 \n",
    "    #epsilon_IA_TATT *= (1 + b_TA * delta)/(1 + 1.0 * np.where(delta >=  -0.999, delta, -0.999))\n",
    "\n",
    "    \n",
    "    \n",
    "    epsilon_int = (epsilon_IA_TATT + epsilon_random)/(1 + np.conjugate(epsilon_IA_TATT) * epsilon_random)\n",
    "    epsilon_obs = (epsilon_int + g)/(1 + epsilon_int * np.conjugate(g))\n",
    "    epsilon_obs_no_noise = (epsilon_IA_TATT + g)/(1 + epsilon_IA_TATT * np.conjugate(g))\n",
    "   \n",
    "    fin_epsilon = Table([ra,dec,z,gamma1, gamma2, weight,delta,epsilon_obs.real,epsilon_obs.imag,epsilon_obs_no_noise.real,epsilon_obs_no_noise.imag,epsilon_IA_TATT.real,epsilon_IA_TATT.imag], names = ('#ra','dec','z','g1', 'g2', 'weight','delta','e_obs_1','e_obs_2', 'e_obs_no_noise_1','e_obs_no_noise_2','e_IA_TATT_1', 'e_IA_TATT_2' ))    \n",
    "    #fin_epsilon.write('../../GalCat/HOD/FromCharlie/fin_epsilon_IA_{}_0.6GpAM.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/HOD/FromCharlie/fin_epsilon_IA_{}_0.6GpAM_V2.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/HOD/FromCharlie/fin_epsilon_IA_{}_0.6GpAM_nofactor2.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    fin_epsilon.write('../../GalCat/HOD/FromCharlie/fin_epsilon_TT_{}_0.6GpAM_nofactor2.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/HOD/FromCharlie/fin_epsilon_IA_{}_0.6GpAM_TT.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    \n",
    "    \n",
    "    print('Done tomo',i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "3d93bf77-cafb-43b9-90d9-259592bbadfd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total 14267616\n",
      "drwxrwsr-x 2 jharno lsst     16384 Aug 20 05:39 \u001b[0m\u001b[01;34mcorrelations\u001b[0m/\n",
      "-rwxrwxr-x 1 jharno lsst 213455758 Oct 26 07:09 \u001b[01;32mepsilon_1_0.6GpAM.asc\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 213726492 Oct 26 07:09 \u001b[01;32mepsilon_2_0.6GpAM.asc\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 213922638 Oct 26 07:09 \u001b[01;32mepsilon_3_0.6GpAM.asc\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 174726993 Oct 26 07:09 \u001b[01;32mepsilon_4_0.6GpAM.asc\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 100204130 Oct 26 07:10 \u001b[01;32mepsilon_5_0.6GpAM.asc\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 383792502 Feb 28  2024 \u001b[01;32mfin_epsilon_IA_1_0.6GpAM_AIA2.dat\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 383101720 Feb 27  2024 \u001b[01;32mfin_epsilon_IA_1_0.6GpAM.dat\u001b[0m*\n",
      "-rw-r--r-- 1 jharno lsst 383719698 Oct 28 00:17 fin_epsilon_IA_1_0.6GpAM_nofactor2.dat\n",
      "-rwxrwxr-x 1 jharno lsst 385968199 Mar  7  2024 \u001b[01;32mfin_epsilon_IA_1_0.6GpAM_TT.dat\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 383795949 Aug 20 05:23 \u001b[01;32mfin_epsilon_IA_1_0.6GpAM_V2.dat\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 382820211 Feb 28  2024 \u001b[01;32mfin_epsilon_IA_2_0.6GpAM_AIA2.dat\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 382297036 Feb 27  2024 \u001b[01;32mfin_epsilon_IA_2_0.6GpAM.dat\u001b[0m*\n",
      "-rw-r--r-- 1 jharno lsst 382760597 Oct 28 00:18 fin_epsilon_IA_2_0.6GpAM_nofactor2.dat\n",
      "-rwxrwxr-x 1 jharno lsst 384591351 Mar  7  2024 \u001b[01;32mfin_epsilon_IA_2_0.6GpAM_TT.dat\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 382820607 Aug 20 05:23 \u001b[01;32mfin_epsilon_IA_2_0.6GpAM_V2.dat\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 382340152 Feb 28  2024 \u001b[01;32mfin_epsilon_IA_3_0.6GpAM_AIA2.dat\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 381886998 Feb 27  2024 \u001b[01;32mfin_epsilon_IA_3_0.6GpAM.dat\u001b[0m*\n",
      "-rw-r--r-- 1 jharno lsst 382285502 Oct 28 00:19 fin_epsilon_IA_3_0.6GpAM_nofactor2.dat\n",
      "-rwxrwxr-x 1 jharno lsst 383853302 Mar  7  2024 \u001b[01;32mfin_epsilon_IA_3_0.6GpAM_TT.dat\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 382342607 Aug 20 05:23 \u001b[01;32mfin_epsilon_IA_3_0.6GpAM_V2.dat\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 358437260 Feb 28  2024 \u001b[01;32mfin_epsilon_IA_4_0.6GpAM_AIA2.dat\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 358752197 Feb 27  2024 \u001b[01;32mfin_epsilon_IA_4_0.6GpAM.dat\u001b[0m*\n",
      "-rw-r--r-- 1 jharno lsst 358402803 Oct 28 00:19 fin_epsilon_IA_4_0.6GpAM_nofactor2.dat\n",
      "-rwxrwxr-x 1 jharno lsst 359153031 Mar  7  2024 \u001b[01;32mfin_epsilon_IA_4_0.6GpAM_TT.dat\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 358436776 Aug 20 05:24 \u001b[01;32mfin_epsilon_IA_4_0.6GpAM_V2.dat\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 314112665 Feb 28  2024 \u001b[01;32mfin_epsilon_IA_5_0.6GpAM_AIA2.dat\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 315748896 Feb 27  2024 \u001b[01;32mfin_epsilon_IA_5_0.6GpAM.dat\u001b[0m*\n",
      "-rw-r--r-- 1 jharno lsst 314116096 Oct 28 00:20 fin_epsilon_IA_5_0.6GpAM_nofactor2.dat\n",
      "-rwxrwxr-x 1 jharno lsst 314118327 Mar  7  2024 \u001b[01;32mfin_epsilon_IA_5_0.6GpAM_TT.dat\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 314117964 Aug 20 05:24 \u001b[01;32mfin_epsilon_IA_5_0.6GpAM_V2.dat\u001b[0m*\n",
      "-rw-r--r-- 1 jharno lsst 382749828 Oct 28 00:26 fin_epsilon_TT_1_0.6GpAM_nofactor2.dat\n",
      "-rw-r--r-- 1 jharno lsst 381894449 Oct 28 00:27 fin_epsilon_TT_2_0.6GpAM_nofactor2.dat\n",
      "-rw-r--r-- 1 jharno lsst 381618638 Oct 28 00:28 fin_epsilon_TT_3_0.6GpAM_nofactor2.dat\n",
      "-rw-r--r-- 1 jharno lsst 358196748 Oct 28 00:28 fin_epsilon_TT_4_0.6GpAM_nofactor2.dat\n",
      "-rw-r--r-- 1 jharno lsst 314114740 Oct 28 00:29 fin_epsilon_TT_5_0.6GpAM_nofactor2.dat\n",
      "-rwxrwxr-x 1 jharno lsst 573971878 Feb 22  2024 \u001b[01;32mGalCat_tomo1_06GpAM_new.dat\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 572894327 Feb 22  2024 \u001b[01;32mGalCat_tomo2_06GpAM_new.dat\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 572622167 Feb 22  2024 \u001b[01;32mGalCat_tomo3_06GpAM_new.dat\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 545852349 Feb 22  2024 \u001b[01;32mGalCat_tomo4_06GpAM_new.dat\u001b[0m*\n",
      "-rwxrwxr-x 1 jharno lsst 495959998 Feb 22  2024 \u001b[01;32mGalCat_tomo5_06GpAM_new.dat\u001b[0m*\n"
     ]
    }
   ],
   "source": [
    "ls -l ../../GalCat/HOD/FromCharlie/\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f534a7ed-2cd8-4d2d-9768-2303df6eac9a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Done tomo 1\n",
      "Done tomo 2\n",
      "Done tomo 3\n"
     ]
    }
   ],
   "source": [
    "for i in range(1,6):\n",
    "#for i in np.arange(5)+1:\n",
    "\n",
    "    #SkySim5000\n",
    "    #V0\n",
    "    #catalog = np.loadtxt('../../GalCat/SRD-Y1/random_pos/GalCat_tomo'+str(i)+'_All_0.6GpAM_RA_Dec_g1_g2_w_z_den_Sij_smooth0.5Mpc.asc')\n",
    "    \n",
    "    #V1\n",
    "    #catalog = Table.read('../../GalCat/StageIV_nz/V1/GalCat_tomo{}_All_0.06GpAM_RA_Dec_g1_g2_w_z_den_Sij_linear_bias.asc'.format(i),\n",
    "    #catalog = Table.read('../../GalCat/StageIV_nz/V1/GalCat_tomo{}_All_0.06GpAM_RA_Dec_g1_g2_w_z_den_Sij_linear_bias_smooth.asc'.format(i),\n",
    "    #catalog = Table.read('../../GalCat/StageIV_nz/V1/GalCat_tomo{}_All_0.06GpAM_RA_Dec_g1_g2_w_z_den_Sij_linear_bias_smooth_1Mpc.asc'.format(i),\n",
    "    #catalog = np.loadtxt('../../GalCat/StageIV_nz/V1/GalCat_tomo'+str(i)+'_All_0.6GpAM_RA_Dec_g1_g2_w_z_den_Sij_linear_bias_smooth.asc')\n",
    "    #catalog = np.loadtxt('../../GalCat/StageIV_nz/V1/GalCat_tomo'+str(i)+'_All_0.6GpAM_RA_Dec_g1_g2_w_z_den_Sij_linear_bias_smooth_bta2.asc')\n",
    "    #catalog = np.loadtxt('../../GalCat/SRD-Y1/lin_bias_pos/GalCat_tomo'+str(i)+'_All_0.6GpAM_RA_Dec_g1_g2_w_z_den_Sij_linear_bias_bta1_smooth0.5Mpc.asc')\n",
    "    #catalog = np.loadtxt('../../GalCat/SRD-Y1/lin_bias_pos/GalCat_tomo'+str(i)+'_All_0.6GpAM_RA_Dec_g1_g2_w_z_den_Sij_linear_bias_bta2_smooth0.1Mpc.asc')\n",
    "    #catalog = np.loadtxt('../../GalCat/SRD-Y1/lin_bias_pos/GalCat_tomo'+str(i)+'_All_0.6GpAM_RA_Dec_g1_g2_w_z_den_Sij_linear_bias_smooth1_bta1_smooth0.5Mpc.asc')\n",
    "    #catalog = np.loadtxt('../../GalCat/SRD-Y1/lin_bias_pos/GalCat_tomo'+str(i)+'_All_0.6GpAM_RA_Dec_g1_g2_w_z_den_Sij_linear_bias_smooth1_bta2_smooth0.5Mpc.asc')\n",
    "    \n",
    "    #V0\n",
    "    #epsilons = Table.read('../../GalCat/StageIV_nz/V0/epsilon_{}_0.06GpAM.asc'.format(i), format='ascii')  #descriptions=['e1','e2','e1_noise','e2_noise']\n",
    "    #epsilons = Table.read('../../GalCat/SRD-Y1/random_pos/epsilon_{}_0.6GpAM_smooth0.5Mpc.asc'.format(i), format='ascii')  #descriptions=['e1','e2','e1_noise','e2_noise']\n",
    "    #epsilons = Table.read('../../GalCat/SRD-Y1/random_pos/epsilon_'+str(i)+'_0.6GpAM_smooth0.5Mpc.asc', format='ascii')  #descriptions=['e1','e2','e1_noise','e2_noise']\n",
    "    #epsilons = Table.read('../../GalCat/SRD-Y1/random_pos/epsilon_'+str(i)+'_0.6GpAM_smooth0.5Mpc_nofactor2.asc', format='ascii')  #descriptions=['e1','e2','e1_noise','e2_noise']\n",
    "\n",
    "    #V1\n",
    "    #epsilons = Table.read('../../GalCat/StageIV_nz/V1/epsilon_{}_0.06GpAM_linear_bias.asc'.format(i), format='ascii')  #descriptions=['e1','e2','e1_noise','e2_noise']\n",
    "    #epsilons = Table.read('../../GalCat/StageIV_nz/V1/epsilon_{}_0.06GpAM_linear_bias_smooth.asc'.format(i), format='ascii')  #descriptions=['e1','e2','e1_noise','e2_noise']\n",
    "    #epsilons = Table.read('../../GalCat/StageIV_nz/V1/epsilon_{}_0.6GpAM_linear_bias_smooth.asc'.format(i), format='ascii')  #descriptions=['e1','e2','e1_noise','e2_noise']\n",
    "    #epsilons = Table.read('../../GalCat/StageIV_nz/V1/epsilon_{}_0.6GpAM_linear_bias_smooth_bta2.asc'.format(i), format='ascii')  #descriptions=['e1','e2','e1_noise','e2_noise']\n",
    "    #epsilons = np.loadtxt('../../GalCat/StageIV_nz/V1/epsilon_'+str(i)+'_0.6GpAM_linear_bias_smooth.asc')  #descriptions=['e1','e2','e1_noise','e2_noise']\n",
    "    #epsilons = Table.read('../../GalCat/StageIV_nz/V1/epsilon_{}_0.06GpAM_linear_bias_smooth_1Mpc.asc'.format(i), format='ascii')  #descriptions=['e1','e2','e1_noise','e2_noise']    \n",
    "    #epsilons = Table.read('../../GalCat/SRD-Y1/lin_bias_pos/epsilon_{}_0.6GpAM_smooth0.5Mpc.asc'.format(i), format='ascii')\n",
    "    #epsilons = Table.read('../../GalCat/SRD-Y1/lin_bias_pos/epsilon_{}_0.6GpAM_smooth_bta2_0.1Mpc.asc'.format(i), format='ascii')\n",
    "    #epsilons = Table.read('../../GalCat/SRD-Y1/lin_bias_pos/epsilon_{}_0.6GpAM_linear_bias_smooth1_bta1_smooth05Mpc_nofactor2.asc'.format(i), format='ascii')\n",
    "    #epsilons = Table.read('../../GalCat/SRD-Y1/lin_bias_pos/epsilon_{}_0.6GpAM_linear_bias_smooth1_bta2_smooth05Mpc.asc'.format(i), format='ascii')\n",
    "\n",
    "    #BCM:\n",
    "    catalog = np.loadtxt('../../../BCM/dmo/GalCat/GalCat_tomo'+str(i)+'_All_0.6GpAM_RA_Dec_g1_g2_w_z_den_Sij_smooth0.5Mpc.asc')\n",
    "    epsilons = Table.read('../../../BCM/dmo/GalCat/epsilon_{}_0.6GpAM_smooth05Mpc.asc'.format(i), format='ascii')\n",
    "   \n",
    "    \n",
    "    ra=np.asarray(catalog[:,0])\n",
    "    dec=np.asarray(catalog[:,1])\n",
    "    weight=np.asarray(catalog[:,4])\n",
    "    z=np.asarray(catalog[:,5])\n",
    "    gamma1 = np.asarray(catalog[:,2])\n",
    "    gamma2 = np.asarray(catalog[:,3])\n",
    "    e1_NLA = np.asarray(epsilons['e1_NLA'])\n",
    "    e2_NLA = np.asarray(epsilons['e2_NLA'])\n",
    "    e1_TT = np.asarray(epsilons['e1_TT'])\n",
    "    e2_TT = np.asarray(epsilons['e2_TT'])\n",
    "    e1_noise = np.asarray(epsilons['e1_noise'])\n",
    "    e2_noise = np.asarray(epsilons['e2_noise'])\n",
    "    delta = np.asarray(catalog[:,6])\n",
    "    epsilon_NLA = np.empty(e1_NLA.shape, dtype = complex)\n",
    "    epsilon_NLA.real = A_IA * e1_NLA  \n",
    "    epsilon_NLA.imag = A_IA * e2_NLA  \n",
    "    epsilon_TT = np.empty(e1_TT.shape, dtype = complex)\n",
    "    epsilon_TT.real = C_2 * e1_TT  \n",
    "    epsilon_TT.imag = C_2 * e2_TT\n",
    "    epsilon_random = np.empty(e1_TT.shape, dtype = complex)\n",
    "    epsilon_random.real = e1_noise\n",
    "    epsilon_random.imag = e2_noise\n",
    "    g = np.empty(gamma1.shape, dtype=complex)\n",
    "    g.real = gamma1\n",
    "    g.imag = gamma2\n",
    "    \n",
    "    \n",
    "    #calibrate all IA:\n",
    "    # Needed after finding the factor of 2.0 bug in the e_IA2\n",
    "    epsilon_NLA*=1.5\n",
    "    epsilon_TT*=1.5\n",
    "    \n",
    "    \n",
    "    #------------------\n",
    "    #calibrate delta:\n",
    "    #delta*=0.0\n",
    "    \n",
    "    #calibrate TT:\n",
    "    # for TT: \n",
    "    epsilon_TT[z<0.5]/=2.5 \n",
    "    \n",
    "    # calibrate delta-TT\n",
    "    #epsilon_TT[z<0.5]/=5.0 \n",
    "    #epsilon_TT/=7 \n",
    "    #------------------\n",
    "    \n",
    "    #And finally, we calculate the observed epsilons using the summation as in equation 11 of the paper\n",
    "    # The V1 mocks already sample the galaxy by linearly tracing the matter field, so they already have the \n",
    "    # delta-weighted term. \n",
    "    # 1-Keep as is:\n",
    "    epsilon_IA_TATT = epsilon_NLA + epsilon_TT\n",
    "    # 2-Undo delta-weight:\n",
    "    #epsilon_IA_TATT = epsilon_NLA / (1 + b_TA * delta) + epsilon_TT\n",
    "    #epsilon_IA_TATT /= (1 + b_TA * delta) \n",
    "    # 3-Change from b_TA = 1 to b_TA = 0 or 2.0, uncomment the following: \n",
    "    #epsilon_IA_TATT *= (1 + b_TA * delta)/(1 + 1.0 * np.where(delta >=  -0.999, delta, -0.999))\n",
    "\n",
    "    \n",
    "    \n",
    "    epsilon_int = (epsilon_IA_TATT + epsilon_random)/(1 + np.conjugate(epsilon_IA_TATT) * epsilon_random)\n",
    "    epsilon_obs = (epsilon_int + g)/(1 + epsilon_int * np.conjugate(g))\n",
    "    epsilon_obs_no_noise = (epsilon_IA_TATT + g)/(1 + epsilon_IA_TATT * np.conjugate(g))\n",
    "   \n",
    "    fin_epsilon = Table([ra,dec,z,gamma1, gamma2, weight,delta,epsilon_obs.real,epsilon_obs.imag,epsilon_obs_no_noise.real,epsilon_obs_no_noise.imag,epsilon_IA_TATT.real,epsilon_IA_TATT.imag], names = ('ra','dec','z','g1', 'g2', 'weight','delta','e_obs_1','e_obs_2', 'e_obs_no_noise_1','e_obs_no_noise_2','e_IA_TATT_1', 'e_IA_TATT_2' ))\n",
    "    \n",
    "    #fin_epsilon.write('/global/cscratch1/sd/jharno/fin_epsilon_IA_{}.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('/global/cscratch1/sd/jharno/fin_epsilon_IA_{}_linear_bias.dat'.format(i), format='ascii', overwrite=True)    \n",
    "\n",
    "    #V0:\n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V0/fin_epsilon_IA_{}_0.06GpAM.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/random_pos/fin_epsilon_IA_{}_0.6GpAM_smooth0.1Mpc.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/random_pos/fin_epsilon_IA_{}_0.6GpAM_smooth0.1Mpc_m1.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/random_pos/fin_epsilon_IA_{}_0.6GpAM_smooth0.5Mpc_TT.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/random_pos/fin_epsilon_IA_{}_0.6GpAM_smooth0.5Mpc_TT_m1.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/NLA/galcat_skysim5000_noIA_tomo{}.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/NLA/galcat_skysim5000_IA_tomo{}_nofactor2.dat'.format(i), format='ascii', overwrite=True)    \n",
    "\n",
    "    \n",
    "    #V1-deltaNLA:\n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V1/fin_epsilon_IA_{}_0.06GpAM_linear_bias.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V1/fin_epsilon_IA_{}_0.06GpAM_linear_bias_smooth.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V1/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V1/fin_epsilon_IA_{}_0.06GpAM_linear_bias_smooth_1Mpc.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/lin_bias_pos/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth0.5Mpc.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/lin_bias_pos/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth0.5Mpc_m1.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/lin_bias_pos/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth0.5Mpc_bta2.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/lin_bias_pos/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth0.1Mpc_bta2_m1.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/lin_bias_pos/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth0.5Mpc_m1.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/lin_bias_pos/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth0.5Mpc_TT.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/lin_bias_pos/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth0.1Mpc_TT_m1.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/deltaNLA/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth1_smooth0.5Mpc_bta1.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/deltaNLA/galcat_skysim5000_noIA_bta1_tomo{}.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/deltaNLA/galcat_skysim5000_AIAp1_bta1_tomo{}_nofactor2.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/lin_bias_pos/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth1_smooth0.5Mpc_bta1_m1.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/lin_bias_pos/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth1_smooth0.5Mpc_bta1_111.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/lin_bias_pos/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth1_smooth0.5Mpc_bta2.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/lin_bias_pos/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth1_smooth0.5Mpc_bta2_TT.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/lin_bias_pos/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth1_smooth0.5Mpc_bta1_TT_V3.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/deltaTT/galcat_skysim5000_C2p1_bta1_tomo{}_nofactor2.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/deltaTT/galcat_skysim5000_C2p1_bta1_tomo{}_V2_nofactor2.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/lin_bias_pos/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth1_smooth0.5Mpc_bta1_TT_V3_m1.dat'.format(i), format='ascii', overwrite=True)    \n",
    "\n",
    "    \n",
    "    #Change bta:\n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/lin_bias_pos/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth0.5Mpc_bta1_to_bta0.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/SRD-Y1/lin_bias_pos/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth1_smooth0.5Mpc_bta1_to_bta0.dat'.format(i), format='ascii', overwrite=True)    \n",
    " \n",
    "    #V1-TT:\n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V1/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth_TT.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V1/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth_TT_V3.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    #fin_epsilon.write('../../GalCat/StageIV_nz/V1/fin_epsilon_IA_{}_0.6GpAM_linear_bias_smooth_TT_nodelta.dat'.format(i), format='ascii', overwrite=True)    \n",
    "    \n",
    "    #BCM:\n",
    "    fin_epsilon.write('../../../BCM/dmo/GalCat/fin_epsilon_NLA_{}_0.6GpAM_smooth05Mpc.dat'.format(i), format='ascii', overwrite=True)\n",
    "    \n",
    "    print('Done tomo',i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "1f3290fa-cd2f-4429-9c9d-6dcbd6d214d4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-rw-r--r-- 1 jharno lsst 2809571239 Oct 26 08:57 ../../GalCat/SRD-Y1/deltaNLA/galcat_skysim5000_AIAp1_bta1_tomo1_nofactor2.dat\n",
      "-rw-r--r-- 1 jharno lsst 2795994172 Oct 26 09:02 ../../GalCat/SRD-Y1/deltaNLA/galcat_skysim5000_AIAp1_bta1_tomo2_nofactor2.dat\n",
      "-rw-r--r-- 1 jharno lsst 2786861244 Oct 26 09:08 ../../GalCat/SRD-Y1/deltaNLA/galcat_skysim5000_AIAp1_bta1_tomo3_nofactor2.dat\n",
      "-rw-r--r-- 1 jharno lsst 2780572643 Oct 26 09:13 ../../GalCat/SRD-Y1/deltaNLA/galcat_skysim5000_AIAp1_bta1_tomo4_nofactor2.dat\n",
      "-rw-r--r-- 1 jharno lsst 2736517581 Oct 26 09:18 ../../GalCat/SRD-Y1/deltaNLA/galcat_skysim5000_AIAp1_bta1_tomo5_nofactor2.dat\n",
      "-rw-r--r-- 1 jharno lsst 2829202683 Oct 26 09:31 ../../GalCat/SRD-Y1/deltaTT/galcat_skysim5000_C2p1_bta1_tomo1_nofactor2.dat\n",
      "-rw-r--r-- 1 jharno lsst 2810905161 Oct 26 09:37 ../../GalCat/SRD-Y1/deltaTT/galcat_skysim5000_C2p1_bta1_tomo2_nofactor2.dat\n",
      "-rw-r--r-- 1 jharno lsst 2798760035 Oct 26 09:42 ../../GalCat/SRD-Y1/deltaTT/galcat_skysim5000_C2p1_bta1_tomo3_nofactor2.dat\n",
      "-rw-r--r-- 1 jharno lsst 2791499024 Oct 26 09:48 ../../GalCat/SRD-Y1/deltaTT/galcat_skysim5000_C2p1_bta1_tomo4_nofactor2.dat\n",
      "-rw-r--r-- 1 jharno lsst 2747558193 Oct 26 09:53 ../../GalCat/SRD-Y1/deltaTT/galcat_skysim5000_C2p1_bta1_tomo5_nofactor2.dat\n",
      "-rw-rw---- 1 jharno lsst 2831972738 Oct 26 08:26 ../../GalCat/SRD-Y1/NLA/galcat_skysim5000_AIAp1_tomo1_nofactor2.dat\n",
      "-rw-rw---- 1 jharno lsst 2816301924 Oct 26 08:32 ../../GalCat/SRD-Y1/NLA/galcat_skysim5000_AIAp1_tomo2_nofactor2.dat\n",
      "-rw-rw---- 1 jharno lsst 2804929447 Oct 26 08:38 ../../GalCat/SRD-Y1/NLA/galcat_skysim5000_AIAp1_tomo3_nofactor2.dat\n",
      "-rw-rw---- 1 jharno lsst 2797717709 Oct 26 08:43 ../../GalCat/SRD-Y1/NLA/galcat_skysim5000_AIAp1_tomo4_nofactor2.dat\n",
      "-rw-rw---- 1 jharno lsst 2752036620 Oct 26 08:49 ../../GalCat/SRD-Y1/NLA/galcat_skysim5000_AIAp1_tomo5_nofactor2.dat\n",
      "-rw-r--r-- 1 jharno lsst 2842359175 Oct 26 11:12 ../../GalCat/SRD-Y1/TT/galcat_skysim5000_c2p1_tomo1_nofactor2.dat\n",
      "-rw-r--r-- 1 jharno lsst 2822942056 Oct 26 11:18 ../../GalCat/SRD-Y1/TT/galcat_skysim5000_c2p1_tomo2_nofactor2.dat\n",
      "-rw-r--r-- 1 jharno lsst 2807301362 Oct 26 11:23 ../../GalCat/SRD-Y1/TT/galcat_skysim5000_c2p1_tomo3_nofactor2.dat\n",
      "-rw-r--r-- 1 jharno lsst 2800197717 Oct 26 11:28 ../../GalCat/SRD-Y1/TT/galcat_skysim5000_c2p1_tomo4_nofactor2.dat\n",
      "-rw-r--r-- 1 jharno lsst 2754512540 Oct 26 11:34 ../../GalCat/SRD-Y1/TT/galcat_skysim5000_c2p1_tomo5_nofactor2.dat\n"
     ]
    }
   ],
   "source": [
    "ls -l  ../../../BCM/dmo/GalCat/SRD-Y1/*/galcat_skysim5000*nofactor2.dat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f5e9a78c-a275-4a5d-91fa-c69d4f374a0f",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "got cats, downsampling...\n",
      "0.5999999710542959\n",
      "got shapes: (111383708, 6)\n",
      "extracted all quantities \n",
      "Done tomo 1\n",
      "got cats, downsampling...\n",
      "0.5999999710542959\n",
      "got shapes: (111390822, 6)\n",
      "extracted all quantities \n",
      "Done tomo 2\n",
      "got cats, downsampling...\n",
      "0.5999999710542959\n",
      "got shapes: (111378350, 6)\n",
      "extracted all quantities \n",
      "Done tomo 3\n",
      "got cats, downsampling...\n",
      "0.5999999710542959\n",
      "got shapes: (111382794, 6)\n",
      "extracted all quantities \n",
      "Done tomo 4\n",
      "got cats, downsampling...\n",
      "0.5999999710542959\n",
      "got shapes: (111394778, 6)\n",
      "extracted all quantities \n",
      "Done tomo 5\n"
     ]
    }
   ],
   "source": [
    "# HACC catalogues\n",
    "from astropy.io import fits\n",
    "\n",
    "nside=8192\n",
    "n_shells=101\n",
    "seed=0\n",
    "# dndz_flag='Euclid_dndz_fu08_bin1-5'\n",
    "dndz_flag='SRDv1_nz'\n",
    "kappa_dir=dndz_flag\n",
    "folder_path='/pscratch/sd/j/jmena/HOS_Sims/HACC150_3/shells_z'+str(n_shells)+'_subsampledensitymax1.5_nside'+str(nside)+'/'+kappa_dir\n",
    "filename='kappa_hacc_seed'+str(seed)+'_nside'+str(nside)+'_'\n",
    "nshells=51\n",
    "\n",
    "for i in range(1,6):\n",
    "    filename1 = folder_path+'/'+filename+'catpoisson_interppixels_b1_imap'+str(i-1)+'.fits'\n",
    "    hdul1=fits.open(filename1)\n",
    "    data1=hdul1[1].data\n",
    "    header1=hdul1[1].header\n",
    "\n",
    "            \n",
    "    data1_dict={}\n",
    "    for column in data1.columns:\n",
    "        if column.name in ['ra','dec', 'z', 'gamma1', 'gamma2']:\n",
    "            data1_dict[column.name]=data1[column.name]\n",
    "\n",
    "    print('got cats, downsampling...')\n",
    "      \n",
    "    # Let's downsample to 0.6 gal/arcmin**2\n",
    "    gal_den=0.6;\n",
    "    import random\n",
    "\n",
    "    selection1=sorted(random.sample(range(len(data1_dict['ra'])),int(gal_den*(4*np.pi*(180/np.pi)**2/8*60**2))))\n",
    "\n",
    "    data1_dict_downsampled={}\n",
    "    for column in data1_dict.keys():\n",
    "        data1_dict_downsampled[column]=data1_dict[column][selection1]\n",
    "\n",
    "    print(len(data1_dict_downsampled['ra'])/(4*np.pi*(180/np.pi)**2/8*60**2)) # this should print 0.6\n",
    "\n",
    "    \n",
    "    #epsilons = Table.read('../../../HACC-Y1/GalCat/SRD-Y1/HACC150_3/IA/epsilon_{}_0.6GpAM.asc'.format(i), format='ascii')\n",
    "    epsilons = np.loadtxt('../../../HACC-Y1/GalCat/SRD-Y1/HACC150_3/IA/epsilon_{}_0.6GpAM.asc'.format(i), skiprows=1)\n",
    "    #epsilons = Table.read('../../GalCat/SRD-Y1/lin_bias_pos/epsilon_{}_0.6GpAM_linear_bias_smooth1_bta2_smooth05Mpc.asc'.format(i), format='ascii')\n",
    "\n",
    "    print('got shapes:', np.shape(epsilons))\n",
    "    \n",
    "    ra=np.asarray(data1_dict_downsampled['ra'])\n",
    "    dec=np.asarray(data1_dict_downsampled['dec'])\n",
    "    weight=np.asarray(np.ones(ra.shape))\n",
    "    z=np.asarray(data1_dict_downsampled['z'])\n",
    "    gamma1 = np.asarray(data1_dict_downsampled['gamma1'])\n",
    "    gamma2 = np.asarray(data1_dict_downsampled['gamma2'])\n",
    "    \n",
    "    # e1_NLA = np.asarray(epsilons['e1_NLA'])\n",
    "    # e2_NLA = np.asarray(epsilons['e2_NLA'])\n",
    "    # e1_TT = np.asarray(epsilons['e1_TT'])\n",
    "    # e2_TT = np.asarray(epsilons['e2_TT'])\n",
    "    # e1_noise = np.asarray(epsilons['e1_noise'])\n",
    "    # e2_noise = np.asarray(epsilons['e2_noise'])\n",
    "    e1_NLA = np.asarray(epsilons[selection1,0])\n",
    "    e2_NLA = np.asarray(epsilons[selection1,1])\n",
    "    e1_TT = np.asarray(epsilons[selection1,2])\n",
    "    e2_TT = np.asarray(epsilons[selection1,3])\n",
    "    e1_noise = np.asarray(epsilons[selection1,4])\n",
    "    e2_noise = np.asarray(epsilons[selection1,5])\n",
    "\n",
    "    \n",
    "    print('extracted all quantities ')\n",
    "\n",
    "    \n",
    "    delta = np.asarray(np.ones(ra.shape))\n",
    "    epsilon_NLA = np.empty(e1_NLA.shape, dtype = complex)\n",
    "    epsilon_NLA.real = A_IA * e1_NLA  \n",
    "    epsilon_NLA.imag = A_IA * e2_NLA  \n",
    "    epsilon_TT = np.empty(e1_TT.shape, dtype = complex)\n",
    "    epsilon_TT.real = C_2 * e1_TT  \n",
    "    epsilon_TT.imag = C_2 * e2_TT\n",
    "    epsilon_random = np.empty(e1_TT.shape, dtype = complex)\n",
    "    epsilon_random.real = e1_noise\n",
    "    epsilon_random.imag = e2_noise\n",
    "    g = np.empty(gamma1.shape, dtype=complex)\n",
    "    g.real = gamma1\n",
    "    g.imag = gamma2\n",
    "    \n",
    "    \n",
    "    #calibrate all IA:\n",
    "    # Needed after finding the factor of 2.0 bug in the e_IA2\n",
    "    epsilon_NLA*=1.5\n",
    "    epsilon_TT*=1.5\n",
    "    \n",
    "    \n",
    "    #------------------\n",
    "    #calibrate delta:\n",
    "    #delta*=0.0\n",
    "    \n",
    "    #calibrate TT:\n",
    "    # for TT: \n",
    "    #epsilon_TT[z<0.5]/=2.5 \n",
    "    \n",
    "    # calibrate delta-TT\n",
    "    epsilon_TT[z<0.5]/=5.0 \n",
    "    epsilon_TT/=7 \n",
    "    #------------------\n",
    "    \n",
    "    #And finally, we calculate the observed epsilons using the summation as in equation 11 of the paper\n",
    "    # The V1 mocks already sample the galaxy by linearly tracing the matter field, so they already have the \n",
    "    # delta-weighted term. \n",
    "    # 1-Keep as is:\n",
    "    epsilon_IA_TATT = epsilon_NLA + epsilon_TT\n",
    "    # 2-Undo delta-weight:\n",
    "    #epsilon_IA_TATT = epsilon_NLA / (1 + b_TA * delta) + epsilon_TT\n",
    "    #epsilon_IA_TATT /= (1 + b_TA * delta) \n",
    "    # 3-Change from b_TA = 1 to b_TA = 0 or 2.0, uncomment the following: \n",
    "    #epsilon_IA_TATT *= (1 + b_TA * delta)/(1 + 1.0 * np.where(delta >=  -0.999, delta, -0.999))\n",
    "\n",
    "    \n",
    "    \n",
    "    epsilon_int = (epsilon_IA_TATT + epsilon_random)/(1 + np.conjugate(epsilon_IA_TATT) * epsilon_random)\n",
    "    epsilon_obs = (epsilon_int + g)/(1 + epsilon_int * np.conjugate(g))\n",
    "    epsilon_obs_no_noise = (epsilon_IA_TATT + g)/(1 + epsilon_IA_TATT * np.conjugate(g))\n",
    "   \n",
    "    fin_epsilon = Table([ra,dec,z,gamma1, gamma2, weight,delta,epsilon_obs.real,epsilon_obs.imag,epsilon_obs_no_noise.real,epsilon_obs_no_noise.imag,epsilon_IA_TATT.real,epsilon_IA_TATT.imag], names = ('ra','dec','z','g1', 'g2', 'weight','delta','e_obs_1','e_obs_2', 'e_obs_no_noise_1','e_obs_no_noise_2','e_IA_TATT_1', 'e_IA_TATT_2' ))\n",
    "    fin_epsilon.write('../../../HACC-Y1/GalCat/SRD-Y1/HACC150_3/NLA/galcat_skysim5000_AIAp1_poisson_tomo{}.dat'.format(i), format='ascii', overwrite=True)        \n",
    "    \n",
    "    print('Done tomo',i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9805f3ff-582a-42d2-a2e3-cbc99a03c5d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.io import fits\n",
    "nside=8192\n",
    "n_shells=101\n",
    "seed=0\n",
    "# dndz_flag='Euclid_dndz_fu08_bin1-5'\n",
    "dndz_flag='SRDv1_nz'\n",
    "kappa_dir=dndz_flag\n",
    "folder_path='/pscratch/sd/j/jmena/HOS_Sims/HACC150_3/shells_z'+str(n_shells)+'_subsampledensitymax1.5_nside'+str(nside)+'/'+kappa_dir\n",
    "filename='kappa_hacc_seed'+str(seed)+'_nside'+str(nside)+'_'\n",
    "nshells=51\n",
    "for i in range(1,2):\n",
    "    filename1 = folder_path+'/'+filename+'catpoisson_interppixels_b1_imap'+str(i-1)+'.fits'\n",
    "    hdul1=fits.open(filename1)\n",
    "    data1=hdul1[1].data\n",
    "    header1=hdul1[1].header\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2bd918af-b88c-43d0-9c08-98da3a09938d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffb889e8-2223-4716-b11e-40e109bc48fb",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "desc-python-dev",
   "language": "python",
   "name": "desc-python-dev"
  },
  "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}