photdm_impl.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e47961e2-2622-457e-845a-f55da75f2a61",
   "metadata": {},
   "source": [
    "## PhotDM Reference Implementation\n",
    "\n",
    "This notebook should how model annotations can be used to plot SEDs from a VOTable.\n",
    "\n",
    "The annotation bind the data cells with the description of their Photometric calibration.\n",
    "- Some brighness are given as magnitudes and some others as fluxes.\n",
    "- The photometric calibrations are given as instances of PhotDM.\n",
    "- The challenge of that code is to build SEDs (wavelength vs fluxes) just by using these annotations. \n",
    "  - There no use of any data not within the VOTable\n",
    "  - All data access are done through the model view of each data row.\n",
    "\n",
    "#### Input Data\n",
    "\n",
    "The input VOTable (NWAYAWGXMM.xml) provided by HEASARC is a subset of the XMMSL2 catalogue.\n",
    "Let's consider the quantities of interest :\n",
    "- Each XMMSL2 source comes with a name, a position and fluxes in EB6,7,8 bands\n",
    "- For each XMMSL2 sources, we get the following counterpart magnitudes:\n",
    "  - Wise W1,2,3,4\n",
    "  - 2MASS Ks,H,J\n",
    "  - GAIA G\n",
    "The others columns are ignored by this proof of concept.\n",
    "\n",
    "#### Annotation Model\n",
    "\n",
    "The goal of this exercice is to validate PhotDM.\n",
    "To make this example working, we need a model that encompasses the photometric calibrations with the brightness data.\n",
    "For this purpose, we have built a mock model named `sed`.\n",
    "- A sed is\n",
    "  - A name\n",
    "  - a position\n",
    "  - a collection a photometric points which are\n",
    "    - A magnitude or a flux\n",
    "    - A pointer to the PhotCal instance\n",
    "    \n",
    "#### The code\n",
    "The code iterates over the data table and build an `sed` instance for each row. \n",
    "The magnitudes are converted in fluxes and one SED is then plotted.\n",
    "\n",
    "#### Watchout: X-Ray Photcal\n",
    "In this data set, we assimilate the X-Ray energy bands to photometric calibrations without zero point and with a square filter profile.\n",
    "This approach is questionnable, but it is workable in this example.\n",
    "The reason is that we use optical/IR calibrations to transform magnitudes into fluxes that are compliant with X-Ray fluxes.\n",
    "X_Ray calibrations are just used to compute the wavelengths.\n",
    "\n",
    "LM 28/02/2022"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a66ddfff-15d0-40ad-a5cf-085376088105",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys\n",
    "import numpy as np\n",
    "import logging\n",
    "import matplotlib\n",
    "import mplcursors\n",
    "import matplotlib.pyplot as plt \n",
    "from astropy.io.votable import parse\n",
    "from astropy import units as u\n",
    "base_dir = os.getcwd()\n",
    "print(base_dir)\n",
    "sys.path.append(os.path.join(base_dir, \"python\") )\n",
    "data_dir = os.path.join(base_dir, \"python\", \"examples\", \"data\") \n",
    "                             \n",
    "from client.photdm.photcal import PhotCal\n",
    "from client.xml_interpreter.model_viewer import ModelViewer\n",
    "\n",
    "matplotlib.font_manager: logging.getLogger('matplotlib.font_manager').setLevel(logging.WARNING)\n",
    "matplotlib.font_manager: logging.getLogger('matplotlib.ticker').setLevel(logging.WARNING)\n",
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23a2d0bf-9721-4a04-9f75-55db38278e0f",
   "metadata": {},
   "source": [
    "Build the model viewer that will be used to explore the data set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d724316-3049-4e40-b020-973b1ab7d830",
   "metadata": {},
   "outputs": [],
   "source": [
    "votable = parse(os.path.join(data_dir, \"NWAYAWGXMM.xml\"))\n",
    "for resource in votable.resources:\n",
    "    mviewer = ModelViewer(resource, votable_path=os.path.join(data_dir, \"NWAYAWGXMM.xml\"))\n",
    "    break;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61122998-d054-4cca-acc5-becdb038084f",
   "metadata": {},
   "source": [
    "### Main loop (limited to 1 rows)\n",
    "For each row,\n",
    "- get the model view\n",
    "- for each point in the row\n",
    "  - build the PhotCal instances\n",
    "  - convert the magnitudes into fluxes\n",
    "  - plot the SED"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04ea705a-7c72-4725-9d27-bb3d507c8bb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "mviewer.connect_table(None)\n",
    "times = []\n",
    "ras = []\n",
    "decs = []\n",
    "cpt=0\n",
    "while True:\n",
    "    row = mviewer.get_next_row() \n",
    "    if row is None:\n",
    "        break  \n",
    "    row_model_view = mviewer.get_model_view(resolve_ref=False)\n",
    "    phot_points = mviewer.get_model_component_by_type(\"sed:PhotPoint\")\n",
    "    wl = []\n",
    "    flx = []\n",
    "    tlt = []\n",
    "    name = None\n",
    "    for ele in row_model_view.xpath('.//ATTRIBUTE[@dmrole=\"sed:Source.name\"]'):\n",
    "        name  = ele.get(\"value\")\n",
    "        break\n",
    "    for phot_point in phot_points:\n",
    "        #XmlUtils.pretty_print(phot_point)\n",
    "        flux = None\n",
    "        unit = None\n",
    "        unit_org = None\n",
    "        phot_cal = None\n",
    "        flux_erg = None\n",
    "        for ele in phot_point.xpath('.//ATTRIBUTE[@dmrole=\"sed:PhotPoint.flux\"]'):\n",
    "            try:\n",
    "                flux = float(ele.get(\"value\"))\n",
    "                unit_org = ele.get(\"unit_org\")\n",
    "                unit = ele.get(\"unit\")\n",
    "            except:\n",
    "                pass\n",
    "            break\n",
    "        for ele in phot_point.xpath('.//INSTANCE[@dmrole=\"sed:PhotPoint.photcal\"]'):\n",
    "            phot_cal = PhotCal(ele)\n",
    "            spec_loc = phot_cal.photometryFilter.spectralLocation\n",
    "            spec_loc_value = None\n",
    "            # convert XMM energy into wavelength\n",
    "            if spec_loc.unitexpression == \"keV\":\n",
    "                spec_loc_value = (spec_loc.value * u.keV).to(u.angstrom, equivalencies=u.spectral()).value\n",
    "            else: \n",
    "                spec_loc_value = spec_loc.value\n",
    "            break\n",
    "        if flux is not None:\n",
    "            if unit_org == \"mag\":\n",
    "                zeroPoint = phot_cal.zeroPoint.flux.value\n",
    "                refMag = phot_cal.zeroPoint.referenceMagnitudeValue\n",
    "                flux_erg = 1e-23 * zeroPoint * np.power( 10, - (flux - refMag)/2.5) * (299792458 / (spec_loc_value*1E-10))\n",
    "            else:\n",
    "                flux_erg = flux\n",
    "            flx.append(flux_erg)\n",
    "            wl.append(spec_loc_value)\n",
    "            tlt.append(phot_cal.identifier)\n",
    "            print(f\"{phot_cal.identifier} {spec_loc.value} {spec_loc.unitexpression} {flux_erg}  {unit}\")\n",
    "\n",
    "    plt.clf()\n",
    "    plt.scatter(wl, flx)\n",
    "    plt.title(\"Source \" + name)\n",
    "    plt.xscale(\"log\")\n",
    "    # Lim must be force to work on Linux Py3.6\n",
    "    plt.ylim([min(flx)/10, max(flx)*10])\n",
    "    plt.yscale(\"log\")\n",
    "    plt.xlabel(\"Wavelength (Angstrom)\")\n",
    "    plt.ylabel(\"flux (erg/cm2/sec)\")\n",
    "    mplcursors.cursor(hover=True).connect(\"add\", lambda sel: sel.annotation.set_text(tlt[sel.index]))\n",
    "    plt.subplots_adjust(bottom=.15, left=.15) \n",
    "    plt.show()\n",
    "    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2550b13b-af6c-45c6-8570-099ecf9ff837",
   "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}