{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PUM actuation during O2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from gwpy.timeseries import TimeSeriesDict\n",
    "from gwpy.segments import DataQualityFlag\n",
    "from gwpy.spectrogram import SpectrogramList\n",
    "from gwpy.plotter import FrequencySeriesPlot\n",
    "from gwpy.frequencyseries import FrequencySeries\n",
    "import numpy as np\n",
    "import scipy.signal as sig\n",
    "import sys\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up for requesting PUM MASTER_OUT data for all IFOs, QUADs, and quadrants during O2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "ifos = ['H1', 'L1']\n",
    "optics = ['ITMX', 'ITMY', 'ETMX', 'ETMY']\n",
    "stages = ['L2',]\n",
    "coils = ['UL', 'UR', 'LL', 'LR']\n",
    "science_seg = '{ifo}:DMT-ANALYSIS_READY:1'\n",
    "mochannelpat = '{ifo}:SUS-{optic}_{stage}_MASTER_OUT_{coil}_DQ'\n",
    "nmchannelpat = '{ifo}:SUS-{optic}_{stage}_NOISEMON_{coil}_OUT_DQ'\n",
    "stchannelpat = '{ifo}:SUS-{optic}_BIO_{stage}_{coil}_STATEREQ'\n",
    "# O2 = 1164556817-1187733618\n",
    "start = 1164556817\n",
    "end = 1187733618"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uncomment this to get the data (a sample of 1% of segments to keep it manageable to start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Segment(1164559438.0, 1164559657.0), Segment(1167860943.0, 1167942130.0), Segment(1170782624.0, 1170783827.0), Segment(1173583963.0, 1173593679.0), Segment(1176705963.0, 1176739235.0), Segment(1182717707.0, 1182726561.0), Segment(1186642718.0, 1186664343.0)]\n",
      "[1164559438.0 ... 1164559657.0)\n",
      "[1167860943.0 ... 1167942130.0)\n",
      "[1170782624.0 ... 1170783827.0)\n",
      "[1173583963.0 ... 1173593679.0)\n",
      "[1176705963.0 ... 1176739235.0)\n",
      "[1182717707.0 ... 1182726561.0)\n",
      "[1186642718.0 ... 1186664343.0)\n",
      "[Segment(1164600796.0, 1164608249.0), Segment(1165990476.0, 1165990565.0), Segment(1168780605.0, 1168784610.0), Segment(1171030725.0, 1171102349.0), Segment(1173577169.0, 1173579561.0), Segment(1176021152.0, 1176021200.0), Segment(1180312587.0, 1180358875.0), Segment(1183488007.0, 1183491392.0), Segment(1186531497.0, 1186561209.0)]\n",
      "[1164600796.0 ... 1164608249.0)\n",
      "[1165990476.0 ... 1165990565.0)\n",
      "[1168780605.0 ... 1168784610.0)\n",
      "[1171030725.0 ... 1171102349.0)\n",
      "[1173577169.0 ... 1173579561.0)\n",
      "[1176021152.0 ... 1176021200.0)\n",
      "[1180312587.0 ... 1180358875.0)\n",
      "[1186531497.0 ... 1186561209.0)\n"
     ]
    }
   ],
   "source": [
    "# data = dict()\n",
    "# for ifo in ifos:\n",
    "#     science = DataQualityFlag.query(science_seg.format(ifo=ifo), start, end)\n",
    "#     segselect = [seg for n, seg in enumerate(science.active) if n % 100 == 0]\n",
    "#     print segselect\n",
    "#     mochannellist = [mochannelpat.format(ifo=ifo, optic=optic, stage=stage, coil=coil)\n",
    "#                      for optic in optics for stage in stages for coil in coils]\n",
    "#     nmchannellist = [nmchannelpat.format(ifo=ifo, optic=optic, stage=stage, coil=coil)\n",
    "#                      for optic in optics for stage in stages for coil in coils]\n",
    "#     stchannellist = [stchannelpat.format(ifo=ifo, optic=optic, stage=stage, coil=coil)\n",
    "#                      for optic in optics for stage in stages for coil in coils]\n",
    "#     data[ifo] = [TimeSeriesDict.get(mochannellist+stchannellist, seg[0], seg[1], frametype='{ifo}_R'.format(ifo=ifo))\n",
    "#                  for seg in segselect if sys.stdout.write(str(seg)+\"\\n\") or True]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uncomment to write time series data to disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for ifo in ifos:\n",
    "#     for n, tsd in enumerate(data[ifo]):\n",
    "#         tsd.write(\"{ifo}_{n}.hdf5\".format(ifo=ifo, n=n))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load time series data from disk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = dict()\n",
    "files = [file for file in os.listdir('.') if file.endswith('.hdf5')]\n",
    "for ifo in ifos:\n",
    "    for file in files:\n",
    "        if not file.startswith(ifo):\n",
    "            continue\n",
    "        data.setdefault(ifo, []).append(TimeSeriesDict.read(file))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute and plot the 99th percentile spectrum from each channel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "p99 = dict()\n",
    "for ifo in ifos:\n",
    "    for optic in optics:\n",
    "        for stage in stages:\n",
    "            for coil in coils:\n",
    "                mochan = mochannelpat.format(ifo=ifo, optic=optic, stage=stage, coil=coil)\n",
    "                data_list = [segdict[mochan] for segdict in data[ifo]]\n",
    "                spec_list = [ts.spectrogram2(fftlength=16, overlap=8) ** (1./2) for ts in data_list]\n",
    "                specgrams = SpectrogramList(*spec_list)\n",
    "                stack_specgrams = specgrams.join(gap='ignore')\n",
    "                p99[mochan] = stack_specgrams.percentile(99)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "p99 = dict()\n",
    "p50 = dict()\n",
    "p01 = dict()\n",
    "for ifo in ifos:\n",
    "    for optic in optics:\n",
    "        for stage in stages:\n",
    "            for coil in coils:\n",
    "                mochan = mochannelpat.format(ifo=ifo, optic=optic, stage=stage, coil=coil)\n",
    "                data_list = [segdict[mochan] for segdict in data[ifo]]\n",
    "                spec_list = [ts.spectrogram2(fftlength=16, overlap=8) ** (1./2) for ts in data_list]\n",
    "                specgrams = SpectrogramList(*spec_list)\n",
    "                stack_specgrams = specgrams.join(gap='ignore')\n",
    "                p99[mochan] = stack_specgrams.percentile(99)\n",
    "                p50[mochan] = stack_specgrams.percentile(50)\n",
    "                p01[mochan] = stack_specgrams.percentile(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H1:SUS-ITMX_BIO_L2_UL_STATEREQ [3, 3, 3, 3, 3, 3, 3]\n",
      "H1:SUS-ITMX_BIO_L2_UR_STATEREQ [3, 3, 3, 3, 3, 3, 3]\n",
      "H1:SUS-ITMX_BIO_L2_LL_STATEREQ [3, 3, 3, 3, 3, 3, 3]\n",
      "H1:SUS-ITMX_BIO_L2_LR_STATEREQ [3, 3, 3, 3, 3, 3, 3]\n",
      "H1:SUS-ITMY_BIO_L2_UL_STATEREQ [3, 3, 3, 3, 3, 3, 3]\n",
      "H1:SUS-ITMY_BIO_L2_UR_STATEREQ [3, 3, 3, 3, 3, 3, 3]\n",
      "H1:SUS-ITMY_BIO_L2_LL_STATEREQ [3, 3, 3, 3, 3, 3, 3]\n",
      "H1:SUS-ITMY_BIO_L2_LR_STATEREQ [3, 3, 3, 3, 3, 3, 3]\n",
      "H1:SUS-ETMX_BIO_L2_UL_STATEREQ [3, 3, 3, 3, 3, 3, 3]\n",
      "H1:SUS-ETMX_BIO_L2_UR_STATEREQ [3, 3, 3, 3, 3, 3, 3]\n",
      "H1:SUS-ETMX_BIO_L2_LL_STATEREQ [3, 3, 3, 3, 3, 3, 3]\n",
      "H1:SUS-ETMX_BIO_L2_LR_STATEREQ [3, 3, 3, 3, 3, 3, 3]\n",
      "H1:SUS-ETMY_BIO_L2_UL_STATEREQ [3, 3, 3, 3, 3, 3, 3]\n",
      "H1:SUS-ETMY_BIO_L2_UR_STATEREQ [3, 3, 3, 3, 3, 3, 3]\n",
      "H1:SUS-ETMY_BIO_L2_LL_STATEREQ [3, 3, 3, 3, 3, 3, 3]\n",
      "H1:SUS-ETMY_BIO_L2_LR_STATEREQ [3, 3, 3, 3, 3, 3, 3]\n",
      "L1:SUS-ITMX_BIO_L2_UL_STATEREQ [1, 1, 1, 1, 1, 3, 3, 1, 1]\n",
      "L1:SUS-ITMX_BIO_L2_UR_STATEREQ [1, 1, 1, 1, 1, 3, 3, 1, 1]\n",
      "L1:SUS-ITMX_BIO_L2_LL_STATEREQ [1, 1, 1, 1, 1, 3, 3, 1, 1]\n",
      "L1:SUS-ITMX_BIO_L2_LR_STATEREQ [1, 1, 1, 1, 1, 3, 3, 1, 1]\n",
      "L1:SUS-ITMY_BIO_L2_UL_STATEREQ [1, 1, 1, 1, 1, 1, 1, 1, 1]\n",
      "L1:SUS-ITMY_BIO_L2_UR_STATEREQ [1, 1, 1, 1, 1, 1, 1, 1, 1]\n",
      "L1:SUS-ITMY_BIO_L2_LL_STATEREQ [1, 1, 1, 1, 1, 1, 1, 1, 1]\n",
      "L1:SUS-ITMY_BIO_L2_LR_STATEREQ [1, 1, 1, 1, 1, 1, 1, 1, 1]\n",
      "L1:SUS-ETMX_BIO_L2_UL_STATEREQ [1, 1, 1, 1, 1, 3, 3, 1, 1]\n",
      "L1:SUS-ETMX_BIO_L2_UR_STATEREQ [1, 1, 1, 1, 1, 3, 3, 1, 1]\n",
      "L1:SUS-ETMX_BIO_L2_LL_STATEREQ [1, 1, 1, 1, 1, 3, 3, 1, 1]\n",
      "L1:SUS-ETMX_BIO_L2_LR_STATEREQ [1, 1, 1, 1, 1, 3, 3, 1, 1]\n",
      "L1:SUS-ETMY_BIO_L2_UL_STATEREQ [1, 1, 1, 1, 1, 1, 1, 1, 1]\n",
      "L1:SUS-ETMY_BIO_L2_UR_STATEREQ [1, 1, 1, 1, 1, 1, 1, 1, 1]\n",
      "L1:SUS-ETMY_BIO_L2_LL_STATEREQ [1, 1, 1, 1, 1, 1, 1, 1, 1]\n",
      "L1:SUS-ETMY_BIO_L2_LR_STATEREQ [1, 1, 1, 1, 1, 1, 1, 1, 1]\n"
     ]
    }
   ],
   "source": [
    "for ifo in ifos:\n",
    "    for optic in optics:\n",
    "        for stage in stages:\n",
    "            for coil in coils:\n",
    "                stchan = stchannelpat.format(ifo=ifo, optic=optic, stage=stage, coil=coil)\n",
    "                data_list = [segdict[stchan] for segdict in data[ifo]]\n",
    "                print stchan, [int(d.mean().value) for d in data_list]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python2.7/site-packages/ipykernel_launcher.py:2: RuntimeWarning: divide by zero encountered in divide\n",
      "  \n"
     ]
    }
   ],
   "source": [
    "f = p99.values()[0].frequencies.value\n",
    "dacnoise = 1e-9*np.sqrt(300.**2*(50/f)**2 + 300.**2)\n",
    "dacnoise = FrequencySeries(dacnoise, frequencies=f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fplotH1 = FrequencySeriesPlot([(((p99[p]*20./2**18)**2)[::-1].cumsum()[::-1])**0.5 for p in p99 if p.startswith('H1')]+[dacnoise*200*np.abs(zpkff),], figsize=(12,4))\n",
    "fplotH1 = FrequencySeriesPlot([p99[p]*20./2**18 for p in sorted(p99) if p.startswith('H1')]+[dacnoise], figsize=(10,8))\n",
    "fplotH1.set_ylabel('ASD [\\,V/$\\sqrt{\\mathrm{Hz}}$\\,]')\n",
    "fplotH1.set_ylim([1e-7, 1e1])\n",
    "fplotH1.set_title('H1 SUS QUAD PUM actuation (99th percentile)')\n",
    "plt.legend([p.split(': ')[0].replace('_',r'\\_') for p in sorted(p99) if p.startswith('H1')] + ['DAC noise model'], bbox_to_anchor=(1,1.015), loc=\"upper left\")\n",
    "plt.draw()\n",
    "plt.savefig('H1_99.pdf', bbox_inches='tight')\n",
    "#plt.show()\n",
    "#fplotL1 = FrequencySeriesPlot([(((p99[p]*20./2**18)**2)[::-1].cumsum()[::-1])**0.5 for p in p99 if p.startswith('L1')]+[dacnoise*200*np.abs(zpkff),], figsize=(12,4))\n",
    "fplotL1 = FrequencySeriesPlot([p99[p]*20./2**18 for p in sorted(p99) if p.startswith('L1')]+[dacnoise], figsize=(10,8))\n",
    "fplotL1.set_ylabel('ASD [\\,V/$\\sqrt{\\mathrm{Hz}}$\\,]')\n",
    "fplotL1.set_ylim([1e-7, 1e1])\n",
    "fplotL1.set_title('L1 SUS QUAD PUM actuation (99th percentile)')\n",
    "plt.legend([p.split(': ')[0].replace('_',r'\\_') for p in sorted(p99) if p.startswith('L1')] + ['DAC noise model'], bbox_to_anchor=(1,1.015), loc=\"upper left\")\n",
    "plt.draw()\n",
    "plt.savefig('L1_99.pdf', bbox_inches='tight')\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "fplotH1 = FrequencySeriesPlot([p50[p]*20./2**18 for p in sorted(p50) if p.startswith('H1')]+[dacnoise], figsize=(10,8))\n",
    "fplotH1.set_ylabel('ASD [\\,V/$\\sqrt{\\mathrm{Hz}}$\\,]')\n",
    "fplotH1.set_ylim([1e-7, 1e1])\n",
    "fplotH1.set_title('H1 SUS QUAD PUM actuation (50th percentile)')\n",
    "plt.legend([p.split(': ')[0].replace('_',r'\\_') for p in sorted(p50) if p.startswith('H1')] + ['DAC noise model'], bbox_to_anchor=(1,1.015), loc=\"upper left\")\n",
    "plt.draw()\n",
    "plt.savefig('H1_50.pdf', bbox_inches='tight')\n",
    "#plt.show()\n",
    "fplotL1 = FrequencySeriesPlot([p50[p]*20./2**18 for p in sorted(p50) if p.startswith('L1')]+[dacnoise], figsize=(10,8))\n",
    "fplotL1.set_ylabel('ASD [\\,V/$\\sqrt{\\mathrm{Hz}}$\\,]')\n",
    "fplotL1.set_ylim([1e-7, 1e1])\n",
    "fplotL1.set_title('L1 SUS QUAD PUM actuation (50th percentile)')\n",
    "plt.legend([p.split(': ')[0].replace('_',r'\\_') for p in sorted(p50) if p.startswith('L1')] + ['DAC noise model'], bbox_to_anchor=(1,1.015), loc=\"upper left\")\n",
    "plt.draw()\n",
    "plt.savefig('L1_50.pdf', bbox_inches='tight')\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "fplotH1 = FrequencySeriesPlot([p01[p]*20./2**18 for p in sorted(p01) if p.startswith('H1')]+[dacnoise], figsize=(10,8))\n",
    "fplotH1.set_ylabel('ASD [\\,V/$\\sqrt{\\mathrm{Hz}}$\\,]')\n",
    "fplotH1.set_ylim([1e-7, 1e1])\n",
    "fplotH1.set_title('H1 SUS QUAD PUM actuation (1st percentile)')\n",
    "plt.legend([p.split(': ')[0].replace('_',r'\\_') for p in sorted(p01) if p.startswith('H1')] + ['DAC noise model'], bbox_to_anchor=(1,1.015), loc=\"upper left\")\n",
    "plt.draw()\n",
    "plt.savefig('H1_01.pdf', bbox_inches='tight')\n",
    "#plt.show()\n",
    "fplotL1 = FrequencySeriesPlot([p01[p]*20./2**18 for p in sorted(p01) if p.startswith('L1')]+[dacnoise], figsize=(10,8))\n",
    "fplotL1.set_ylabel('ASD [\\,V/$\\sqrt{\\mathrm{Hz}}$\\,]')\n",
    "fplotL1.set_ylim([1e-7, 1e1])\n",
    "fplotL1.set_title('L1 SUS QUAD PUM actuation (1st percentile)')\n",
    "plt.legend([p.split(': ')[0].replace('_',r'\\_') for p in sorted(p01) if p.startswith('L1')] + ['DAC noise model'], bbox_to_anchor=(1,1.015), loc=\"upper left\")\n",
    "plt.draw()\n",
    "plt.savefig('L1_01.pdf', bbox_inches='tight')\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2.7.5 [Long Term - SL7]",
   "language": "python",
   "name": "python2-lta"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}