Zr_DESY_2019.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": [],
   "source": [
    "\n",
    "#name = '02_0.15dpa_Zr4'\n",
    "#dpa_val = 0.15\n",
    "\n",
    "#lookup_val = 'y1'\n",
    "#zlookup = [[-0.4, -0.2, 0, 0.2, 0.4],[13.635, 13.633, 13.632, 13.631, 13.630]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wavelength: 0.32541784 A\n"
     ]
    }
   ],
   "source": [
    "### Experiment specific settings\n",
    "\n",
    "base = \"/mnt/manchester_rds/201910_DESY/\"                            # This is the base directory where the data is stored\n",
    "\n",
    "calib_file = [base + 'processed/25Oct_PE1_calib.poni', base + 'processed/25Oct_PE2_calib.poni']             # These are the calibration files (.poni) for the detectors\n",
    "\n",
    "ais = [pyFAI.load(calib) for calib in calib_file]\n",
    "wavelength = ais[0].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",
    "## Background\n",
    "baseline= [5.734,5.929,6.112,6.289,6.42, 6.973,7.962,8.164,9.246, 10.378,10.7, 11, 11.912, \n",
    "        12.21,12.42,13.15,14.107,14.2182,14.783,14.968,15.588,16.236,\n",
    "        16.595,16.927,17.5007,18.329,18.9106,19.65, 17.88, 8.5, 8.8]            # What 2theta points to calculte the background spline from\n",
    "\n",
    "baseline_interpolate=[]                            # What 2theta points to add in extra point\n",
    "baseline_interpolate_factor=[]\n",
    "\n",
    "limits=[5.5, 20]                                                       # 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": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Sample specific settings\n",
    "\n",
    "name = '0.1 dpa'\n",
    "\n",
    "directories = [\"00_BESTDATA/02_Zr4_0p1_realScan2_top_loop1\", \"00_BESTDATA/02_Zr4_0p1_realScan2_top_loop2\",\n",
    "    \"00_BESTDATA/02_Zr4_0p1_realScan2_top_loop3\", \"00_BESTDATA/02_Zr4_0p1_realScan2_bottom_loop1\",\n",
    "    \"00_BESTDATA/02_Zr4_0p1_realScan2_bottom_loop2\"]\n",
    "\n",
    "outputdir = \"/home/rhys/Documents/CMWP-211102/2019_DESY/0p1/\"       # This is the output directory (normally a sub-folder in your CMWP dir)\n",
    "\n",
    "\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 = None\n",
    "y_range = None\n",
    "spp=False; hyd=False;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Integrate lab6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "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": 31,
   "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": 35,
   "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": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "name: 0.1 dpa\tz from 13.620 to 13.700\ty from -0.400 to 0.400\n"
     ]
    }
   ],
   "source": [
    "# Initialise lists\n",
    "idtz2_list = []; idty1_list = []; file_list = []; channel_list = [];\n",
    "\n",
    "# Get darks\n",
    "pe1_darks_list = glob.glob(base + 'raw/' + directories[0] + '/_PE1/dark*.tif')\n",
    "pe2_darks_list = glob.glob(base + 'raw/' + directories[0] + '/_PE2/dark*.tif')\n",
    "\n",
    "pe1_dark_image = load_tifs(pe1_darks_list)\n",
    "pe2_dark_image = load_tifs(pe2_darks_list)\n",
    "\n",
    "# Loop over directories\n",
    "for directory in directories:\n",
    "\n",
    "    # Read in log file\n",
    "    idtz2, curr = np.loadtxt(fname = base + 'raw/' + directory + \"/log.log\", delimiter=\" \", \n",
    "                             skiprows=7, usecols = (0, 4), unpack=True)\n",
    "    fio_filename2 = np.loadtxt(fname = base + 'raw/'+ directory + \"/log.log\", delimiter=\" \", \n",
    "                               skiprows=7, usecols = (5), dtype=str, unpack=True)\n",
    "\n",
    "    fio_filename = []\n",
    "    for k in range(len(fio_filename2)):\n",
    "        fio_filename.append(base + fio_filename2[k][-25:])\n",
    "\n",
    "    # Read in each fio file and write to lists\n",
    "    # Loop over idtz2\n",
    "    for i in range(len(idtz2)):\n",
    "        idty1,channel = np.loadtxt(fname=fio_filename[i], skiprows=29, usecols = (1, 2), unpack=True)\n",
    "        data_filename, clearing = np.loadtxt(fname=fio_filename[i], skiprows=29, usecols = (3, 4), \n",
    "                                             dtype=str, unpack=True)\n",
    "\n",
    "        # Loop over idty1\n",
    "        for j in range(len(idty1)):\n",
    "            if clearing[j] != \"clearing\":\n",
    "                    idtz2_list.append(idtz2[i])\n",
    "                    idty1_list.append(idty1[j])\n",
    "\n",
    "                    # Correct filenames (pad number with zeroes to 5 digits) and add full path\n",
    "                    if channel[j]==1:\n",
    "                        file_list.append(base + 'raw/'+directory+\"/_PE1/sample_\" + data_filename[j]\n",
    "                                         .split(\"_\")[1].zfill(9))\n",
    "                    if channel[j]==2:\n",
    "                        file_list.append(base + 'raw/'+directory+\"/_PE2/sample_\" + data_filename[j]\n",
    "                                         .split(\"_\")[1].zfill(9))\n",
    "                    channel_list.append(channel[j])\n",
    "\n",
    "# Convert to pandas dataframe\n",
    "df = pd.DataFrame(data={'channel': channel_list, 'idtz2': idtz2_list, 'idty1': idty1_list, 'filename': file_list})\n",
    "\n",
    "z_values = df['idtz2'].unique()\n",
    "y_values = df['idty1'].unique()\n",
    "\n",
    "pref = '\\n';\n",
    "print(pref+'name: {0}\\tz from {1:.3f} to {2:.3f}\\ty from {3:.3f} to {4:.3f}'.format(name, np.min(z_values), np.max(z_values), np.min(y_values), np.max(y_values)))\n",
    "\n",
    "# Group image names by channel \\ idtz2 \\ idty1\n",
    "new = df.groupby(by=['idtz2', 'idty1', 'channel'])['filename'].apply(list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "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>channel</th>\n",
       "      <th>idtz2</th>\n",
       "      <th>idty1</th>\n",
       "      <th>filename</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.0</td>\n",
       "      <td>13.62</td>\n",
       "      <td>-0.4</td>\n",
       "      <td>/mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>13.62</td>\n",
       "      <td>-0.4</td>\n",
       "      <td>/mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0</td>\n",
       "      <td>13.62</td>\n",
       "      <td>-0.2</td>\n",
       "      <td>/mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.0</td>\n",
       "      <td>13.62</td>\n",
       "      <td>-0.2</td>\n",
       "      <td>/mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1.0</td>\n",
       "      <td>13.62</td>\n",
       "      <td>0.0</td>\n",
       "      <td>/mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2855</th>\n",
       "      <td>2.0</td>\n",
       "      <td>13.70</td>\n",
       "      <td>0.2</td>\n",
       "      <td>/mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2856</th>\n",
       "      <td>2.0</td>\n",
       "      <td>13.70</td>\n",
       "      <td>0.4</td>\n",
       "      <td>/mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2857</th>\n",
       "      <td>2.0</td>\n",
       "      <td>13.70</td>\n",
       "      <td>0.4</td>\n",
       "      <td>/mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2858</th>\n",
       "      <td>2.0</td>\n",
       "      <td>13.70</td>\n",
       "      <td>0.4</td>\n",
       "      <td>/mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2859</th>\n",
       "      <td>2.0</td>\n",
       "      <td>13.70</td>\n",
       "      <td>0.4</td>\n",
       "      <td>/mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2860 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      channel  idtz2  idty1                                           filename\n",
       "0         1.0  13.62   -0.4  /mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...\n",
       "1         1.0  13.62   -0.4  /mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...\n",
       "2         1.0  13.62   -0.2  /mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...\n",
       "3         1.0  13.62   -0.2  /mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...\n",
       "4         1.0  13.62    0.0  /mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...\n",
       "...       ...    ...    ...                                                ...\n",
       "2855      2.0  13.70    0.2  /mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...\n",
       "2856      2.0  13.70    0.4  /mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...\n",
       "2857      2.0  13.70    0.4  /mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...\n",
       "2858      2.0  13.70    0.4  /mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...\n",
       "2859      2.0  13.70    0.4  /mnt/manchester_rds/201910_DESY/raw/00_BESTDAT...\n",
       "\n",
       "[2860 rows x 4 columns]"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Print the dataframe\n",
    "df"
   ]
  },
  {
   "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": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "z from 14.917 to 14.975\n",
      "y from 20.695 to 21.695\n",
      "Current:   118 / 118\ty = 20.695 / 21.695\tz = 14.975 / 14.975     [hyd_00533.cbf]       0 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 hyd 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 spp 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": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot integrated intensities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "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')\n",
    "\n",
    "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": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.917.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.917.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.918.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.918.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.919.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.919.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.920.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.920.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.921.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.921.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.922.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.922.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.923.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.923.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.924.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.924.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.925.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.925.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.926.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.926.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.927.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.927.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.928.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.928.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.929.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.929.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.930.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.930.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.931.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.931.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.932.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.932.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.933.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.933.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.934.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.934.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.935.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.935.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.936.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.936.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.937.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.937.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.938.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.938.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.939.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.939.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.940.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.940.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.941.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.941.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.942.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.942.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.943.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.943.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.944.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.944.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.945.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.945.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.946.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.946.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.947.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.947.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.948.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.948.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.949.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.949.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.950.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.950.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.951.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.951.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.952.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.952.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.953.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.953.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.954.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.954.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.955.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.955.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.956.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.956.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.957.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.957.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.958.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.958.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.959.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.959.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.960.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.960.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.961.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.961.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.962.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.962.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.963.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.963.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.964.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.964.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.965.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.965.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.966.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.966.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.967.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.967.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.968.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.968.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.969.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.969.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.970.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.970.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.971.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.971.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.972.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.972.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.973.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.973.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.974.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.695_z_14.974.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_21.195_z_14.975.dat auto\n",
      "./evaluate 2020_11_DESY/hyd_3p_80_100/y_20.695_z_14.975.dat auto\n"
     ]
    }
   ],
   "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
}