epic_xml_files/epic_sixte_xml.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abandoned-vermont",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "import numpy as np\n",
    "from astropy.table import Table\n",
    "from astropy.io import fits\n",
    "from bs4 import BeautifulSoup"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "anticipated-whale",
   "metadata": {},
   "source": [
    "Definition of SIXTE's XML Instrument Configuration for XMM-Newton EPIC-PN camera\n",
    "================================================================================\n",
    "\n",
    "The current XMM files distributed by the SIXTE team include XIC (XML Instrument Configuration) files for a single PN CCD and for the whole camera. However, the last one just reproduces the size of the whole array. For a more realistic simulation, we need to define independent XIC files for each CCD.\n",
    "\n",
    "The geometry of the camera is as follows:\n",
    "\n",
    "![PN geometry](epicpn.png \"EPIC-PN camera and spatial coordinate systems\")\n",
    "\n",
    "* CC1 corresponds to CAMCOORD1, the camera reference system with origin in the center of the camera.\n",
    "* CC2 corresponds to CAMCOORD2, the camera reference system with origin in the optical axis of the corresponding telescope (XRT3).\n",
    "* X/Y axis defined for each CCD corresponds to PIXCOORD1, node-oriented CCD pixel coordinate system (RAWX/RAWY in the event list).\n",
    "* SC corresponds to SACCOORD, the spacecraft reference system. -Xsc axis corresponds to the viewing direction of the telescope. -Zsc axis, is North-aligned under a position angle of 0 deg, i.e. the declination axis.\n",
    "\n",
    "See the XMM-Newton's [Calibration Access and Data Handbook](https://xmm-tools.cosmos.esa.int/external/xmm_calibration/calib/documentation/CALHB/node7.html) for more details about these coordinate systems.\n",
    "\n",
    "According to [SIXTE's manual](https://www.sternwarte.uni-erlangen.de/research/sixte/data/simulator_manual_v1.3.14.pdf) (Sect 4.3): \"The instrument coordinate system (X,Y) is defined as a planar, Cartesian system located in the focal plane behind the instrument’s optics. Its units are meter, with the origin being the point of intersection between the optical axis and the focal plane. If an attitude file is used which has the keyword ALIGNMEN set to MOTION, the X-axis of the instrument coordinate system points along the direction of motion of the pointing vector. In any other case, the X-axis points towards the celestial north pole. The Y -axis then points towards the west. If additionally a Roll Angle is defined in the attitude file, the X- and Y -axes as described before are rotated around the optical axis by this angle.\"\n",
    "\n",
    "So, if I understand all this correctly, the SIXTE's instrument coordinate system (SXT) would be equivalent to the CAMCOORD2 reference system. The position of each CCD is given in the XMM calibration files (EPN_LINCOORD_0009.CCF is the most recent version as of 2021/09/24) using the center of the CCD as the reference pixel in PIXCOORD1 (RAWX0, RAWY0) in all CCD and giving its position in the CAMCOORD1 system (X0, Y0). Note that the orientation of PIXCOORD1 is the same as CAMCOORD1/2 for CCDs 7-12, but it is rotated 180 deg for CCDs 1-6. I need to convert the CC1 positions of the reference pixels into CC2 coordinates, and use that in the XIC files. The conversion factors for CC1->CC2 are given in the header of the calibration file and, as shown in the diagram above, is just a coordinates shift.\n",
    "\n",
    "XIC files use the following parameters to define the geometry and position of a detector:\n",
    "* xwidth, ywidth: Number of pixels in X_SXT and Y_SXT direction.\n",
    "* xrpix, yrpix: The reference pixel.\n",
    "* xrval, yrval: The location of the reference pixel in the focal plane (in m).\n",
    "* xdelt, ydelt: The dimensions of the pixels (in m)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "functioning-blogger",
   "metadata": {},
   "outputs": [],
   "source": [
    "filter = \"thin\"\n",
    "add_diehlbkg = False\n",
    "add_flatbkg = False\n",
    "add_cfbkg = True\n",
    "bkg_level = \"low\" # high, med, low\n",
    "\n",
    "\n",
    "# We use the fullframe XIC file as a reference for defining the new files for each CCD\n",
    "sixte_instruments_epicpn_path = Path(\"/\", \"storage\", \"sixte_instrument_files\", \"xmm\", \"epicpn\")\n",
    "epicpn_fullframe_xic = sixte_instruments_epicpn_path.joinpath(f\"fullframe_{filter}filter.xml\")\n",
    "\n",
    "# We select the most recent xmm calibration\n",
    "xmmccf_path = Path(\"/\", \"storage\", \"xmmccf\")\n",
    "epicpn_lincoord_ccf = sorted(xmmccf_path.glob(\"EPN_LINCOORD_*.CCF\"))[-1]\n",
    "lincoord = Table.read(epicpn_lincoord_ccf, hdu=1)\n",
    "    \n",
    "for ccd in lincoord:\n",
    "    with epicpn_fullframe_xic.open() as fp:\n",
    "        xic = BeautifulSoup(fp, \"xml\")\n",
    "\n",
    "    xic.instrument[\"instrume\"] = \"EPN\"\n",
    "    xic.fov[\"diameter\"] = 0.5\n",
    "    \n",
    "    xic.dimensions[\"xwidth\"] = 64\n",
    "    xic.dimensions[\"ywidth\"] = 200\n",
    "    \n",
    "    xic.wcs[\"xrpix\"] = ccd[\"RAWX0\"] + 0.5  # The pixel coordinates in the CCF correspond to the center of the pixel\n",
    "    xic.wcs[\"yrpix\"] = ccd[\"RAWY0\"] + 0.5\n",
    "\n",
    "    xic.wcs[\"xrval\"] = (ccd[\"X0\"] - lincoord.meta[\"CC12_TX\"]) * 1e-3  # Shift to CC2 system and convert from mm to meters\n",
    "    xic.wcs[\"yrval\"] = (ccd[\"Y0\"] - lincoord.meta[\"CC12_TY\"]) * 1e-3\n",
    "\n",
    "    xic.wcs[\"xdelt\"] = lincoord.meta[\"MM_PX_X\"] * 1e-3  # Convert from mm to meters\n",
    "    xic.wcs[\"ydelt\"] = lincoord.meta[\"MM_PX_Y\"] * 1e-3\n",
    "    \n",
    "    # TODO: check if a change in the readout direction is needed\n",
    "    if ccd['CCD_ID'] <= 6:\n",
    "        xic.wcs[\"rota\"] = 180.\n",
    "    else:\n",
    "        xic.wcs[\"rota\"] = 0.\n",
    "        \n",
    "    xic.rmf[\"filename\"] = f\"pn-{filter}-10.rmf\"\n",
    "    xic.arf[\"filename\"] = f\"pn-{filter}-10.arf\"\n",
    "    \n",
    "    xic_path = Path(\"xml\", f\"fullframe_ccd{ccd['CCD_ID']:02}_{filter}filter.xml\")\n",
    "    \n",
    "    if add_diehlbkg:\n",
    "        bkgtag = xic.new_tag(\"phabackground\", filename=f\"pntffu_{bkg_level}_background_spectrum.fits\")\n",
    "        xic.arf.insert_after(bkgtag)\n",
    "\n",
    "        xic_path = xic_path.parent.joinpath(\n",
    "            f\"{xic_path.stem}_{bkg_level}diehlbkg{xic_path.suffix}\"\n",
    "        )\n",
    "    \n",
    "    elif add_flatbkg:\n",
    "        bkgtag = xic.new_tag(\"phabackground\", filename=f\"pnflat_{bkg_level}_background_spectrum.fits\")\n",
    "        xic.arf.insert_after(bkgtag)\n",
    "\n",
    "        xic_path = xic_path.parent.joinpath(\n",
    "            f\"{xic_path.stem}_{bkg_level}flatbkg{xic_path.suffix}\"\n",
    "        )\n",
    "\n",
    "    elif add_cfbkg:\n",
    "        bkgtag = xic.new_tag(\"phabackground\", filename=f\"pnclosedfilter_{bkg_level}_background_spectrum.fits\")\n",
    "        xic.arf.insert_after(bkgtag)\n",
    "\n",
    "        xic_path = xic_path.parent.joinpath(\n",
    "            f\"{xic_path.stem}_{bkg_level}cfbkg{xic_path.suffix}\"\n",
    "        )\n",
    "\n",
    "    with xic_path.open(\"w\") as fp:\n",
    "        fp.write(xic.prettify())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "hidden-abortion",
   "metadata": {},
   "source": [
    "EPIC-MOS\n",
    "----------\n",
    "\n",
    "Now I need to define the XIC files for the MOS cameras. The geometry of the detectors is as follows:\n",
    "\n",
    "![MOS geometry](epicmos.png \"EPIC-MOS1/2 cameras and spatial coordinate systems\")\n",
    "\n",
    "For MOS1, the SIXTE reference system is oriented such as that X_SXT coincides with the Y axis of CC2, and Y_SXT is -X_CC2, i.e. SXT and CC2 are rotated 90 deg with respect to each other. Then I have to rotate each detector so the RAWX/Y matches the correct orientation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "romantic-bumper",
   "metadata": {},
   "outputs": [],
   "source": [
    "filter = \"thin\"\n",
    "add_flatbkg = False\n",
    "add_cfbkg = False\n",
    "bkg_level = \"low\"\n",
    "\n",
    "# We use the fullframe XIC file as a reference for defining the new files for each CCD\n",
    "sixte_instruments_epicmos_path = Path(\"/\", \"storage\", \"sixte_instrument_files\", \"xmm\", \"epicmos\")\n",
    "epicmos1_fullframe_xic = sixte_instruments_epicmos_path.joinpath(f\"mos1_center_{filter}filter.xml\")\n",
    "\n",
    "# We select the most recent xmm calibration\n",
    "epicmos1_lincoord_ccf = sorted(xmmccf_path.glob(\"EMOS1_LINCOORD_*.CCF\"))[-1]\n",
    "lincoord = Table.read(epicmos1_lincoord_ccf, hdu=1)\n",
    "\n",
    "for ccd in lincoord:\n",
    "    if ccd[\"NODE_ID\"] == 1:\n",
    "        continue\n",
    "    \n",
    "    with epicmos1_fullframe_xic.open() as fp:\n",
    "        xic = BeautifulSoup(fp, \"xml\")\n",
    "        \n",
    "    xic.instrument[\"instrume\"] = \"EMOS1\"\n",
    "    xic.instrument[\"telescop\"] = \"XMM\"\n",
    "    xic.fov[\"diameter\"] = 0.5\n",
    "    \n",
    "    vigtag = xic.new_tag(\"vignetting\", filename=\"xmm_mos1_vignet.fits\")\n",
    "    xic.psf.insert_after(vigtag)\n",
    "    \n",
    "    xic.dimensions[\"xwidth\"] = 600\n",
    "    xic.dimensions[\"ywidth\"] = 600\n",
    "    \n",
    "    xic.wcs[\"xrpix\"] = ccd[\"RAWX0\"] + 0.5  # The pixel coordinates in the CCF correspond to the center of the pixel\n",
    "    xic.wcs[\"yrpix\"] = ccd[\"RAWY0\"] + 0.5\n",
    "\n",
    "    xic.wcs[\"xrval\"] = (ccd[\"Y0\"] - lincoord.meta[\"CC12_TY\"]) * 1e-3  # Shift to CC2 system and convert from mm to meters\n",
    "    xic.wcs[\"yrval\"] = -(ccd[\"X0\"] - lincoord.meta[\"CC12_TX\"]) * 1e-3\n",
    "\n",
    "    xic.wcs[\"xdelt\"] = lincoord.meta[\"MM_PX_X\"] * 1e-3  # Convert from mm to meters\n",
    "    xic.wcs[\"ydelt\"] = lincoord.meta[\"MM_PX_Y\"] * 1e-3\n",
    "    \n",
    "    # TODO: check if a change in the readout direction is needed\n",
    "    if ccd['CCD_ID'] == 1:\n",
    "        xic.wcs[\"rota\"] = 90.\n",
    "    \n",
    "    elif ccd['CCD_ID'] >= 2 and ccd['CCD_ID'] <=4:\n",
    "        xic.wcs[\"rota\"] = 0.\n",
    "    \n",
    "    elif ccd['CCD_ID'] >= 5:\n",
    "        xic.wcs[\"rota\"] = 180.\n",
    "        \n",
    "    xic.rmf[\"filename\"] = f\"mos1-{filter}-10.rmf\"\n",
    "    xic.arf[\"filename\"] = f\"mos1-{filter}-10.arf\"\n",
    "    \n",
    "    xic_path = Path(\"xml\", f\"mos1_fullframe_ccd{ccd['CCD_ID']:02}_{filter}filter.xml\")\n",
    "    \n",
    "    if add_flatbkg:\n",
    "        bkgtag = xic.new_tag(\"phabackground\", filename=f\"mos1flat_{bkg_level}_background_spectrum.fits\")\n",
    "        xic.arf.insert_after(bkgtag)\n",
    "\n",
    "        xic_path = xic_path.parent.joinpath(\n",
    "            f\"{xic_path.stem}_{bkg_level}flatbkg{xic_path.suffix}\"\n",
    "        )\n",
    "        \n",
    "    elif add_cfbkg:\n",
    "        bkgtag = xic.new_tag(\"phabackground\", filename=f\"mos1closedfilter_{bkg_level}_background_spectrum.fits\")\n",
    "        xic.arf.insert_after(bkgtag)\n",
    "\n",
    "        xic_path = xic_path.parent.joinpath(\n",
    "            f\"{xic_path.stem}_{bkg_level}cfbkg{xic_path.suffix}\"\n",
    "        )\n",
    "\n",
    "    with xic_path.open(\"w\") as fp:\n",
    "        fp.write(xic.prettify())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "loved-organ",
   "metadata": {},
   "source": [
    "For MOS2 the SIXTE reference system is oriented such as that X_SXT coincides with the -X axis of CC2, and Y_SXT is -Y_CC2, i.e. SXT and CC2 are rotated 180 deg with respect to each other. Then I have to rotate each detector so the RAWX/Y matches the correct orientation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "capable-point",
   "metadata": {},
   "outputs": [],
   "source": [
    "filter = \"thin\"\n",
    "add_flatbkg = False\n",
    "add_cfbkg = False\n",
    "bkg_level = \"low\"\n",
    "\n",
    "# We use the fullframe XIC file as a reference for defining the new files for each CCD\n",
    "sixte_instruments_epicmos_path = Path(\"/\", \"storage\", \"sixte_instrument_files\", \"xmm\", \"epicmos\")\n",
    "epicmos2_fullframe_xic = sixte_instruments_epicmos_path.joinpath(f\"mos2_center_{filter}filter.xml\")\n",
    "\n",
    "# We select the most recent xmm calibration\n",
    "epicmos2_lincoord_ccf = sorted(xmmccf_path.glob(\"EMOS2_LINCOORD_*.CCF\"))[-1]\n",
    "lincoord = Table.read(epicmos2_lincoord_ccf, hdu=1)\n",
    "\n",
    "for ccd in lincoord:\n",
    "    if ccd[\"NODE_ID\"] == 1:\n",
    "        continue\n",
    "    \n",
    "    with epicmos1_fullframe_xic.open() as fp:\n",
    "        xic = BeautifulSoup(fp, \"xml\")\n",
    "        \n",
    "    xic.instrument[\"instrume\"] = \"EMOS2\"\n",
    "    xic.instrument[\"telescop\"] = \"XMM\"\n",
    "    xic.fov[\"diameter\"] = 0.5\n",
    "    \n",
    "    vigtag = xic.new_tag(\"vignetting\", filename=\"xmm_mos2_vignet.fits\")\n",
    "    xic.psf.insert_after(vigtag)\n",
    "    \n",
    "    xic.dimensions[\"xwidth\"] = 600\n",
    "    xic.dimensions[\"ywidth\"] = 600\n",
    "    \n",
    "    xic.wcs[\"xrpix\"] = ccd[\"RAWX0\"] + 0.5  # The pixel coordinates in the CCF correspond to the center of the pixel\n",
    "    xic.wcs[\"yrpix\"] = ccd[\"RAWY0\"] + 0.5\n",
    "\n",
    "    xic.wcs[\"xrval\"] = -(ccd[\"X0\"] - lincoord.meta[\"CC12_TX\"]) * 1e-3  # Shift to CC2 system and convert from mm to meters\n",
    "    xic.wcs[\"yrval\"] = -(ccd[\"Y0\"] - lincoord.meta[\"CC12_TY\"]) * 1e-3\n",
    "\n",
    "    xic.wcs[\"xdelt\"] = lincoord.meta[\"MM_PX_X\"] * 1e-3  # Convert from mm to meters\n",
    "    xic.wcs[\"ydelt\"] = lincoord.meta[\"MM_PX_Y\"] * 1e-3\n",
    "    \n",
    "    # TODO: check if a change in the readout direction is needed\n",
    "    if ccd['CCD_ID'] == 1:\n",
    "        xic.wcs[\"rota\"] = 180.\n",
    "    \n",
    "    elif ccd['CCD_ID'] >= 2 and ccd['CCD_ID'] <=4:\n",
    "        xic.wcs[\"rota\"] = -90.\n",
    "    \n",
    "    elif ccd['CCD_ID'] >= 5:\n",
    "        xic.wcs[\"rota\"] = 90.\n",
    "        \n",
    "    xic.rmf[\"filename\"] = f\"mos2-{filter}-10.rmf\"\n",
    "    xic.arf[\"filename\"] = f\"mos2-{filter}-10.arf\"\n",
    "    \n",
    "    xic_path = Path(\"xml\", f\"mos2_fullframe_ccd{ccd['CCD_ID']:02}_{filter}filter.xml\")\n",
    "    \n",
    "    if add_flatbkg:\n",
    "        bkg = \"low\"\n",
    "        bkgtag = xic.new_tag(\"phabackground\", filename=f\"mos1flat_{bkg}_background_spectrum.fits\")\n",
    "        xic.arf.insert_after(bkgtag)\n",
    "\n",
    "        xic_path = xic_path.parent.joinpath(\n",
    "            f\"{xic_path.stem}_{bkg_level}flatbkg{xic_path.suffix}\"\n",
    "        )\n",
    "            \n",
    "    elif add_cfbkg:\n",
    "        bkg = \"low\"\n",
    "        bkgtag = xic.new_tag(\"phabackground\", filename=f\"mos2closedfilter_{bkg}_background_spectrum.fits\")\n",
    "        xic.arf.insert_after(bkgtag)\n",
    "\n",
    "        xic_path = xic_path.parent.joinpath(\n",
    "            f\"{xic_path.stem}_{bkg_level}cfbkg{xic_path.suffix}\"\n",
    "        )\n",
    "\n",
    "    with xic_path.open(\"w\") as fp:\n",
    "        fp.write(xic.prettify())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "strong-valentine",
   "metadata": {},
   "source": [
    "MOS vignetting\n",
    "--------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "instrumental-flush",
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://github.com/SamSweere/xmm-epicpn-simulator/blob/master/other/vignet_creation/create_xmm_vinget_sixte.ipynb\n",
    "detector = \"MOS1\"\n",
    "xrt_xareaef_ccf = sorted(xmmccf_path.glob(f\"XRT{detector[-1]}_XAREAEF_*.CCF\"))[-1]\n",
    "print(xrt_xareaef_ccf)\n",
    "\n",
    "with fits.open(xrt_xareaef_ccf) as hdu:\n",
    "    # load the energy bins\n",
    "    energy = hdu['VIGNETTING'].data['ENERGY']\n",
    "\n",
    "    # Transform them to the sixte format\n",
    "    energy_lo = np.array(energy[:-1]/1000.0).reshape(1, len(energy[:-1])) # E to keV \n",
    "    energy_hi = np.array(energy[1:]/1000.0).reshape(1, len(energy[1:])) # # E to keV\n",
    "    \n",
    "    # Create the theta step_size bins\n",
    "    theta_stepsize = hdu['VIGNETTING'].header['D_THETA']\n",
    "    theta_bins = hdu['VIGNETTING'].data['VIGNETTING_FACTOR'].shape[1]\n",
    "    theta = np.array([i*theta_stepsize for i in range(theta_bins)], ndmin=2)\n",
    "\n",
    "    # We do not have a phi, thus an array of zero\n",
    "    phi = np.zeros((1,1))\n",
    "    \n",
    "    # Create the vignet array\n",
    "    vignet = np.zeros((1,theta.shape[1]-1,energy_lo.shape[1])) # , dtype=[('VIGNET', '>f4')]\n",
    "\n",
    "    ccf_vignet_data = hdu['VIGNETTING'].data['VIGNETTING_FACTOR']\n",
    "\n",
    "    for i in range(ccf_vignet_data.shape[0]-1):\n",
    "        # i indexes the theta variable\n",
    "        # print(ccf_vignet_data[i][0])\n",
    "        for j in range(ccf_vignet_data.shape[1]-1):\n",
    "            # j intexes the energy bin variable\n",
    "            # Convert to keV\n",
    "            vignet[0][j][i] = ccf_vignet_data[i][j]\n",
    "\n",
    "            \n",
    "# Save vignetting into a fits file            \n",
    "# The format is based on: https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_021/cal_gen_92_021.html\n",
    "energy_lo_col = fits.Column('ENERG_LO', format=f'{energy_lo.shape[1]}E', unit='keV', dim=f'{energy_lo.shape[1]}', array=energy_lo) \n",
    "energy_hi_col = fits.Column('ENERG_HI', format=f'{energy_hi.shape[1]}E', unit='keV', dim=f'{energy_hi.shape[1]}', array=energy_hi)\n",
    "theta_col = fits.Column('THETA', format=f'{theta.shape[1]}E', unit='degree', dim=f'{theta.shape[1]}', array=theta)\n",
    "phi_col = fits.Column('PHI', format=f'{phi.shape[1]}E', unit='degree', dim=f'{phi.shape[1]}', array=phi)\n",
    "vignet_col = fits.Column('VIGNET', format=f'{vignet.shape[0]*vignet.shape[1]*vignet.shape[2]}E', dim=f'{vignet.shape[::-1]}', array=vignet)\n",
    "\n",
    "columns = fits.ColDefs([energy_lo_col, energy_hi_col, theta_col, phi_col, vignet_col])\n",
    "\n",
    "# Include the required headers\n",
    "vignet_hdr = fits.Header()\n",
    "vignet_hdr['HDUCLASS'] = 'OGIP'\n",
    "vignet_hdr['HDUCLAS1'] = 'RESPONSE'\n",
    "vignet_hdr['HDUVERS1'] = '1.0.0'\n",
    "vignet_hdr['HDUCLAS2'] = 'VIGNET'\n",
    "vignet_hdr['HDUVERS2'] = '1.1.0'\n",
    "vignet_hdr['VERSION'] = '20171016'\n",
    "vignet_hdr['MISSION'] = 'XMM'\n",
    "vignet_hdr['TELESCOP'] = 'XMM'\n",
    "vignet_hdr['DETNAM'] = 'XMM'\n",
    "vignet_hdr['INSTRUME'] = f'E{detector}'\n",
    "\n",
    "vignet_hdr['TUNIT1'] = 'keV'\n",
    "vignet_hdr['TUNIT2'] = 'keV'\n",
    "vignet_hdr['TUNIT3'] = 'degree'\n",
    "vignet_hdr['TUNIT4'] = 'degree'\n",
    "\n",
    "vignet_hdr.add_history(f'Produced by Angel Ruiz (NOA) according to data from the XMM-{detector} calibration file {xrt_xareaef_ccf.name}')\n",
    "\n",
    "vignet_hdu = fits.BinTableHDU.from_columns(columns, header=vignet_hdr)\n",
    "vignet_hdu.name = 'VIGNET'\n",
    "\n",
    "# Create the final hdul\n",
    "primary_hdu = fits.PrimaryHDU()\n",
    "hdul = fits.HDUList([primary_hdu, vignet_hdu])\n",
    "\n",
    "hdul_path = Path(\"xml\", f\"xmm_{detector.lower()}_vignet.fits\")\n",
    "hdul.writeto(str(hdul_path), overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "parliamentary-carrier",
   "metadata": {},
   "source": [
    "Particle background models\n",
    "============================\n",
    "\n",
    "We created two different particle background model to be used with SIXTE for XMM-Newton simulations.\n",
    "\n",
    "The first one is a simple flat background spectrum (as is used in Athena):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "colonial-accommodation",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We just modify an existing fits file to makes things easier\n",
    "bkg_level = \"low\"\n",
    "spec_fits_file = Path(\"xml\", f\"pnflat_{bkg_level}_background_spectrum.fits\")\n",
    "\n",
    "with fits.open(spec_fits_file, \"update\") as hdu:\n",
    "    hdu[1].data[\"CHANNEL\"] = np.arange(len(hdu[1].data[\"RATE\"]))\n",
    "    hdu[1].data[\"RATE\"] = 1. * np.ones(len(hdu[1].data[\"RATE\"]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "banned-president",
   "metadata": {},
   "outputs": [],
   "source": [
    "bkg_level = \"low\"\n",
    "spec_fits_file = Path(\"xml\", f\"mos1flat_{bkg_level}_background_spectrum.fits\")\n",
    "\n",
    "spec = Table()\n",
    "spec[\"CHANNEL\"] = np.arange(2400, dtype=np.int32)\n",
    "spec[\"CHANNEL\"].unit = \"\"\n",
    "spec[\"RATE\"] = 0.1 * np.ones(2400)\n",
    "spec[\"RATE\"].unit = \"cts/s/bin/m/m\"\n",
    "\n",
    "spec_hdu = fits.BinTableHDU(spec)\n",
    "primary_hdu = fits.PrimaryHDU()\n",
    "\n",
    "hdul = fits.HDUList([primary_hdu, spec_hdu])\n",
    "hdul[1].name = \"SPECTRUM\"\n",
    "hdul.writeto(spec_fits_file.as_posix(), overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "atomic-breakdown",
   "metadata": {},
   "source": [
    "The second one is a more realistic spectrum, based on actual data from closed-filter observations.\n",
    "\n",
    "* I am using the merged event file produced by merging various closed filter EPIC-PN observations. The corresponding Merged Event File can be downloaded from here https://www.cosmos.esa.int/web/xmm-newton/filter-closed\n",
    "\n",
    "* This file lists for every event the PI (Pulse Invariant channel) in units of eV. I have used this quantitity to identify the correspond EPIC PN channel using the EBOUNDS extension of an rmf (epn_e1_ff20_sdY0_v20.0.rmf - but it does not matter which specific file is used) which gives the energy limits (EMIN EMAX) of each of the 4096 channels. The PHA (Pulse Height Analyser), as far as I understand, is not the correct quantity to use for identifying channels.\n",
    "\n",
    "* The RATE value (in counts s-1 m-2) that needs to be included in a SIXTE particle background is calculated building the histogram of the PI in the events file. I defined the histogram bins using a list of the EMIN, EMAX of the channels taken from the RMF. I have also divided by the exposure time of the merged event file (755.7e3 s) and the area in m2 of the EPIC-PN (36e-4 m2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "developmental-lithuania",
   "metadata": {},
   "outputs": [],
   "source": [
    "spec_file = Path(\"..\", \"xspec\", \"xmm_particle_backspec_new.txt\")\n",
    "spec0 = Table.read(spec_file, format=\"ascii\")\n",
    "\n",
    "# We just modify an existing fits file to makes things easier\n",
    "bkg_level = \"low\"\n",
    "spec_fits_file = Path(\"xml\", f\"pnclosedfilter_{bkg_level}_background_spectrum.fits\")\n",
    "\n",
    "\n",
    "spec = Table()\n",
    "spec[\"CHANNEL\"] = np.arange(4096, dtype=np.int32)\n",
    "spec[\"CHANNEL\"].unit = \"\"\n",
    "spec[\"RATE\"] = 2 * spec0[\"col2\"]\n",
    "spec[\"RATE\"].unit = \"cts/s/bin/m/m\"\n",
    "\n",
    "spec_hdu = fits.BinTableHDU(spec)\n",
    "primary_hdu = fits.PrimaryHDU()\n",
    "\n",
    "hdul = fits.HDUList([primary_hdu, spec_hdu])\n",
    "hdul[1].name = \"SPECTRUM\"\n",
    "hdul.writeto(spec_fits_file.as_posix(), overwrite=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pyxsel (source detection)",
   "language": "python",
   "name": "pyxsel"
  },
  "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
}