Zr_DESY_2020_close.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DESY depth profiling data analysis workbook [2020]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a workbook used to integrate synchrotron data from DESY then generate the bg-spline and peak-index files required for running CMWP."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import stuff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib qt\n",
    "import os\n",
    "import glob\n",
    "import pyFAI, pyFAI.azimuthalIntegrator\n",
    "import csv\n",
    "from shutil import copyfile\n",
    "\n",
    "from src.cmwp_tools import load_tifs\n",
    "from src.xrd_tools import getReflections\n",
    "from src.desy_tools import load_fio\n",
    "from src.cmwp_tools import getBaseline, getPeaks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wavelength: 0.18505104 A\n"
     ]
    }
   ],
   "source": [
    "### Experiment specific settings\n",
    "\n",
    "base = \"/mnt/manchester_rds/202011_DESY/\"                            # This is the base directory where the data is stored\n",
    "\n",
    "calib_file = base + 'raw/close_detector/calib_new.poni'             # These are the calibration files (.poni) for the detectors\n",
    "\n",
    "ai = pyFAI.load(calib_file)\n",
    "wavelength = ai.wavelength * 1e10\n",
    "print('Wavelength: {0:.8f} A'.format(wavelength))\n",
    "\n",
    "intpoints = 2000                                                  # This is how many integration points to use for pyFAI\n",
    "\n",
    "## Zr peaks\n",
    "peak_name, peak_pos = getReflections(crystalType='hcp', a=3.232, c=5.147, wavelength=0.1848, printReflections=False)\n",
    "\n",
    "## SPP peaks\n",
    "spp_peak_name, spp_peak_pos = getReflections(crystalType='hcp', a=5.028, c=8.248, wavelength=0.1848, printReflections=False)\n",
    "spp_peak_name = spp_peak_name[[5,7,8,9,26]]\n",
    "spp_peak_pos = spp_peak_pos[[5,7,8,9,26]]\n",
    "\n",
    "## hydride peaks\n",
    "delta_peak_name, delta_peak_pos = getReflections(crystalType='fcc', a=4.777, wavelength=0.1848, printReflections=False)\n",
    "delta_peak_name = delta_peak_name[[0,2]]\n",
    "delta_peak_pos = delta_peak_pos[[0,2]]\n",
    "\n",
    "## Background\n",
    "baseline=[3.1, 3.5, 4.7, 5.2, 5.9, 6.9, 8.8, 9.3, 11.55]             # What 2theta points to calculte the background spline from\n",
    "baseline_interpolate=[4.1, 7.85, 10.475]                            # What 2theta points to add in extra point\n",
    "baseline_interpolate_factor=[0.97, 1, 1.15]\n",
    "\n",
    "limits=[3, 11.6]                                                       # 2theta bounds of integration\n",
    "\n",
    "searchrange = int(0.005 * intpoints)                         # How many data points each side of the approximate 2theta peak position to search for the true peak\n",
    "\n",
    "method = 'full_csr'\n",
    "azimuth_range=None\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### 0.1 dpa hydrided\n",
    "               \n",
    "directory = base + \"raw/close_detector/hyd/hyd/hyd/\"                                  \n",
    "outputdir = \"/home/rhys/Documents/CMWP-211102/2020_11_DESY/hyd/\"       # This is the output directory (normally a sub-folder in your CMWP dir)\n",
    "fio_file = base + \"raw/close_detector/hyd/eh3scan1_00144.fio\"\n",
    "darks = base + \"raw/close_detector/hyd/dark\"\n",
    "templates = base + \"templates/hyd/\"                                      # This is where the template files are stored. These are copied for each integration.\n",
    "z_range = (14.916, 14.976)\n",
    "y_range = None\n",
    "spp=False; hyd=False;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### 0.1 dpa hydrided **3phase*\n",
    "               \n",
    "directory = base + \"raw/close_detector/hyd/hyd/hyd/\"                                  \n",
    "outputdir = \"/home/rhys/Documents/CMWP-211102/2020_11_DESY/hyd_3p/\"       # This is the output directory (normally a sub-folder in your CMWP dir)\n",
    "fio_file = base + \"raw/close_detector/hyd/eh3scan1_00144.fio\"\n",
    "darks = base + \"raw/close_detector/hyd/dark\"\n",
    "templates = base + \"templates/hyd_3phase/\"                                      # This is where the template files are stored. These are copied for each integration.\n",
    "z_range = (14.916, 14.976)\n",
    "y_range = None\n",
    "spp=True; hyd=True;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### 0.1 dpa hydrided **3phase* *limazim*\n",
    "               \n",
    "directory = base + \"raw/close_detector/hyd/hyd/hyd/\"                                  \n",
    "outputdir = \"/home/rhys/Documents/CMWP-211102/2020_11_DESY/hyd_3p_80_100/\"       # This is the output directory (normally a sub-folder in your CMWP dir)\n",
    "fio_file = base + \"raw/close_detector/hyd/eh3scan1_00144.fio\"\n",
    "darks = base + \"raw/close_detector/hyd/dark\"\n",
    "templates = base + \"templates/hyd_3phase_bgfit/\"                                      # This is where the template files are stored. These are copied for each integration.\n",
    "z_range = (14.916, 14.976)\n",
    "y_range = None\n",
    "spp=True; hyd=True;\n",
    "\n",
    "intpoints=1500\n",
    "\n",
    "method = 'IntegrationMethod(2d int, pseudo split, histogram, cython)'\n",
    "azimuth_range=[[70,110], [250, 290]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### 2 dpa\n",
    "\n",
    "fio_file = base + \"raw/close_detector/dpa2/eh3scan1_00096.fio\"\n",
    "directory = base + \"raw/close_detector/dpa2/dpa2/\"\n",
    "outputdir = \"/home/rhys/Documents/CMWP-211102/2020_11_DESY/dpa2/\"\n",
    "darks = base + \"raw/close_detector/dpa2/dark\"\n",
    "templates = base + \"templates/2dpa/\"                                      # This is where the template files are stored. These are copied for each integration.\n",
    "z_range = (14.830, 14.890)\n",
    "y_range = None\n",
    "spp=False; hyd=False;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### 2 dpa 2p\n",
    "\n",
    "fio_file = base + \"raw/close_detector/dpa2/eh3scan1_00096.fio\"\n",
    "directory = base + \"raw/close_detector/dpa2/dpa2/\"\n",
    "outputdir = \"/home/rhys/Documents/CMWP-211102/2020_11_DESY/dpa2_2p/\"\n",
    "darks = base + \"raw/close_detector/dpa2/dark\"\n",
    "templates = base + \"templates/hyd_3phase_bgfit/\"                                         # This is where the template files are stored. These are copied for each integration.\n",
    "z_range = (14.830, 14.890)\n",
    "y_range = None\n",
    "spp=True; hyd=False;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### 1 dpa\n",
    "\n",
    "fio_file = base + \"raw/close_detector/dpa1/eh3scan1_00077.fio\"\n",
    "directory = base + \"raw/close_detector/dpa1/dpa1/\"\n",
    "outputdir = \"/home/rhys/Documents/CMWP-211102/2020_11_DESY/dpa1/\"\n",
    "darks = base + \"raw/close_detector/dpa1/dark\"\n",
    "templates = base + \"templates/2dpa/\"  \n",
    "z_range = (14.810, 14.870)\n",
    "y_range = (3.2, 4.3)\n",
    "spp=False; hyd=False;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fio_file = base + \"raw/close_detector/dpa05/dpa05_02.fio\"\n",
    "directory = base + \"raw/close_detector/dpa05/dpa05_02/\"\n",
    "outputdir = \"/home/rhys/Documents/CMWP-211102/2020_11_DESY/dpa05/\"\n",
    "darks = base + \"raw/close_detector/dpa05/dark\"\n",
    "templates = base + \"templates/2dpa/\"   \n",
    "z_range = (14.780, 14.860)\n",
    "y_range = None\n",
    "spp=False; hyd=False;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "fio_file = base + \"raw/close_detector/dpa05/dpa05_02.fio\"\n",
    "directory = base + \"raw/close_detector/dpa05/dpa05_02/\"\n",
    "outputdir = \"/home/rhys/Documents/CMWP-211102/2020_11_DESY/dpa05_2p/\"\n",
    "darks = base + \"raw/close_detector/dpa05/dark\"\n",
    "templates = base + \"templates/hyd_3phase_bgfit/\"\n",
    "z_range = (14.780, 14.860)\n",
    "y_range = (-5.7, -3.7)\n",
    "spp=True; hyd=False;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Integrate lab6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#lab6 = glob.glob('/mnt/manchester_rds/202011_DESY/raw/close_detector/LaB6/LaB6_00*.tif')[2:]\n",
    "#lab6_dark = glob.glob('/mnt/manchester_rds/202011_DESY/raw/close_detector/LaB6/dark*.tif')\n",
    "\n",
    "#ai = pyFAI.load(calib_file)\n",
    "\n",
    "#data=load_tifs(lab6)-load_tifs(lab6_dark)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##output = ai.integrate1d(data, npt=intpoints, radial_range=[2.2,11.3],\n",
    "        correctSolidAngle=True, method='full_csr',polarization_factor=0.99, unit='2th_deg')\n",
    "#np.savetxt('/mnt/manchester_rds/202011_DESY/raw/close_detector/LaB6.dat', np.array(output).T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#output = ai.integrate1d(data, npt=intpoints, radial_range=[2.2,11.3], azimuth_range=[90,110],\n",
    "#        method='IntegrationMethod(2d int, pseudo split, histogram, cython)', correctSolidAngle=True, polarization_factor=0.99, unit='2th_deg')\n",
    "#np.savetxt('/mnt/manchester_rds/202011_DESY/raw/close_detector/LaB6_limazim.dat', np.array(output).T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read log file and make pandas table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The .fio file contains a table of motor positions (ie idtz2, idty1), image names and whether the image was a clearing frame or actual exposure. The section reads in the file into a Pandas dataframe for ease of use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = load_fio(filename=fio_file)\n",
    "\n",
    "if z_range is not None:\n",
    "    df = df[np.logical_and(df['idtz2']>z_range[0], df['idtz2']<z_range[1])]\n",
    "if y_range is not None:\n",
    "    df = df[(df['idty2(encoder)']>y_range[0]) & (df['idty2(encoder)']<y_range[1])]\n",
    "\n",
    "df.reset_index(inplace=True)\n",
    "df.drop('index', axis=1, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>idty2(encoder)</th>\n",
       "      <th>end pos</th>\n",
       "      <th>idtz2</th>\n",
       "      <th>channel</th>\n",
       "      <th>filename</th>\n",
       "      <th>type</th>\n",
       "      <th>unix time</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-4.99994</td>\n",
       "      <td>-4.99994</td>\n",
       "      <td>14.782</td>\n",
       "      <td>2</td>\n",
       "      <td>dpa05_00282.cbf</td>\n",
       "      <td>exposure</td>\n",
       "      <td>1.605998e+09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-4.49994</td>\n",
       "      <td>-4.49994</td>\n",
       "      <td>14.782</td>\n",
       "      <td>2</td>\n",
       "      <td>dpa05_00283.cbf</td>\n",
       "      <td>exposure</td>\n",
       "      <td>1.605998e+09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-3.99994</td>\n",
       "      <td>-3.99994</td>\n",
       "      <td>14.782</td>\n",
       "      <td>2</td>\n",
       "      <td>dpa05_00284.cbf</td>\n",
       "      <td>exposure</td>\n",
       "      <td>1.605998e+09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-3.99994</td>\n",
       "      <td>-3.99994</td>\n",
       "      <td>14.784</td>\n",
       "      <td>2</td>\n",
       "      <td>dpa05_00287.cbf</td>\n",
       "      <td>exposure</td>\n",
       "      <td>1.605998e+09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-4.49994</td>\n",
       "      <td>-4.49994</td>\n",
       "      <td>14.784</td>\n",
       "      <td>2</td>\n",
       "      <td>dpa05_00288.cbf</td>\n",
       "      <td>exposure</td>\n",
       "      <td>1.605998e+09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>131</th>\n",
       "      <td>-4.99994</td>\n",
       "      <td>-4.99994</td>\n",
       "      <td>14.856</td>\n",
       "      <td>2</td>\n",
       "      <td>dpa05_00469.cbf</td>\n",
       "      <td>exposure</td>\n",
       "      <td>1.605998e+09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>132</th>\n",
       "      <td>-5.49994</td>\n",
       "      <td>-5.49994</td>\n",
       "      <td>14.856</td>\n",
       "      <td>2</td>\n",
       "      <td>dpa05_00470.cbf</td>\n",
       "      <td>exposure</td>\n",
       "      <td>1.605998e+09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>133</th>\n",
       "      <td>-4.99994</td>\n",
       "      <td>-4.99994</td>\n",
       "      <td>14.858</td>\n",
       "      <td>2</td>\n",
       "      <td>dpa05_00472.cbf</td>\n",
       "      <td>exposure</td>\n",
       "      <td>1.605999e+09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>134</th>\n",
       "      <td>-4.49994</td>\n",
       "      <td>-4.49994</td>\n",
       "      <td>14.858</td>\n",
       "      <td>2</td>\n",
       "      <td>dpa05_00473.cbf</td>\n",
       "      <td>exposure</td>\n",
       "      <td>1.605999e+09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>135</th>\n",
       "      <td>-3.99994</td>\n",
       "      <td>-3.99994</td>\n",
       "      <td>14.858</td>\n",
       "      <td>2</td>\n",
       "      <td>dpa05_00474.cbf</td>\n",
       "      <td>exposure</td>\n",
       "      <td>1.605999e+09</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>136 rows × 7 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     idty2(encoder)  end pos   idtz2  channel         filename      type  \\\n",
       "0          -4.99994 -4.99994  14.782        2  dpa05_00282.cbf  exposure   \n",
       "1          -4.49994 -4.49994  14.782        2  dpa05_00283.cbf  exposure   \n",
       "2          -3.99994 -3.99994  14.782        2  dpa05_00284.cbf  exposure   \n",
       "3          -3.99994 -3.99994  14.784        2  dpa05_00287.cbf  exposure   \n",
       "4          -4.49994 -4.49994  14.784        2  dpa05_00288.cbf  exposure   \n",
       "..              ...      ...     ...      ...              ...       ...   \n",
       "131        -4.99994 -4.99994  14.856        2  dpa05_00469.cbf  exposure   \n",
       "132        -5.49994 -5.49994  14.856        2  dpa05_00470.cbf  exposure   \n",
       "133        -4.99994 -4.99994  14.858        2  dpa05_00472.cbf  exposure   \n",
       "134        -4.49994 -4.49994  14.858        2  dpa05_00473.cbf  exposure   \n",
       "135        -3.99994 -3.99994  14.858        2  dpa05_00474.cbf  exposure   \n",
       "\n",
       "        unix time  \n",
       "0    1.605998e+09  \n",
       "1    1.605998e+09  \n",
       "2    1.605998e+09  \n",
       "3    1.605998e+09  \n",
       "4    1.605998e+09  \n",
       "..            ...  \n",
       "131  1.605998e+09  \n",
       "132  1.605998e+09  \n",
       "133  1.605999e+09  \n",
       "134  1.605999e+09  \n",
       "135  1.605999e+09  \n",
       "\n",
       "[136 rows x 7 columns]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Print the dataframe\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df = df.iloc[100:102]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integration, bg-spline and peak-index creation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section takes all the images in the above table and integrates them according to the calibrations defined previously for all 4 detectors. \n",
    "This is saved as a .dat file with a prefix containing the motor positions. Then, the files in the template directory are copied with the same prefix. \n",
    "A background spline is created from the baseline points specified above and saved with the .bg-spline.dat suffix.\n",
    "The a peak-index.dat file is made based on the Zr indexes specifed above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "z from 14.782 to 14.858\n",
      "y from -5.500 to -4.000\n",
      "Current:   136 / 136\ty = -4.000 / -4.000\tz = 14.858 / 14.858     [dpa05_00474.cbf]       1 file error(s)       "
     ]
    }
   ],
   "source": [
    "# Get list of unique motor position values\n",
    "z_values = df['idtz2'].unique()\n",
    "y_values = df['idty2(encoder)'].unique()\n",
    "\n",
    "if not os.path.exists(outputdir):\n",
    "    os.makedirs(outputdir)\n",
    "if not os.path.exists(outputdir + '/0plots'):\n",
    "    os.makedirs(outputdir + '/0plots')\n",
    "    \n",
    "pe1_darks_list = glob.glob(darks + '/*.cbf')[2:]\n",
    "dark = load_tifs(pe1_darks_list)\n",
    "\n",
    "print('z from {0:.3f} to {1:.3f}'.format(np.min(z_values), np.max(z_values)))\n",
    "print('y from {0:.3f} to {1:.3f}'.format(np.min(y_values), np.max(y_values)))\n",
    "\n",
    "int_list = []\n",
    "num_error = 0\n",
    "\n",
    "for index, row in df.iterrows():\n",
    "    y=row['idty2(encoder)']\n",
    "    z=row['idtz2']\n",
    "    \n",
    "    prefix = 'y_{1:.3f}_z_{0:.3f}'.format(z,y)\n",
    "\n",
    "    file = (df[(df['idtz2']==z) & (df['idty2(encoder)']==y)])['filename'].values[0]\n",
    "    print('\\rCurrent:   {0} / {1}\\ty = {2:.3f} / {3:.3f}\\tz = {4:.3f} / {5:.3f}     [{6}]       {7} file error(s)       '\n",
    "          .format(index+1, len(df.index), y, np.max(y_values), z, np.max(z_values), file, num_error), end='')\n",
    "\n",
    "    file = directory + file\n",
    "\n",
    "    ################### Integrate #############################################\n",
    "\n",
    "    try:\n",
    "        data = load_tifs(file) - dark\n",
    "    except:\n",
    "        num_error +=1\n",
    "    else:\n",
    "        outputs = []\n",
    "        if azimuth_range == None:\n",
    "            output = ai.integrate1d(data, npt=intpoints, correctSolidAngle=True, method=method, \n",
    "                                    polarization_factor=0.99, unit='2th_deg', radial_range=(limits[0], limits[1]))\n",
    "            xvals = output[0]\n",
    "            yvals = output[1]\n",
    "            \n",
    "        else:\n",
    "            ## If azimuth range is specified\n",
    "            for azim in azimuth_range:\n",
    "                outputs.append(ai.integrate1d(data, npt=intpoints, correctSolidAngle=True, method=method, azimuth_range=azim,\n",
    "                                    polarization_factor=0.99, unit='2th_deg', radial_range=(limits[0], limits[1])))\n",
    "\n",
    "            xvals=outputs[0][0]; \n",
    "            yvals=np.sum([output[1] for output in outputs], axis=0); \n",
    "        \n",
    "        yvals = (yvals - np.min(yvals)) + 20\n",
    "\n",
    "        for templateName in glob.glob(templates + '*'):\n",
    "            copyfile(templateName, outputdir + prefix + templateName.split('/')[-1][8:])\n",
    "\n",
    "        ####################### Save integrated data ######################\n",
    "\n",
    "        with open(outputdir + prefix + '.dat', 'w+') as f:\n",
    "             np.savetxt(fname = f, X=np.transpose([xvals, yvals]), fmt = ('%1.5f'))\n",
    "\n",
    "        ########################### Make figure ###########################\n",
    "\n",
    "        plt.ioff()\n",
    "        fig, (ax2) = plt.subplots(1, 1, figsize=(16,8))\n",
    "        ax2.set_title('Integrated data');\n",
    "        ax2.set_xlabel('2theta (deg)'); \n",
    "        ax2.set_ylabel('Intensity');\n",
    "        ax2.plot(xvals, yvals)\n",
    "        ax2.set_xlim(np.min(xvals)+0.01, np.max(xvals)-0.01)\n",
    "\n",
    "        ########################### Make bg-spline ###########################\n",
    "      \n",
    "        baseline_pos, baseline_int, cs = getBaseline(xvals, yvals, baseline, baseline_interpolate, baseline_interpolate_factor,\n",
    "                                             write_to = outputdir + prefix + '.bg-spline.dat')\n",
    "\n",
    "        ax2.plot(xvals, cs(xvals))\n",
    "        ax2.plot(baseline_pos, baseline_int, 'o',c='r')\n",
    "\n",
    "        ########################### Make Zr peak-index ###########################\n",
    "\n",
    "        peak_pos_new, peak_name_new, peak_int_new = getPeaks(xvals, yvals, \n",
    "                                             peak_pos, peak_name, cs, searchrange, \n",
    "                                                    ax=ax2, plotcolour='r',\n",
    "                                             write_to = outputdir + prefix + '.peak-index.dat')\n",
    "        \n",
    "        ########################\n",
    "        if spp is True:\n",
    "\n",
    "            peak_pos_new, peak_name_new, peak_int_new = getPeaks(xvals, yvals, \n",
    "                                                 spp_peak_pos, spp_peak_name, cs, searchrange, \n",
    "                                                    mode='a', phasenumber=2, \n",
    "                                                    ax=ax2, plotcolour='g',\n",
    "                                                 write_to = outputdir + prefix + '.peak-index.dat')\n",
    "            peak_pos_new, peak_name_new, peak_int_new = getPeaks(xvals, yvals, \n",
    "                                                 np.array([4.1956]), np.array(['110']), cs, 1, \n",
    "                                                    mode='a', phasenumber=2, \n",
    "                                                ax=ax2, plotcolour='g',\n",
    "                                                 write_to = outputdir + prefix + '.peak-index.dat')\n",
    "        if hyd is True:\n",
    "            peak_pos_new, peak_name_new, peak_int_new = getPeaks(xvals, yvals, \n",
    "                                                 delta_peak_pos, delta_peak_name, cs, 3, \n",
    "                                                    mode='a', phasenumber=3, \n",
    "                                                ax=ax2, plotcolour='k',\n",
    "                                                 write_to = outputdir + prefix + '.peak-index.dat')\n",
    "        ########################\n",
    "        \n",
    "    #Append to list\n",
    "    int_list.append(np.sum(yvals - cs(xvals)))\n",
    "\n",
    "    ## Save plots #################################\n",
    "\n",
    "    plt.savefig(outputdir + '/0plots/plot_' + prefix + '.pdf')\n",
    "    plt.yscale('log')\n",
    "    plt.savefig(outputdir + '/0plots/log_' + prefix + '.pdf')\n",
    "    plt.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot integrated intensities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close()\n",
    "\n",
    "for val in df['idty2(encoder)'].unique():\n",
    "    yval = np.array(int_list)[np.array(df['idty2(encoder)'].tolist() == val)]\n",
    "    zval = np.array(df['idtz2'].tolist())[np.array(df['idty2(encoder)'].tolist() == val)]\n",
    "    \n",
    "    # Take differential and calculate edge position\n",
    "    edge = zval[np.argmax(np.gradient(yval))]\n",
    "    \n",
    "    plt.plot(zval, yval, label = 'y = {0:.3f}   edge = {1:.3f}'.format(val, edge))\n",
    "    \n",
    "plt.xlabel('Z position (mm)')\n",
    "plt.ylabel('Integrated intensity')\n",
    "plt.show()\n",
    "plt.legend(loc='lower right')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.savefig(outputdir + '/0plots/0integrated_intensity.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make a bash script"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save this text into a .sh file and run it - this will execute CMWP for each file sequentially"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmwpfolder = \"/home/rhys/Documents/CMWP-211102/\"\n",
    "\n",
    "for index, row in df.iterrows():\n",
    "    y=row['idty2(encoder)']; z=row['idtz2'];    \n",
    "    \n",
    "    print('./evaluate ' + outputdir.split(cmwpfolder)[-1] + 'y_{1:.3f}_z_{0:.3f}.dat auto'.format(z,y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2D integration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = [load_tifs(file) for file in files]\n",
    "\n",
    "output2d = ais.integrate2d(data, npt_rad=intpoints, npt_azim=3600, correctSolidAngle=True, polarization_factor=0.99)\n",
    "output1d = ais.integrate1d(data, npt=intpoints, correctSolidAngle=True, polarization_factor=0.99)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(2,1, figsize=(14,12))\n",
    "ax1.imshow(np.where(output2d[0]==0, np.nan, output2d[0]), vmin=0, vmax=0.5e11)\n",
    "\n",
    "ax2.plot(output1d[0], output1d[1])\n",
    "ax2.set_xlim(np.min(output1d[0]),np.max(output1d[0]))\n",
    "ax2.set_yscale('log')\n",
    "\n",
    "for name, pos in zip(peak_name, peak_pos):\n",
    "    \n",
    "    if np.min(output1d[0]) < pos < np.max(output1d[0]):\n",
    "\n",
    "        frac = 1.42*intpoints*(pos-np.min(output1d[0]))/(np.max(output1d[0]))\n",
    "        \n",
    "        # draw line and print name\n",
    "        ax2.axvline(pos, alpha=0.3, c='r')\n",
    "        ax1.axvline(frac, alpha=0.3, c='r')\n",
    "        \n",
    "        ax1.text(frac, 3500, name, horizontalalignment = 'center', c='r')\n",
    "        ax2.text(pos, np.max(output1d[1]), name, horizontalalignment = 'center', c='r')\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}