Zr_DESY_2021.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DESY depth profiling data analysis workbook [2021]"
   ]
  },
  {
   "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": null,
   "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 globs\n",
    "import pyFAI, pyFAI.azimuthalIntegrator\n",
    "from pyFAI.multi_geometry import MultiGeometry\n",
    "import csv\n",
    "from scipy.interpolate import CubicSpline\n",
    "from shutil import copyfile\n",
    "from src.cmwp_tools import load_tifs\n",
    "from src.xrd_tools import getReflections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Experiment specific settings\n",
    "\n",
    "base = \"/mnt/manchester_rds/202107_DESY_DP/\"                            # This is the base directory where the data is stored\n",
    "\n",
    "calib_files = [base + 'det1_calib.poni',\n",
    "               base + 'det2_calib.poni',\n",
    "               base + 'det3_calib.poni',\n",
    "               base + 'det4_calib.poni']                                # These are the calibration files (.poni) for the detectors\n",
    "\n",
    "mask_files = [base + 'det1_mask.edf',\n",
    "               base + 'det2_mask.edf',\n",
    "               base + 'det3_mask.edf',\n",
    "               base + 'det4_mask.edf']                                  # These are the mask files for each detector\n",
    "\n",
    "templates = base + \"templates/\"                                         # This is where the template files are stored. These are copied for each integration.\n",
    "\n",
    "intpoints = 10000                                                       # This is how many integration points to use for pyFAI\n",
    "\n",
    "wavelength = 0.1836803                                                  # Wavelength in angstrom\n",
    "a=3.232                                                                 # a lattice paramater in angstrom\n",
    "c=5.147                                                                 # c lattice paramater in angstrom\n",
    "\n",
    "baseline=[3.45, 4.6, 5.1, 6.2, 6.8, 8.3, 9.3, 10.32, 11.45]             # What 2theta points to calculte the background spline from\n",
    "\n",
    "limits=[3.4,11.5]                                                       # 2theta bounds of integration\n",
    "\n",
    "searchrange = 100                                                       # How many data points each side of the approximate 2theta peak position to search for the true peak\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Sample specific settings\n",
    "               \n",
    "directory = base + \"depth_prof/sam1a\"                                   # Directory containing the data (should have folders 1, 2, 3, 4 and a .fio file within)\n",
    "\n",
    "outputdir = \"/home/rhys/Documents/CMWP-210315/DESY_202107/sam1a/\"       # This is the output directory (normally a sub-folder in your CMWP dir)"
   ]
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get .fio file path\n",
    "names = glob.glob(directory+'/*.fio')\n",
    "if len(names) == 1:\n",
    "    fio_file = names[0]\n",
    "else:\n",
    "    raise Exception('Either there are no .fio files in the directory, or there are more than 1')\n",
    "    \n",
    "# Read in column names and data type\n",
    "colnames = []; formats = [];\n",
    "with open(fio_file) as input_data:\n",
    "    for i, line in enumerate(input_data):\n",
    "        if ' Col' in line:\n",
    "            colnames.append(' '.join(line.split(' ')[3:-1]))\n",
    "            skip = i+1\n",
    "            if 'DOUBLE' in line.split(' ')[-1]: formats.append('f4')\n",
    "            if 'INTEGER' in line.split(' ')[-1]: formats.append('i4')\n",
    "            if 'STRING' in line.split(' ')[-1]: formats.append('str')\n",
    "\n",
    "# Read in log file into dataframe and remove clearing frames\n",
    "df = pd.read_csv(fio_file, names = colnames, skiprows=skip, sep=' ', skipinitialspace=True) \n",
    "df = df[df.type != 'clearing']\n",
    "\n",
    "# Get list of unique motor position values\n",
    "z_values = df['idtz2'].unique()\n",
    "y_values = df['idty1(encoder)'].unique()\n",
    "    \n",
    "df.reset_index(inplace=True)\n",
    "df.drop('index', axis=1, inplace=True)\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "peak_name, peak_pos = getReflections(crystalType='hcp', a=a, c=c, wavelength=wavelength, printReflections=False)\n",
    "\n",
    "ais = [pyFAI.load(calib_file) for calib_file in calib_files]\n",
    "mask = [load_tifs(mask_file) for mask_file in mask_files]\n",
    "ais = MultiGeometry(ais, unit='2th_deg', radial_range=(limits[0], limits[1]))\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",
    "for index, row in df.iterrows():\n",
    "    y=row['idty1(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['idty1(encoder)']==y)])['filename'].values[0]\n",
    "    print('\\rCurrent:   {0} / {1}\\ty = {2:.3f} / {3:.3f}\\tz = {4:.3f} / {5:.3f}     [{6}] '\n",
    "          .format(index+1, len(df.index), y, np.max(y_values), z, np.max(z_values), file), end='')\n",
    "\n",
    "    files = [directory + '/1/' + file, directory + '/2/' + file, directory + '/3/' + file, directory + '/4/' + file]\n",
    "    \n",
    "    data = [load_tifs(file) for file in files]\n",
    "\n",
    "    output = ais.integrate1d(data, npt=intpoints, lst_mask=mask, correctSolidAngle=True, polarization_factor=0.99)\n",
    "\n",
    "    xvals=output[0]; yvals=output[1];\n",
    "    yvals = yvals / 100000\n",
    "    #yvals = yvals - np.min(yvals)\n",
    "    #yvals += 10000\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",
    "    ########################### Save 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",
    "    #ax2.set_ylim(np.min(yvals)-2, np.max(yvals)*1.1)\n",
    "\n",
    "    x_plot_list = []\n",
    "    y_plot_list = []\n",
    "\n",
    "    ########################### Make bg-spline ###########################\n",
    "    baseline_int = []\n",
    "    for j in baseline:\n",
    "        num_index=np.argmin(np.abs(xvals-j))\n",
    "\n",
    "        baseline_int.append(np.mean(yvals[num_index-5:num_index+5]))\n",
    "\n",
    "    baseline, baseline_int = (list(t) for t in zip(*sorted(zip(baseline, baseline_int))))\n",
    "\n",
    "    with open(outputdir + prefix + '.bg-spline.dat', 'w+') as f:\n",
    "        np.savetxt(fname = f, X=np.transpose([baseline, baseline_int]), fmt = ('%1.5f'))\n",
    "        \n",
    "    cs = CubicSpline(baseline, baseline_int)\n",
    "\n",
    "    ax2.plot(xvals, cs(xvals))\n",
    "    ax2.plot(baseline, baseline_int, 'o',c='r')\n",
    "\n",
    "    ########################### Make peak-index ###########################\n",
    "\n",
    "    if len(peak_pos) != len(peak_name):\n",
    "        raise ValueError('peak_pos and peak_name arrays should be the same size')\n",
    "\n",
    "    with open(outputdir + prefix + '.peak-index.dat', 'w', newline='') as csvfile:\n",
    "        writer = csv.writer(csvfile)\n",
    "        \n",
    "        for name, pos in zip(peak_name, peak_pos):\n",
    "\n",
    "            if np.min(xvals) < pos < np.max(xvals):         # if the peak is within the data\n",
    "                \n",
    "                # approx peak postion\n",
    "                approx_peak_index = np.argmin(np.abs(xvals-pos))                    # get the index of the peak\n",
    "\n",
    "                # get actual peak position\n",
    "                peak_index = np.argmax(yvals[approx_peak_index-searchrange:approx_peak_index+searchrange])+approx_peak_index-searchrange\n",
    "                yval = yvals[peak_index]\n",
    "\n",
    "                # draw line and print name\n",
    "                ax2.axvline(xvals[peak_index], alpha=0.1, c='r')\n",
    "                ax2.text(xvals[peak_index], yval+10, name, horizontalalignment = 'center', c='r')\n",
    "\n",
    "                intensity = yval - cs(xvals)[peak_index]\n",
    "\n",
    "                writer.writerow(['{0:.4f} {1:.1f} {2} 0'.format(xvals[peak_index], intensity, name)])\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": "markdown",
   "metadata": {},
   "source": [
    "## Plot integrated intensities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for val in df['idty1(encoder)'].unique():\n",
    "    yval = np.array(int_list)[np.array(df['idty1(encoder)'].tolist() == val)]\n",
    "    zval = np.array(df['idtz2'].tolist())[np.array(df['idty1(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 + '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-210315/\"\n",
    "\n",
    "for index, row in df.iterrows():\n",
    "    y=row['idty1(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
}