crab.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "reference_instrument=\"spi\"\n",
    "subcases_pattern=\"_0605\"\n",
    "osa_version=\"OSA10.2\"\n",
    "nscw=3\n",
    "chi2_limit=1.2\n",
    "ng_sig_limit=3.\n",
    "systematic_fraction=0.01\n",
    "freeze_grbm_tem=-1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_name=\"Crab\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import matplotlib.pylab as plt\n",
    "\n",
    "\n",
    "import numpy as np\n",
    "import importlib\n",
    "from astroquery.simbad import Simbad\n",
    "from astropy import coordinates as coord\n",
    "\n",
    "from IPython.display import Image\n",
    "from IPython.display import display \n",
    "import xspec\n",
    "import shutil\n",
    "\n",
    "import os\n",
    "import shutil\n",
    "import oda\n",
    "\n",
    "from scipy import stats\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if reference_instrument.lower() == 'spi':\n",
    "    r=oda.evaluate(\"kb\", \"http://odahub.io/cc/isgri-oda-spi-reference\",\n",
    "               subcases_pattern=subcases_pattern,\n",
    "               reference_location=os.getcwd(), \n",
    "               source_name=source_name, \n",
    "               osa_version=osa_version,\n",
    "               nscw=nscw,\n",
    "               _ntries=1)\n",
    "elif reference_instrument.lower() == 'nustar':\n",
    "    #r=odakb.evaluate(\"git@gitlab.astro.unige.ch:integral/cc-workflows/cc-isgri-oda-nustar-reference.git\",\n",
    "    r=oda.evaluate(\"kb\", \"http://odahub.io/cc/isgri-oda-nustar-reference\",\n",
    "               subcases_pattern=subcases_pattern,\n",
    "               reference_location=os.getcwd(), \n",
    "               source_name=source_name, \n",
    "               osa_version=osa_version,nscw=nscw) \n",
    "               #,nbname=\"isgri-oda-nustar-reference\")\n",
    "else:\n",
    "    print(\"Unknown reference instrument\")\n",
    "    raise ValueError"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k,v in r.items():\n",
    "    if not k.endswith(\"_content\"):\n",
    "        print(k,v)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# I add response and arf keywords for easier reading\n",
    "import astropy.io.fits as pf\n",
    "\n",
    "if reference_instrument.lower() == 'spi':\n",
    "    ff=pf.open(r[\"reference_spec\"],  mode='update')\n",
    "    ff[2].header['RESPFILE']=r[\"reference_rmf\"]\n",
    "    ff.flush()\n",
    "    ff.close()\n",
    "elif reference_instrument.lower() == 'nustar':\n",
    "    for x in ['A', 'B']:\n",
    "        \n",
    "        ff=pf.open(r[\"reference_spec_\"+x],  mode='update')\n",
    "        ff[1].header['RESPFILE']=r[\"reference_rmf_\"+x]\n",
    "        ff[1].header['ANCRFILE']=r[\"reference_arf_\"+x]\n",
    "        ff[1].header['BACKFILE']=r[\"reference_bkg_\"+x]\n",
    "        ff.flush()\n",
    "        ff.close()\n",
    "else:\n",
    "    print(\"Undefined reference instrument\")\n",
    "    raise ValueError\n",
    "\n",
    "\n",
    "ff=pf.open(r[\"isgri_spec\"],  mode='update')\n",
    "ff[1].header['RESPFILE']=r[\"isgri_rmf\"]\n",
    "ff[1].header['ANCRFILE']=r[\"isgri_arf\"]\n",
    "ff.flush()\n",
    "ff.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xspec\n",
    "importlib.reload(xspec)\n",
    "\n",
    "#init dictionaries\n",
    "fit_by_lt = {}\n",
    "fn_by_lt={}\n",
    "\n",
    "xspec.AllModels.systematic=systematic_fraction   \n",
    "\n",
    "#Here, I do not care for the test about the low thresholds, I care about the comparison of the cyclotron line energy\n",
    "#I keep for consistency the loop over one element\n",
    "for c_emin in [15, 20, 28]: #np.linspace(17,40,5):\n",
    "    #freeze parameters apart from norm and E_cyc\n",
    "    \n",
    "    xspec.AllData.clear()\n",
    "\n",
    "    if reference_instrument == 'spi':\n",
    "        model='cflux*grbm'\n",
    "    elif reference_instrument == 'nustar':\n",
    "        model='cflux*grbm'\n",
    "            \n",
    "    \n",
    "    m1=xspec.Model(model)\n",
    "    \n",
    "    \n",
    "    if reference_instrument == 'spi':\n",
    "        xspec.AllData(r['isgri_spec']+\" 2:2 \" + r[\"reference_spec\"] )\n",
    "    elif reference_instrument == 'nustar':\n",
    "        m1(6).link=\"1\"\n",
    "        m1(7).values = 2.0\n",
    "        m1(8).values = 1.0\n",
    "        xspec.AllData(r['isgri_spec']+\" 2:2 \" + r[\"reference_spec_A\"] + \" 3:3 \" + r[\"reference_spec_B\"] )\n",
    "    else:\n",
    "        raise ValueError\n",
    "    \n",
    "\n",
    "    s     = xspec.AllData(1)\n",
    "    s_ref = xspec.AllData(2)\n",
    "    \n",
    "    isgri = xspec.AllModels(1)\n",
    "    ref   = xspec.AllModels(2)\n",
    "    \n",
    "    print(m1.nParameters)\n",
    "    \n",
    "    xspec.AllData.ignore('bad')\n",
    "    xspec.AllData.ignore('500.0-**')\n",
    "    \n",
    "    ig=\"**-%.2f,500.-**\"%c_emin\n",
    "    print(ig)\n",
    "    s.ignore(ig)\n",
    "    \n",
    "    if reference_instrument == 'nustar':\n",
    "        ig=\"**-3.0,70.0-**\"\n",
    "        print(ig)\n",
    "        s_ref.ignore(ig)\n",
    "        \n",
    "        s_ref_b=xspec.AllData(3)\n",
    "        ref_b=xspec.AllModels(3)\n",
    "        \n",
    "        s_ref_b.ignore(ig)\n",
    "    \n",
    "    #Key for output\n",
    "    lt_key='%.10lg'%c_emin\n",
    "\n",
    "        \n",
    "    #ref.gabs.Strength.link = \"%d\"%(ref_b.gabs.Strength.index)\n",
    "    \n",
    "    \n",
    "    ref.cflux.lg10Flux.untie()\n",
    "    ref.cflux.lg10Flux.frozen=False\n",
    "    \n",
    "    isgri.cflux.lg10Flux=-8\n",
    "    ref.cflux.lg10Flux=-8\n",
    "    \n",
    "    if reference_instrument == 'nustar':\n",
    "        ref_b.cflux.lg10Flux.untie()\n",
    "        ref_b.cflux.lg10Flux.frozen=False\n",
    "        ref_b.gabs.LineE.link=\"%d\"%(m1.nParameters+ref_b.gabs.LineE.index)\n",
    "        ref_b.cflux.lg10Flux=-8\n",
    "    \n",
    "    \n",
    "    isgri.cflux.Emin=20.\n",
    "    isgri.cflux.Emax=80.\n",
    "    isgri.grbm.norm=1.0\n",
    "    isgri.grbm.norm.frozen=True\n",
    "    isgri.grbm.tem=600 if freeze_grbm_tem < 0 else freeze_grbm_tem\n",
    "    isgri.grbm.tem.frozen=True\n",
    "        \n",
    "    xspec.Fit.perform()\n",
    "    \n",
    "    isgri.grbm.tem.frozen = freeze_grbm_tem > 0\n",
    "    xspec.Fit.perform()\n",
    "    \n",
    "    if reference_instrument == 'nustar':\n",
    "        isgri.highecut.foldE.frozen=False\n",
    "        isgri.highecut.cutoffE.frozen=False\n",
    "    \n",
    "    xspec.Fit.perform()\n",
    "        \n",
    "    xspec.Fit.error(\"1.0 max 5.0 1-%d\"%(2*m1.nParameters))\n",
    "\n",
    "    models={}\n",
    "    \n",
    "    fit_by_lt[lt_key]=dict(\n",
    "            emin=c_emin,\n",
    "            chi2_red=xspec.Fit.statistic/xspec.Fit.dof,                                \n",
    "            chi2=xspec.Fit.statistic,\n",
    "            ndof=xspec.Fit.dof,                                    \n",
    "            models=models,\n",
    "        )\n",
    "        \n",
    "    for m,ss in (isgri, 'isgri'), (ref, 'ref'):\n",
    "        #initialize dictionaries\n",
    "        models[ss]={}\n",
    "        #models[ss]['flux']={}\n",
    "        \n",
    "        #fills dictionaries\n",
    "        for i in range(1,m.nParameters+1): \n",
    "            if (not m(i).frozen) and (not bool(m(i).link)):\n",
    "                #use the name plus position because there could be parameters with same name from multiple \n",
    "                #model components (e.g., several gaussians)\n",
    "                print(m(i).name, \"%.2f\"%(m(i).values[0]), m(i).frozen,bool(m(i).link) )\n",
    "                models[ss][m(i).name+\"_%02d\"%(i)]=[ m(i).values[0], m(i).error[0], m(i).error[1] ]              \n",
    "        \n",
    "        \n",
    "\n",
    "\n",
    "#     for flux_e1, flux_e2 in [(30,80), (80,200)]:\n",
    "#         xspec.AllModels.calcFlux(\"{} {} err\".format(flux_e1, flux_e2))\n",
    "#         #print ( xspec.AllData(1).flux)\n",
    "#         models['isgri']['flux'][(flux_e1, flux_e2)] = xspec.AllData(1).flux\n",
    "#         models['ref']['flux'][(flux_e1, flux_e2)] = xspec.AllData(2).flux\n",
    "\n",
    "        \n",
    "\n",
    "    xspec.Plot.device=\"/png\"\n",
    "    #xspec.Plot.addCommand(\"setplot en\")\n",
    "    xspec.Plot.xAxis=\"keV\"\n",
    "    xspec.Plot(\"ldata del\")\n",
    "    xspec.Plot.device=\"/png\"\n",
    "    \n",
    "    fn=\"fit_lt%.5lg.png\"%c_emin\n",
    "    \n",
    "    fn_by_lt[lt_key] = fn\n",
    "    \n",
    "    shutil.move(\"pgplot.png_2\", fn)\n",
    "\n",
    "    _=display(Image(filename=fn,format=\"png\"))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# marginally more adcanced choice of low thresholds\n",
    "print(\"lt\\tchi2_r\\tprob\\tsigmas\")\n",
    "for lt,d in fit_by_lt.items():\n",
    "    d['nh_prob']=stats.chi2(d['ndof']).sf(d['chi2'])\n",
    "    d['nh_sig']=stats.norm().isf(d['nh_prob'])\n",
    "    print(\"%.1f\\t%.2f\\t%.2f\\t%.2f\"%(float(lt), d['chi2_red'], d['nh_prob'], d['nh_sig']) )\n",
    "    \n",
    "\n",
    "acceptable = [p for p in fit_by_lt.items() if p[1]['nh_sig']<ng_sig_limit]\n",
    "\n",
    "if len(acceptable)>0:\n",
    "    good_lt = min(acceptable, key=lambda x:x)\n",
    "    good_lt_next = min(acceptable, key=lambda x:x)\n",
    "else:\n",
    "    good_lt = None\n",
    "    good_lt_next = None\n",
    "    print(\"WARNING: no acceptable LT found!\")\n",
    "\n",
    "\n",
    "good_lt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parameter_comparison={}\n",
    "#compare parameters of best fit\n",
    "\n",
    "if good_lt:\n",
    "    best_result=good_lt[1]\n",
    "    for par_name, par_values_ref in best_result['models']['ref'].items():\n",
    "        par_values_isgri=best_result['models']['isgri'][par_name]\n",
    "        isgri_value= par_values_isgri[0]\n",
    "        ref_value  = par_values_ref[0]\n",
    "\n",
    "        if ref_value > isgri_value:\n",
    "            err_ref  =  ref_value - par_values_ref[1]\n",
    "            err_isgri = par_values_isgri[2] - isgri_value\n",
    "        else:\n",
    "            err_ref  =  - (ref_value - par_values_ref[2] )\n",
    "            err_isgri = - (par_values_isgri[1] - isgri_value )\n",
    "\n",
    "        par_diff_sigmas = (isgri_value - ref_value) /np.sqrt(err_ref**2 + err_isgri**2)\n",
    "\n",
    "        print(par_name+\" %.2f +/- %.2f ; %.2f +/- %.2f ; %.1f\"%(isgri_value, err_isgri, ref_value, err_ref, par_diff_sigmas))\n",
    "\n",
    "        #We do not compare the normalization with NuSTAR strictly, because of non-simultaneity issues\n",
    "        if ('log10Flux' in par_name) and (reference_instrument == 'nustar') :\n",
    "            success=True\n",
    "        else:\n",
    "            success = bool(np.abs(par_diff_sigmas) < ng_sig_limit)\n",
    "\n",
    "\n",
    "        parameter_comparison[par_name] = dict(\n",
    "\n",
    "            ref = [ref_value, err_ref],\n",
    "            isgri = [isgri_value, err_isgri ],\n",
    "\n",
    "            significance = np.abs(par_diff_sigmas),\n",
    "\n",
    "            success=success\n",
    "        )        \n",
    "else:\n",
    "        parameter_comparison={}\n",
    "        \n",
    "parameter_comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if good_lt:\n",
    "    fn_picked=fn_by_lt[good_lt[0]]\n",
    "    fn_next_picked=fn_by_lt[good_lt_next[0]]\n",
    "else:\n",
    "    fn_picked=\"\"\n",
    "    fn_next_picked=\"\"\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary={\n",
    "    'status': 'OK' if good_lt else \"NOK\",\n",
    "    'emin_good': good_lt,\n",
    "    'parameter_comparison': parameter_comparison\n",
    "}\n",
    "summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "outputs"
    ]
   },
   "outputs": [],
   "source": [
    "summary=summary\n",
    "fit_results=fit_by_lt\n",
    "good_fit_png=fn_picked\n",
    "next_good_fit_png=fn_next_picked"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}