PyHyperScattering_Batch_SST1_JupyterHub.ipynb

{
 "cells": [
  {
   "cell_type": "raw",
   "id": "6cc84247",
   "metadata": {},
   "source": [
    "!pip install --upgrade pyhyperscattering"
   ]
  },
  {
   "cell_type": "raw",
   "id": "71450613",
   "metadata": {},
   "source": [
    "pip install -e git+https://github.com/phonghnguyen/PyHyperScattering.git#egg=PyHyperScattering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dafbbea4-2f41-4917-b191-aae994f4e7aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80f0526c-4e08-4c2b-8d5c-7e0290f19ba5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import PyHyperScattering\n",
    "import os\n",
    "import shutil\n",
    "import psutil\n",
    "import pickle\n",
    "import imageio\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import matplotlib as mpl\n",
    "import matplotlib.ticker as ticker\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import font_manager\n",
    "from matplotlib.colors import LogNorm\n",
    "from matplotlib.colorbar import Colorbar as colorbar\n",
    "from IPython.display import clear_output\n",
    "from tqdm.auto import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe6b89a1-d182-4972-bc70-f84ac69d8e77",
   "metadata": {},
   "outputs": [],
   "source": [
    "# variables to store directories, enables easy switching between systems\n",
    "font_dirs = ['/nsls2/users/pnguyen2/Fonts/']\n",
    "\n",
    "base_data_dirs = ['/nsls2/data/sst/legacy/RSoXS/suitcased_data/users/2022-3/309585-UCSB/auto/Blends/2022-09-27/',\n",
    "             '/nsls2/data/sst/legacy/RSoXS/suitcased_data/users/2022-3/309585-UCSB/auto/Dopant/2022-09-27/',\n",
    "             '/nsls2/data/sst/legacy/RSoXS/suitcased_data/users/2022-3/309585-UCSB/auto/Dopant/2022-09-28/']\n",
    "\n",
    "base_save_dir = '/nsls2/users/pnguyen2/Processing Outputs/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d7f0cd6-1d41-46e9-8b1a-34ce46b02390",
   "metadata": {},
   "outputs": [],
   "source": [
    "font_files = font_manager.findSystemFonts(fontpaths = font_dirs)\n",
    "\n",
    "for font_file in font_files:\n",
    "    font_manager.fontManager.addfont(font_file)\n",
    "\n",
    "# font_names = [f.name for f in font_manager.fontManager.ttflist]\n",
    "# print(font_names)\n",
    "plt.rcParams['font.family'] = 'sans-serif'\n",
    "plt.rcParams['font.sans-serif'] = 'Avenir'\n",
    "plt.rcParams['font.size'] = 18\n",
    "plt.rcParams['axes.linewidth'] = 2\n",
    "\n",
    "cmap = mpl.cm.viridis\n",
    "cmap.set_bad((68/255, 1/255, 84/255), 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39a08fdc-8f38-46f3-902f-b52283c3f934",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_immediate_subdirectories(a_dir):\n",
    "    \n",
    "    return [name for name in os.listdir(a_dir)\n",
    "            if os.path.isdir(os.path.join(a_dir, name))]\n",
    "\n",
    "def fmt(x, pos):\n",
    "    a, b = '{:.0e}'.format(x).split('e')\n",
    "    b = int(b)\n",
    "    \n",
    "    return r'${} \\times 10^{{{}}}$'.format(a, b)\n",
    "\n",
    "def find_edge(scan):\n",
    "    if np.min(scan.energy) == 270:\n",
    "        current_edge = 'C 1s'\n",
    "        current_edge_folder = 'C1s'\n",
    "    elif np.min(scan.energy) == 385:\n",
    "        current_edge = 'N 1s'\n",
    "        current_edge_folder = 'N1s'\n",
    "    elif np.min(scan.energy) == 670:\n",
    "        current_edge = 'F 1s'\n",
    "        current_edge_folder = 'F1s'\n",
    "    else:\n",
    "        current_edge = 'unknown'\n",
    "        current_edge_folder = 'unknown'\n",
    "        \n",
    "    return [current_edge, current_edge_folder]\n",
    "\n",
    "def match_pols(scans):\n",
    "    '''\n",
    "    Iterates through list of scans in scans (made by appending scans to an initally empty list). \n",
    "    First scan in list is used as a reference and the corresponding scan of orthogonal polarization is found.\n",
    "    As matches are found, items from the list are removed and returned, along with the scans that were matched.\n",
    "    The intention is that the list is fed back to this function to find the next pair to process until the list is empty.\n",
    "    '''\n",
    "    # Case for finding 0 pol at list beginning, find matching 90 pol\n",
    "    if scans[0].unstack('system').polarization == 0:\n",
    "        current_scan_0pol = scans[0]\n",
    "        del scans[0]\n",
    "        [current_edge_0pol, current_edge_folder_0pol] = find_edge(current_scan_0pol)\n",
    "        print(f'Scan {current_scan_0pol.sampleid} matches {current_scan_0pol.sample_name} {current_edge_0pol} ({current_scan_0pol.rsoxs_config}) with 0°  polarity')\n",
    "\n",
    "        if len(scans) != 0:\n",
    "            for i, scan in enumerate(scans):\n",
    "                if scan.rsoxs_config == current_scan_0pol.rsoxs_config:\n",
    "                    [current_edge_90pol, current_edge_folder_90pol] = find_edge(scan)\n",
    "                    if current_edge_90pol == current_edge_0pol:\n",
    "                        if scan.sample_name == current_scan_0pol.sample_name:\n",
    "                            if scan.unstack('system').polarization == 90:\n",
    "                                current_scan_90pol = scan\n",
    "                                print(f'Scan {scan.sampleid} matches {scan.sample_name} {current_edge_90pol} ({current_scan_90pol.rsoxs_config}) with 90° polarity')\n",
    "                                del scans[i]\n",
    "                                return [current_scan_0pol, current_scan_90pol, scans]\n",
    "            current_scan_90pol = []\n",
    "        else:\n",
    "            print('Reached end, unable to pair last scan')\n",
    "            \n",
    "    # Case for finding 90 pol at list beginning, find matching 0 pol\n",
    "    if scans[0].unstack('system').polarization == 90:\n",
    "        current_scan_90pol = scans[0]\n",
    "        del scans[0]\n",
    "        \n",
    "        [current_edge_90pol, current_edge_folder_90pol] = find_edge(current_scan_90pol)\n",
    "            \n",
    "        print(f'Scan {current_scan_90pol.sampleid} matches {current_scan_90pol.sample_name} {current_edge_90pol} ({current_scan_90pol.rsoxs_config}) with 90° polarity')\n",
    "\n",
    "        if len(scans) != 0:\n",
    "            for i, scan in enumerate(scans):\n",
    "                if scan.rsoxs_config == current_scan_90pol.rsoxs_config:\n",
    "                    [current_edge_0pol, current_edge_folder_0pol] = find_edge(scan)\n",
    "                    if current_edge_0pol == current_edge_90pol:   \n",
    "                        if scan.sample_name == current_scan_90pol.sample_name:\n",
    "                            if scan.unstack('system').polarization == 0:\n",
    "                                current_scan_0pol = scan\n",
    "                                [current_edge, current_edge_folder] = find_edge(current_scan_0pol)\n",
    "                                print(f'Scan {scan.sampleid} matches {scan.sample_name} {current_edge_0pol} ({current_scan_0pol.rsoxs_config}) with 0°  polarity')\n",
    "                                del scans[i]\n",
    "                                return [current_scan_0pol, current_scan_90pol, scans]\n",
    "            current_scan_0pol = []\n",
    "        else:\n",
    "            print('Reached end, unable to pair last scan')\n",
    "            \n",
    "    return [current_scan_0pol, current_scan_90pol, scans]\n",
    "\n",
    "def assign_mask(scan):\n",
    "    if scan.rsoxs_config == 'waxs':\n",
    "        mask_method = 'nika'\n",
    "        mask = 'SST1_Sept2022_WAXS_mask.hdf'\n",
    "    elif scan.rsoxs_config == 'saxs':\n",
    "        mask_method = 'none'\n",
    "        mask = []\n",
    "    else:\n",
    "        mask_method = 'none'\n",
    "        mask = []\n",
    "    \n",
    "    return [mask_method, mask]\n",
    "\n",
    "def make_output_dir(base_save_dir, scan, tail):\n",
    "# tail is a string that describes the data file being exported, used as folder prefix\n",
    "    [current_edge, current_edge_folder] = find_edge(scan)\n",
    "    \n",
    "    # create directories for output files\n",
    "    save_dir = os.path.join(base_save_dir, f'{tail}_{sample_name_root(scan)}')\n",
    "    os.makedirs(save_dir, exist_ok = True)\n",
    "    \n",
    "    return save_dir\n",
    "\n",
    "def clear_dir(dir):\n",
    "    for file_name in os.listdir(dir):\n",
    "        file_path = os.path.join(dir, file_name)\n",
    "        try:\n",
    "            if os.path.isfile(file_path) or os.path.islink(file_path):\n",
    "                os.unlink(file_path)\n",
    "            elif os.path.isdir(file_path):\n",
    "                shutil.rmtree(file_path)\n",
    "        except Exception as e:\n",
    "            print('Failed to delete %s. Reason: %s' % (file_path, e))\n",
    "            \n",
    "def sample_name_root(scan):\n",
    "    [current_edge, current_edge_folder] = find_edge(scan)\n",
    "    return f'{current_edge_folder}_{scan.sample_name}'\n",
    "\n",
    "def scan_loader_SST1(base_data_dirs, exclude = [], ll = [], ul = []):\n",
    "    '''\n",
    "    Args:\n",
    "        base_data_dirs (list of str): path to the directories to load\n",
    "        exlude (list of str): list of scan IDs to exclude from loading into scans\n",
    "        ll (str): scan ID below which to exclude from loading into scans\n",
    "        ul (str) scan ID above which to exclude from loading into scans\n",
    "    '''\n",
    "    base_data_dir_tracker = 0\n",
    "    scan_tracker = 0\n",
    "    scan_skip_corr = 0\n",
    "    scans = []\n",
    "    for base_data_dir in base_data_dirs:\n",
    "        base_data_dir_tracker += 1\n",
    "        scan_tracker = 0\n",
    "        scan_skip_corr = 0\n",
    "        all_scans = pd.Series(get_immediate_subdirectories(base_data_dir))\n",
    "        for current_scan in np.sort(all_scans):\n",
    "            scan_tracker += 1\n",
    "            if ((len(ll) == 0) and (len(ul) == 0)):\n",
    "                if current_scan in exclude:\n",
    "                    print(f'Skipping scan {current_scan} - ID in list to exclude')\n",
    "                    scan_skip_corr += 1\n",
    "                else:\n",
    "                    clear_output(wait = True)\n",
    "                    print(f'RAM usage: {psutil.virtual_memory().available * 100 / psutil.virtual_memory().total :.2f} %')\n",
    "                    print(f'CPU usage: {psutil.cpu_percent() :.2f} %')\n",
    "                    print(f'Currently pulling data from \\'{base_data_dir}\\' ({base_data_dir_tracker}/{len(base_data_dirs)})')\n",
    "                    print(f'Currently loading scan {current_scan} ({scan_tracker - scan_skip_corr}/{len(all_scans) - scan_skip_corr})')\n",
    "                    scan = file_loader.loadFileSeries(base_data_dir + current_scan, ['energy', 'polarization'])\n",
    "                    scans.append(scan)\n",
    "            elif ((len(ll) > 0) and (len(ul) == 0)):\n",
    "                if current_scan >= ll:\n",
    "                    if current_scan in exclude:\n",
    "                        print(f'Skipping scan {current_scan} - ID in list to exclude')\n",
    "                        scan_skip_corr += 1\n",
    "                    else:\n",
    "                        clear_output(wait = True)\n",
    "                        print(f'RAM usage: {psutil.virtual_memory().available * 100 / psutil.virtual_memory().total :.2f} %')\n",
    "                        print(f'CPU usage: {psutil.cpu_percent() :.2f} %')\n",
    "                        print(f'Currently pulling data from \\'{base_data_dir}\\' ({base_data_dir_tracker}/{len(base_data_dirs)})')\n",
    "                        print(f'Currently loading scan {current_scan} ({scan_tracker - scan_skip_corr}/{len(all_scans) - scan_skip_corr})')\n",
    "                        scan = file_loader.loadFileSeries(base_data_dir + current_scan, ['energy', 'polarization'])\n",
    "                        scans.append(scan)\n",
    "            elif ((len(ll) == 0) and (len(ul) > 0)):\n",
    "                if current_scan <= ul:\n",
    "                    if current_scan in exclude:\n",
    "                        print(f'Skipping scan {current_scan} - ID in list to exclude')\n",
    "                        scan_skip_corr += 1\n",
    "                    else:\n",
    "                        clear_output(wait = True)\n",
    "                        print(f'RAM usage: {psutil.virtual_memory().available * 100 / psutil.virtual_memory().total :.2f} %')\n",
    "                        print(f'CPU usage: {psutil.cpu_percent() :.2f} %')\n",
    "                        print(f'Currently pulling data from \\'{base_data_dir}\\' ({base_data_dir_tracker}/{len(base_data_dirs)})')\n",
    "                        print(f'Currently loading scan {current_scan} ({scan_tracker - scan_skip_corr}/{len(all_scans) - scan_skip_corr})')\n",
    "                        scan = file_loader.loadFileSeries(base_data_dir + current_scan, ['energy', 'polarization'])\n",
    "                        scans.append(scan)\n",
    "            elif ((len(ll) > 0) and (len(ul) > 0)):\n",
    "                if ((current_scan >= ll) and (current_scan <= ul)):\n",
    "                    if current_scan in exclude:\n",
    "                        print(f'Skipping scan {current_scan} - ID in list to exclude')\n",
    "                        scan_skip_corr += 1\n",
    "                    else:\n",
    "                        clear_output(wait = True)\n",
    "                        print(f'RAM usage: {psutil.virtual_memory().available * 100 / psutil.virtual_memory().total :.2f} %')\n",
    "                        print(f'CPU usage: {psutil.cpu_percent() :.2f} %')\n",
    "                        print(f'Currently pulling data from \\'{base_data_dir}\\' ({base_data_dir_tracker}/{len(base_data_dirs)})')\n",
    "                        print(f'Currently loading scan {current_scan} ({scan_tracker - scan_skip_corr}/{len(all_scans) - scan_skip_corr})')\n",
    "                        scan = file_loader.loadFileSeries(base_data_dir + current_scan, ['energy', 'polarization'])\n",
    "                        scans.append(scan)\n",
    "\n",
    "    return scans\n",
    "\n",
    "def setup_plot():\n",
    "    plt.rcParams.update({\"font.size\": 18})\n",
    "    plt.rcParams.update({\"axes.linewidth\": 2})\n",
    "\n",
    "    fig, axes = plt.subplots(1, 1, figsize = (3.5, 3.5))\n",
    "\n",
    "    # Edit the major and minor ticks of the x and y axes\n",
    "    axes.xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', top = 'on')\n",
    "    axes.xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', top = 'on')\n",
    "    axes.yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "    axes.yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "    \n",
    "    return [fig, axes]\n",
    "    \n",
    "def save_fig(fig, scan, current_energy, save_dir, tail, res):\n",
    "    if current_energy > 0:\n",
    "        save_path = os.path.join(save_dir, f'{tail}_{sample_name_root(scan)}_{current_energy}eV.png')\n",
    "        fig.savefig(save_path, bbox_inches = 'tight', pad_inches = 0.2, dpi = res, transparent = False)\n",
    "    else:\n",
    "        save_path = os.path.join(save_dir, f'{tail}_{sample_name_root(scan)}.png')\n",
    "        fig.savefig(save_path, bbox_inches = 'tight', pad_inches = 0.2, dpi = res, transparent = False)\n",
    "        \n",
    "def make_gif(scan, input_dir, output_dir, tail):\n",
    "    image_folder = os.fsencode(input_dir)\n",
    "\n",
    "    file_names = []\n",
    "\n",
    "    for file in os.listdir(image_folder):\n",
    "        file_name = os.fsdecode(file)\n",
    "        file_names.append(file_name)\n",
    "\n",
    "    file_names.sort() # this iteration technique has no built in order, so sort the frames\n",
    "\n",
    "    # Make it pause at the end so that the viewers can ponder\n",
    "    for _ in range(10):\n",
    "        file_names.append(file_names[-1])\n",
    "\n",
    "    images = list(map(lambda file_name: imageio.imread(os.path.join(input_dir, file_name)), file_names))\n",
    "\n",
    "    imageio.mimsave(os.path.join(output_dir, f'{tail}_{sample_name_root(scan)}.gif'), images, duration = 0.5) # modify duration as needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "794f9c0d-04b9-4e1a-96fd-baf3a434872b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# sacrificial empty plot because the first plot does not properly update params\n",
    "[fig, axes] = setup_plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f1f8533-6cbc-4c81-8ba6-6b956628507e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Determined from Nika\n",
    "sample_to_det_dist = 43.11\n",
    "BC_x = 395.87\n",
    "BC_y = 552.74"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "174a5d25-afa6-4a37-adc5-238cd3daf210",
   "metadata": {},
   "outputs": [],
   "source": [
    "file_loader = PyHyperScattering.load.SST1RSoXSLoader(corr_mode = 'none')\n",
    "\n",
    "integ = PyHyperScattering.integrate.PFEnergySeriesIntegrator(geomethod = \"nika\",\n",
    "                                                             NIdistance = sample_to_det_dist, \n",
    "                                                             NIbcx = BC_x, NIbcy = BC_y,\n",
    "                                                             NIpixsizex = 0.06, NIpixsizey = 0.06,\n",
    "                                                             correctSolidAngle = True, \n",
    "                                                             integration_method = 'csr_ocl',\n",
    "                                                             npts = 200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cadff66-0205-4065-b078-360f4aaf2d8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "scans = []\n",
    "del scans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4770a221-d86a-4520-9e70-2148818fcda5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# scans = scan_loader_SST1(base_data_dirs, exclude = ['44068', '44109'], ll = '44065')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b20fc70e-1d3e-4790-9c10-52692124f38f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # save scans to pickle\n",
    "# with open('scans.pkl', 'wb') as file:\n",
    "#     pickle.dump(scans, file, protocol = -1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd4500ed-ef92-4db8-9e01-964c5da00bf5",
   "metadata": {},
   "outputs": [],
   "source": [
    "scans = []\n",
    "del scans\n",
    "\n",
    "# load scans from pickle\n",
    "with open('scans.pkl', 'rb') as file:\n",
    "    scans = pickle.load(file)\n",
    "    \n",
    "# variables to append sample ID and AR xarrays\n",
    "sample_names = []\n",
    "edges = []\n",
    "\n",
    "ARs = []\n",
    "AR_scales = []\n",
    "\n",
    "paras = []\n",
    "para_scales = []\n",
    "\n",
    "perps = []\n",
    "perp_scales = []\n",
    "\n",
    "circs = []\n",
    "circ_scales = []\n",
    "\n",
    "# check whether some progress has already been made:\n",
    "save_path = os.path.join(os.getcwd(), 'PyHyperScattering_Batch_SST1_JupyterHub.pkl')\n",
    "tail = os.path.split(save_path)[1]\n",
    "\n",
    "if tail in os.listdir(os.getcwd()):   \n",
    "    with open('PyHyperScattering_Batch_SST1_JupyterHub.pkl', 'rb') as file:\n",
    "        [sample_names, edges, ARs, AR_scales, paras, para_scales, perps, perp_scales, circs, circ_scales] = pickle.load(file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e779d71a-9637-47ea-a75a-387276044e27",
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2afd7b92-5cc1-4987-828e-e9af86d00733",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print(f'{len(scans)} scans imported')\n",
    "\n",
    "\n",
    "\n",
    "while len(scans) > 0:\n",
    "    \n",
    "    # reset energy\n",
    "    current_energy = 0\n",
    "        \n",
    "    # Find samples with matching experimental configurations, elemental edges, names, and corresponding polarizations\n",
    "    # also takes care of removing already processed scans from the list of scans remaining\n",
    "    current_scan_pol0 = []\n",
    "    current_scan_pol90 = []\n",
    "    \n",
    "    while ((len(current_scan_pol0) == 0) or (len(current_scan_pol90) == 0) and (len(scans) > 0)):\n",
    "        [current_scan_pol0, current_scan_pol90, scans] = match_pols(scans)\n",
    "    \n",
    "    # make sure final processing output directory exists\n",
    "    reduced_save_dir = make_output_dir(base_save_dir, current_scan_pol0, '0_reduced')\n",
    "    \n",
    "    # check whether sample has been fully processed\n",
    "    save_path = os.path.join(reduced_save_dir, f'abs_linecuts_{sample_name_root(current_scan_pol0)}.gif')\n",
    "    tail = os.path.split(save_path)[1]\n",
    "\n",
    "    # if tail not in os.listdir(reduced_save_dir):\n",
    "    if True:\n",
    "    \n",
    "        # Assign the mask depending on experimental config\n",
    "        [mask_method, mask] = assign_mask(current_scan_pol0)\n",
    "\n",
    "        # import mask into PyHyperScattering integrator object\n",
    "        integ.loadNikaMask(filetoload = mask, rotate_image = False)\n",
    "\n",
    "        # rotate mask to match experimental orientation\n",
    "        integ.mask = np.flipud(np.rot90(integ.mask))\n",
    "\n",
    "        # # compare experimental and mask orientation\n",
    "        # f, axarr = plt.subplots(1,2)\n",
    "        # axarr[0].imshow(scan_unstacked.sel(energy=np.max(scan.energy)),\n",
    "        #                 # norm = LogNorm(1,np.nanpercentile(scan_unstacked.sel(energy=current_energy),99)\n",
    "        #                 vmin = 0,\n",
    "        #                 vmax = np.nanpercentile(scan_unstacked.sel(energy=np.max(scan.energy)), 90),\n",
    "        #                 cmap = 'viridis',\n",
    "        #                 aspect='1')\n",
    "        # axarr[1].imshow(integ.mask, norm = LogNorm(1,np.nanpercentile(scan_unstacked.sel(energy=np.max(scan.energy)),90)))\n",
    "    \n",
    "        # Find the elemental edge of the current scan\n",
    "        [current_edge, current_edge_folder] = find_edge(current_scan_pol0)\n",
    "\n",
    "        # log info for mosaics\n",
    "        sample_names.append(current_scan_pol0.sample_name)\n",
    "        edges.append(current_edge)\n",
    "    \n",
    "        # Make sure directories exist\n",
    "        rel_detImage_save_dir = make_output_dir(base_save_dir, current_scan_pol0, 'rel_detImage')\n",
    "        abs_detImage_save_dir = make_output_dir(base_save_dir, current_scan_pol0, 'abs_detImage')\n",
    "\n",
    "        rel_linecuts_save_dir = make_output_dir(base_save_dir, current_scan_pol0, 'rel_linecuts')\n",
    "        abs_linecuts_save_dir = make_output_dir(base_save_dir, current_scan_pol0, 'abs_linecuts')\n",
    "    \n",
    "        # Remove old files\n",
    "        clear_dir(rel_detImage_save_dir)\n",
    "        clear_dir(abs_detImage_save_dir)\n",
    "\n",
    "        clear_dir(rel_linecuts_save_dir)\n",
    "        clear_dir(abs_linecuts_save_dir)\n",
    "\n",
    "        clear_dir(reduced_save_dir)\n",
    "    \n",
    "        print(f'RAM usage: {psutil.virtual_memory().available * 100 / psutil.virtual_memory().total :.2f} %')\n",
    "        print(f'CPU usage: {psutil.cpu_percent() :.2f} %')\n",
    "        \n",
    "        print(f'Integrating scan {current_scan_pol0.sampleid}:')\n",
    "              \n",
    "        integrated_data_pol0 = integ.integrateImageStack(current_scan_pol0)\n",
    "        integrated_data_pol0_unstacked = integrated_data_pol0.unstack('system')\n",
    "        \n",
    "        print(f'RAM usage: {psutil.virtual_memory().available * 100 / psutil.virtual_memory().total :.2f} %')\n",
    "        print(f'CPU usage: {psutil.cpu_percent() :.2f} %')\n",
    "        \n",
    "        print(f'Integrating scan {current_scan_pol90.sampleid}:')\n",
    "        \n",
    "        integrated_data_pol90 = integ.integrateImageStack(current_scan_pol90)\n",
    "        integrated_data_pol90_unstacked = integrated_data_pol90.unstack('system')\n",
    "        \n",
    "        # Set the energy and q range of interest depending on the elemental edge of the current scan\n",
    "        if current_edge == 'C 1s':\n",
    "            start_en = 280\n",
    "            end_en = 295\n",
    "            start_q = 0.01\n",
    "            end_q = 0.1\n",
    "        else:\n",
    "            start_en = np.nanmin(integrated_data_pol0.energy)\n",
    "            end_en = np.nanmax(integrated_data_pol0.energy)\n",
    "            start_q = 0.01\n",
    "            end_q = np.nanmax(integrated_data_pol0.q)\n",
    "            \n",
    "        # calculate total intensities to normalize by\n",
    "        total_intensity_pol0 = integrated_data_pol0_unstacked.sel(q = slice(start_q, end_q,), polarization = 0).mean('q').mean('energy').mean('chi')\n",
    "        total_intensity_pol90 = integrated_data_pol90_unstacked.sel(q = slice(start_q, end_q,), polarization = 90).mean('q').mean('energy').mean('chi')\n",
    "        \n",
    "        total_intensity_pol0_corr = (total_intensity_pol0 + total_intensity_pol90)/2/total_intensity_pol0\n",
    "        total_intensity_pol90_corr = (total_intensity_pol0 + total_intensity_pol90)/2/total_intensity_pol90\n",
    "        \n",
    "        clear_output(wait = True)\n",
    "        \n",
    "\n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "           \n",
    "        ## Make plot of the AR\n",
    "        pol0_para = integrated_data_pol0_unstacked.rsoxs.slice_chi(90, chi_width = 22.5).sel(q = slice(start_q, end_q,), polarization = 0)*total_intensity_pol0_corr\n",
    "        pol0_perp = integrated_data_pol0_unstacked.rsoxs.slice_chi(180, chi_width = 22.5).sel(q = slice(start_q, end_q), polarization = 0)*total_intensity_pol0_corr\n",
    "\n",
    "        pol90_para = integrated_data_pol90_unstacked.rsoxs.slice_chi(180, chi_width = 22.5).sel(q = slice(start_q, end_q), polarization = 90)*total_intensity_pol90_corr\n",
    "        pol90_perp = integrated_data_pol90_unstacked.rsoxs.slice_chi(90, chi_width = 22.5).sel(q = slice(start_q, end_q), polarization = 90)*total_intensity_pol90_corr\n",
    "\n",
    "        para = pol0_para + pol90_para\n",
    "        perp = pol0_perp + pol90_perp\n",
    "\n",
    "        AR_pol0 = (pol0_para - pol0_perp)/(pol0_para + pol0_para)\n",
    "        AR_pol90 = (pol90_para - pol90_perp)/(pol90_para + pol90_para)\n",
    "        \n",
    "        AR = (AR_pol0 + AR_pol90)/2\n",
    "\n",
    "        scale = np.nanmax([abs(np.nanpercentile(AR.sel(energy = slice(280, 295)), 1)), abs(np.nanpercentile(AR.sel(energy = slice(start_en, end_en)), 99))])\n",
    "        # scale = np.nanmax([abs(np.nanmin(AR.sel(energy = slice(280, 295)))), abs(np.nanmax(AR.sel(energy = slice(start_en, end_en))))])\n",
    "        \n",
    "        AR_scales.append(scale)\n",
    "    \n",
    "        # log info for mosaic\n",
    "        ARs.append(AR)    \n",
    "\n",
    "        [fig, axes] = setup_plot()\n",
    "\n",
    "        im_cbar = AR.sel(energy = slice(start_en, end_en)).plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axes,\n",
    "            vmin = -scale,\n",
    "            vmax = scale,\n",
    "            cmap = 'RdBu',\n",
    "            add_colorbar = False)\n",
    "\n",
    "        cbar_ax = fig.add_axes([0.975, 0.145, 0.03, 0.715])\n",
    "        cbar = fig.colorbar(im_cbar, cax = cbar_ax)\n",
    "\n",
    "        # plt.setp(cbar.ax.get_yticklabels(), ha = \"right\")\n",
    "        # cbar.ax.tick_params(pad = 37.5)\n",
    "        cbar.ax.set_ylabel('Anisotropy')\n",
    "\n",
    "        fig.suptitle(f'{current_scan_pol0.sample_name}, {current_edge} edge', fontsize = 20, y = 1, x = 0.5)\n",
    "        axes.set_xlabel('$\\it{q}$ (Å$^{-1}$)')\n",
    "        axes.set_ylabel('Energy (eV)')\n",
    "        axes.set_xlim([start_q, end_q])\n",
    "        # axes.set_xticks([0.01, 0.03, 0.05, 0.07, 0.09], ['', '0.03', '0.05', '0.07', '0.09'])\n",
    "        axes.yaxis.set_major_formatter(ticker.FormatStrFormatter('%.2f'))\n",
    "        axes.set_ylim([start_en, end_en])\n",
    "\n",
    "        axes.xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', top = 'on')\n",
    "        axes.xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', top = 'on')\n",
    "        axes.yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "        axes.yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "        plt.show()\n",
    "        \n",
    "        save_fig(fig, current_scan_pol0, current_energy, reduced_save_dir, 'AR', 450)\n",
    "\n",
    "\n",
    "        [fig, axes] = setup_plot()\n",
    "\n",
    "        im_cbar = para.plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axes,\n",
    "            vmin = np.nanpercentile(para, 97.5),\n",
    "            vmax = scale,\n",
    "            cmap = cmap,\n",
    "            add_colorbar = False)\n",
    "        \n",
    "        paras.append(para)\n",
    "        para_scales.append(np.nanpercentile(para, 97.5))\n",
    "\n",
    "        cbar_ax = fig.add_axes([0.975, 0.145, 0.03, 0.715])\n",
    "        fig.colorbar(im_cbar, cax = cbar_ax)\n",
    "\n",
    "        fig.suptitle(f'{current_scan_pol0.sample_name}, {current_edge} edge', fontsize = 20, y = 1, x = 0.5)\n",
    "        axes.set_title('')\n",
    "        axes.text(0.975, 0.91, f'Parallel', color = 'white', fontsize = 15, transform = axes.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "        axes.set_xlabel('$\\it{q}$ (Å$^{-1}$)')\n",
    "        axes.set_ylabel('Energy (eV)')\n",
    "        axes.set_xlim([start_q, end_q])\n",
    "        # axes.set_xticks([0.01, 0.03, 0.05, 0.07, 0.09], ['', '0.03', '0.05', '0.07', '0.09'])\n",
    "        axes.yaxis.set_major_formatter(ticker.FormatStrFormatter('%.2f'))\n",
    "        axes.set_ylim([start_en, end_en])\n",
    "\n",
    "        axes.xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', top = 'on')\n",
    "        axes.xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', top = 'on')\n",
    "        axes.yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "        axes.yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "        plt.show()\n",
    "        \n",
    "        save_fig(fig, current_scan_pol0, current_energy, reduced_save_dir, 'para', 450)\n",
    "\n",
    "\n",
    "        [fig, axes] = setup_plot()\n",
    "\n",
    "        im_cbar = perp.plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axes,\n",
    "            vmin = np.nanpercentile(perp, 97.5),\n",
    "            vmax = scale,\n",
    "            cmap = cmap,\n",
    "            add_colorbar = False)\n",
    "        \n",
    "        perps.append(perp)\n",
    "        perp_scales.append(np.nanpercentile(perp, 97.5))\n",
    "\n",
    "        cbar_ax = fig.add_axes([0.975, 0.145, 0.03, 0.715])\n",
    "        fig.colorbar(im_cbar, cax = cbar_ax)\n",
    "\n",
    "        fig.suptitle(f'{current_scan_pol0.sample_name}, {current_edge} edge', fontsize = 20, y = 1, x = 0.5)\n",
    "        axes.set_title('')\n",
    "        axes.text(0.975, 0.91, f'Perpendicular', color = 'white', fontsize = 15, transform = axes.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "        axes.set_xlabel('$\\it{q}$ (Å$^{-1}$)')\n",
    "        axes.set_ylabel('Energy (eV)')\n",
    "        axes.set_xlim([start_q, end_q])\n",
    "        # axes.set_xticks([0.01, 0.03, 0.05, 0.07, 0.09], ['', '0.03', '0.05', '0.07', '0.09'])\n",
    "        axes.yaxis.set_major_formatter(ticker.FormatStrFormatter('%.2f'))\n",
    "        axes.set_ylim([start_en, end_en])\n",
    "\n",
    "        axes.xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', top = 'on')\n",
    "        axes.xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', top = 'on')\n",
    "        axes.yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "        axes.yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "        plt.show()\n",
    "        \n",
    "        save_fig(fig, current_scan_pol0, current_energy, reduced_save_dir, 'perp', 450)\n",
    "\n",
    "\n",
    "        para_perp = para + perp\n",
    "\n",
    "        [fig, axes] = setup_plot()\n",
    "\n",
    "        im_cbar = para_perp.plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axes,\n",
    "            vmin = np.nanpercentile(para_perp, 97.5),\n",
    "            vmax = scale,\n",
    "            cmap = cmap,\n",
    "            add_colorbar = False)\n",
    "        \n",
    "        circs.append(para_perp)\n",
    "        circ_scales.append(np.nanpercentile(para_perp, 97.5))\n",
    "        \n",
    "        cbar_ax = fig.add_axes([0.975, 0.145, 0.03, 0.715])\n",
    "        fig.colorbar(im_cbar, cax = cbar_ax)\n",
    "\n",
    "        fig.suptitle(f'{current_scan_pol0.sample_name}, {current_edge} edge', fontsize = 20, y = 1, x = 0.5)\n",
    "        axes.set_title('')\n",
    "        axes.text(0.975, 0.91, f'Circular', color = 'white', fontsize = 15, transform = axes.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "        axes.set_xlabel('$\\it{q}$ (Å$^{-1}$)')\n",
    "        axes.set_ylabel('Energy (eV)')\n",
    "        axes.set_xlim([start_q, end_q])\n",
    "        # axes.set_xticks([0.01, 0.03, 0.05, 0.07, 0.09], ['', '0.03', '0.05', '0.07', '0.09'])\n",
    "        axes.yaxis.set_major_formatter(ticker.FormatStrFormatter('%.2f'))\n",
    "        axes.set_ylim([start_en, end_en])\n",
    "\n",
    "        axes.xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', top = 'on')\n",
    "        axes.xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', top = 'on')\n",
    "        axes.yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "        axes.yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "        plt.show()\n",
    "        \n",
    "        save_fig(fig, current_scan_pol0, current_energy, reduced_save_dir, 'circ', 450)\n",
    "\n",
    "\n",
    "        integrated_data_unstacked = (integrated_data_pol0_unstacked.sel(polarization = 0)*total_intensity_pol0_corr) + (integrated_data_pol90_unstacked.sel(polarization = 90)*total_intensity_pol90_corr)\n",
    "\n",
    "        [fig, axes] = setup_plot()\n",
    "\n",
    "        FY_NEXAFS = integrated_data_unstacked.mean('chi').sel(q = slice(0.055, 0.075), energy = slice(start_en, end_en)).mean('q')\n",
    "        \n",
    "        FY_NEXAFS.plot(\n",
    "            ax = axes,\n",
    "            color = 'k')\n",
    "\n",
    "        fig.suptitle(f'{current_scan_pol0.sample_name}, {current_edge} edge', fontsize = 20, y = 1, x = 0.5)\n",
    "        axes.set_ylabel('Intensity (a.u.)')\n",
    "        axes.set_xlabel('Energy (eV)')\n",
    "        axes.set_xlim([start_en, end_en])\n",
    "        axes.set_ylim([0.9*np.nanmin(FY_NEXAFS), \n",
    "                       1.1*np.nanmax(FY_NEXAFS)])\n",
    "\n",
    "        plt.show()\n",
    "        \n",
    "        save_fig(fig, current_scan_pol0, current_energy, reduced_save_dir, 'FY_NEXAFS', 450)\n",
    "\n",
    "\n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        # variables to hist energy-dependent contrasts and plot limits\n",
    "        current_vmax_detImage = []\n",
    "        current_linecuts_ylim_min = []\n",
    "        current_linecuts_ylim_max = []\n",
    "        \n",
    "        ### Loop for iterating through energies\n",
    "        for current_energy in np.sort(integrated_data_pol0_unstacked.sel(energy = slice(start_en, end_en)).energy):\n",
    "\n",
    "            ## Make plots of detector image            \n",
    "            current_scan_pol0_unstacked = current_scan_pol0.unstack('system')\n",
    "            current_scan_pol90_unstacked = current_scan_pol90.unstack('system')\n",
    "\n",
    "            current_vmax_detImage.append(np.min([np.nanpercentile(current_scan_pol0_unstacked.sel(energy = current_energy)*total_intensity_pol0_corr, 97.5), \n",
    "                                        np.nanpercentile(current_scan_pol90_unstacked.sel(energy = current_energy)*total_intensity_pol90_corr, 97.5)]))\n",
    "\n",
    "            current_vmax_cake = current_vmax_detImage\n",
    "\n",
    "            fig, axes = plt.subplots(2, 2, figsize = (7, 7))\n",
    "\n",
    "            # im_pol0 = axes[0, 0].imshow(current_scan_pol0_unstacked.sel(energy = current_energy, polarization = 0).where(~np.array(integ.mask)),\n",
    "            im_pol0 = axes[0, 0].imshow(current_scan_pol0_unstacked.sel(energy = current_energy, polarization = 0)*total_intensity_pol0_corr,\n",
    "                            # norm = LogNorm(1, np.nanmax(current_vmax_detImage)),\n",
    "                            vmin = 0,\n",
    "                            vmax = current_vmax_detImage[-1],\n",
    "                            cmap = cmap,\n",
    "                            aspect = 1)\n",
    "\n",
    "            axes[0, 0].set_xticks([])\n",
    "            axes[0, 0].set_yticks([])\n",
    "\n",
    "            # im_pol90 = axes[0, 1].imshow(current_scan_pol90_unstacked.sel(energy = current_energy, polarization = 90).where(~np.array(integ.mask)),\n",
    "            im_pol90 = axes[0, 1].imshow(current_scan_pol90_unstacked.sel(energy = current_energy, polarization = 90)*total_intensity_pol90_corr,\n",
    "                            # norm = LogNorm(1, np.nanmax(current_vmax_detImage)),\n",
    "                            vmin = 0,\n",
    "                            vmax = current_vmax_detImage[-1],\n",
    "                            cmap = cmap,\n",
    "                            aspect = 1)\n",
    "\n",
    "            axes[0, 1].set_xticks([])\n",
    "            axes[0, 1].set_yticks([])\n",
    "\n",
    "            cake_pol0 = integrated_data_pol0_unstacked.sel(energy = current_energy, q = slice(start_q, end_q))*total_intensity_pol0_corr\n",
    "            \n",
    "            im = cake_pol0.plot(\n",
    "                x = 'chi', y = 'q',\n",
    "                ax = axes[1, 0],\n",
    "                # norm = LogNorm(1, current_vmax_detImage[-1]),\n",
    "                vmin = 0,\n",
    "                vmax = current_vmax_detImage[-1],\n",
    "                cmap = cmap,\n",
    "                add_colorbar = False)\n",
    "\n",
    "            axes[1, 0].set_ylabel('$\\it{q}$ (Å$^{-1}$)')\n",
    "            axes[1, 0].set_xlabel('Azimuth, $\\it{χ}$')\n",
    "            # axes[1, 0].set_yticks([0.01, 0.03, 0.05, 0.07, 0.09])\n",
    "            axes[1, 0].set_xticks([-180, -90, 0, 90, 180], ['', '-90°', '0°', '90°', ''])\n",
    "            # axes[1, 0].set_xticks([-180, -90, 0, 90, 180], ['-π', '', '0', '', 'π'])\n",
    "\n",
    "            axes[1, 0].xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', top = 'on')\n",
    "            axes[1, 0].xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', top = 'on')\n",
    "            axes[1, 0].yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'inout', right = 'on')\n",
    "            axes[1, 0].yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'inout', right = 'on')\n",
    "            \n",
    "            \n",
    "            cake_pol90 = integrated_data_pol90_unstacked.sel(energy = current_energy, q = slice(start_q, end_q))*total_intensity_pol90_corr\n",
    "\n",
    "            im = cake_pol90.plot(\n",
    "                x = 'chi', y = 'q',\n",
    "                ax = axes[1, 1],\n",
    "                # norm = LogNorm(1, current_vmax_detImage[-1]),\n",
    "                vmin = 0,\n",
    "                vmax = current_vmax_detImage[-1],\n",
    "                cmap = cmap,\n",
    "                add_colorbar = False)\n",
    "\n",
    "            axes[1, 1].set_xlabel('Azimuth, $\\it{χ}$')\n",
    "            axes[1, 1].set_ylabel('')\n",
    "            axes[1, 1].set_xticks([-180, -90, 0, 90, 180], ['', '-90°', '0°', '90°', '180°'])\n",
    "            # axes[1, 1].set_xticks([-180, -90, 0, 90, 180], ['-π', '', '0', '', 'π'])\n",
    "            # axes[1, 1].set_yticks([0.01, 0.03, 0.05, 0.07, 0.09])\n",
    "            axes[1, 1].set_yticklabels([])\n",
    "\n",
    "            axes[1, 1].xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', top = 'on')\n",
    "            axes[1, 1].xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', top = 'on')\n",
    "            axes[1, 1].yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'inout', right = 'on')\n",
    "            axes[1, 1].yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'inout', right = 'on')\n",
    "\n",
    "            cbar_ax = fig.add_axes([0.975, 0.145, 0.02, 0.725])\n",
    "            cbar = fig.colorbar(im_pol90, cax = cbar_ax)\n",
    "            # cbar.ax.set_ylabel('Intensity (a.u.)')\n",
    "\n",
    "            fig.suptitle(f'{current_scan_pol0.sample_name}, {current_edge} edge', fontsize = 17.5, y = 0.95, x = 0.55)\n",
    "            axes[0, 0].set_title('')\n",
    "            axes[0, 1].set_title('')\n",
    "            axes[1, 0].set_title('')\n",
    "            axes[1, 1].set_title('')\n",
    "\n",
    "            axes[0, 0].text(0.975, 0.91, f'{current_energy:.2f} eV, 0°', color = 'white', fontsize = 15, transform = axes[0, 0].transAxes, horizontalalignment = 'right')\n",
    "            axes[0, 1].text(0.975, 0.91, f'{current_energy:.2f} eV, 90°', color = 'white', fontsize = 15, transform = axes[0, 1].transAxes, horizontalalignment = 'right')\n",
    "            axes[1, 0].text(0.975, 0.91, f'{current_energy:.2f} eV, 0°', color = 'white', fontsize = 15, transform = axes[1, 0].transAxes, horizontalalignment = 'right')\n",
    "            axes[1, 1].text(0.975, 0.91, f'{current_energy:.2f} eV, 90°', color = 'white', fontsize = 15, transform = axes[1, 1].transAxes, horizontalalignment = 'right')\n",
    "\n",
    "            plt.tight_layout()\n",
    "\n",
    "            # plt.show()\n",
    "\n",
    "            save_fig(fig, current_scan_pol0, current_energy, rel_detImage_save_dir, 'detImage', 450)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "            [fig, axes] = plt.subplots(2, 2, figsize = (7, 7))\n",
    "            fig.suptitle(f'{current_scan_pol0.sample_name}, {current_edge} edge', fontsize = 20, y = 0.975, x = 0.5)\n",
    "            axes[1, 1].axis('off')\n",
    "\n",
    "            pol0_para = integrated_data_pol0_unstacked.rsoxs.slice_chi(90, chi_width = 22.5).sel(energy = current_energy, q = slice(start_q, end_q,), polarization = 0)*total_intensity_pol0_corr\n",
    "            pol0_perp = integrated_data_pol0_unstacked.rsoxs.slice_chi(180, chi_width = 22.5).sel(energy = current_energy, q = slice(start_q, end_q), polarization = 0)*total_intensity_pol0_corr\n",
    "\n",
    "            pol90_para = integrated_data_pol90_unstacked.rsoxs.slice_chi(180, chi_width = 22.5).sel(energy = current_energy, q = slice(start_q, end_q), polarization = 90)*total_intensity_pol90_corr\n",
    "            pol90_perp = integrated_data_pol90_unstacked.rsoxs.slice_chi(90, chi_width = 22.5).sel(energy = current_energy, q = slice(start_q, end_q), polarization = 90)*total_intensity_pol90_corr\n",
    "\n",
    "            para = pol0_para + pol90_para\n",
    "            perp = pol0_perp + pol90_perp\n",
    "\n",
    "            axes[0, 0].errorbar(np.asarray(para.q).flatten(), np.asarray(para).flatten(), np.sqrt(para.values).flatten(),\n",
    "                        linestyle ='none', color = (57/255, 115/255, 180/255), alpha=0.25)\n",
    "\n",
    "            axes[0, 0].plot(para.q, para, '-', color = (57/255, 115/255, 180/255), label = 'Parallel')\n",
    "\n",
    "            axes[0, 0].errorbar(np.asarray(perp.q).flatten(), np.asarray(perp).flatten(), np.sqrt(perp.values).flatten(),\n",
    "                        linestyle ='none', color = (170/255, 21/255, 40/255), alpha=0.25)\n",
    "\n",
    "            axes[0, 0].plot(perp.q, perp, '-', color = (170/255, 21/255, 40/255), label = 'Perpendicular')\n",
    "\n",
    "            axes[0, 0].legend(loc = 'upper right', frameon = False, fontsize = 15)\n",
    "\n",
    "            axes[0, 0].set_title('')\n",
    "            axes[0, 0].text(0.395, 0.6, f'{current_energy:.2f} eV', color = 'black', fontsize = 15, transform = axes[0, 0].transAxes, horizontalalignment = 'left')\n",
    "\n",
    "            axes[0, 0].set_ylabel('Intensity (a.u.)')\n",
    "            axes[0, 0].set_xlabel('')\n",
    "            axes[0, 0].set_xlim([start_q, end_q])\n",
    "            axes[0, 0].set_ylim([np.nanmin([np.nanpercentile(para, 1), np.nanpercentile(perp, 1)]),\n",
    "                     np.nanmax([np.nanpercentile(para, 99.9), np.nanpercentile(perp, 99.9)])])\n",
    "            axes[0, 0].set_yticklabels([])\n",
    "            axes[0, 0].set_xticklabels([])\n",
    "            # axes[0, 0].set_xticks([0.01, 0.03, 0.05, 0.07, 0.09])\n",
    "\n",
    "            axes[0, 0].yaxis.set_major_formatter(ticker.FuncFormatter(fmt))\n",
    "            axes[0, 0].set_yscale('log')\n",
    "\n",
    "            axes[0, 0].xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'inout', top = 'on')\n",
    "            axes[0, 0].xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'inout', top = 'on')\n",
    "            axes[0, 0].yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "            axes[0, 0].yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "            current_linecuts_ylim_min.append(np.nanmin(axes[0, 0].get_ylim()))\n",
    "            current_linecuts_ylim_max.append(np.nanmax(axes[0, 0].get_ylim()))\n",
    "\n",
    "\n",
    "            FY_NEXAFS.plot(\n",
    "                ax = axes[0, 1],\n",
    "                color = 'k',\n",
    "                alpha = 0.2)\n",
    "\n",
    "            integrated_data_unstacked.mean('chi').sel(q = slice(0.055, 0.075), energy = slice(start_en, current_energy)).mean('q').plot(\n",
    "                ax = axes[0, 1],\n",
    "                color = 'k',\n",
    "                alpha = 1)\n",
    "\n",
    "            axes[0, 1].plot(current_energy, integrated_data_unstacked.mean('chi').sel(q = slice(0.055, 0.075)).mean('q').sel(energy = current_energy), 'o', color = 'black')\n",
    "\n",
    "            axes[0, 1].set_title('Fluorescence NEXAFS', color = 'black', fontsize = 18)\n",
    "            # axes[0, 1].text(0.975, 0.05, f'FY NEXAFS', color = 'black', fontsize = 15, transform = axes[0, 1].transAxes, horizontalalignment = 'right')\n",
    "\n",
    "            # axes[0, 1].set_ylabel('Intensity (a.u.)')\n",
    "            axes[0, 1].set_xlabel('Energy (eV)')\n",
    "            axes[0, 1].set_xlim([start_en, end_en])\n",
    "            axes[0, 1].set_ylim([0.9*np.nanmin(FY_NEXAFS), \n",
    "                                 1.1*np.nanmax(FY_NEXAFS)])\n",
    "            axes[0, 1].set_yticklabels([])\n",
    "            # axes[0, 1].set_yscale('log')\n",
    "\n",
    "            axes[0, 1].xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', top = 'on')\n",
    "            axes[0, 1].xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', top = 'on')\n",
    "            axes[0, 1].yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "            axes[0, 1].yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "\n",
    "            pol0_para = integrated_data_pol0_unstacked.rsoxs.slice_chi(90, chi_width = 22.5).sel(energy = slice(start_en, current_energy), q = slice(start_q, end_q,), polarization = 0)*total_intensity_pol0_corr\n",
    "            pol0_perp = integrated_data_pol0_unstacked.rsoxs.slice_chi(180, chi_width = 22.5).sel(energy = slice(start_en, current_energy), q = slice(start_q, end_q), polarization = 0)*total_intensity_pol0_corr\n",
    "\n",
    "            pol90_para = integrated_data_pol90_unstacked.rsoxs.slice_chi(180, chi_width = 22.5).sel(energy = slice(start_en, current_energy), q = slice(start_q, end_q), polarization = 90)*total_intensity_pol90_corr\n",
    "            pol90_perp = integrated_data_pol90_unstacked.rsoxs.slice_chi(90, chi_width = 22.5).sel(energy = slice(start_en, current_energy), q = slice(start_q, end_q), polarization = 90)*total_intensity_pol90_corr\n",
    "\n",
    "            AR_pol0 = (pol0_para - pol0_perp)/(pol0_para + pol0_para)\n",
    "            AR_pol90 = (pol90_para - pol90_perp)/(pol90_para + pol90_para)\n",
    "\n",
    "            AR = (AR_pol0 + AR_pol90)/2\n",
    "            \n",
    "            if len(AR.energy) < 2:\n",
    "                pol0_para = integrated_data_pol0_unstacked.rsoxs.slice_chi(90, chi_width = 22.5).sel(energy = slice(start_en, end_en), q = slice(start_q, end_q), polarization = 0)*total_intensity_pol0_corr\n",
    "                pol0_perp = integrated_data_pol0_unstacked.rsoxs.slice_chi(180, chi_width = 22.5).sel(energy = slice(start_en, end_en), q = slice(start_q, end_q), polarization = 0)*total_intensity_pol0_corr\n",
    "\n",
    "                pol90_para = integrated_data_pol90_unstacked.rsoxs.slice_chi(180, chi_width = 22.5).sel(energy = slice(start_en, end_en), q = slice(start_q, end_q), polarization = 90)*total_intensity_pol90_corr\n",
    "                pol90_perp = integrated_data_pol90_unstacked.rsoxs.slice_chi(90, chi_width = 22.5).sel(energy = slice(start_en, end_en), q = slice(start_q, end_q), polarization = 90)*total_intensity_pol90_corr\n",
    "\n",
    "                AR_pol0 = (pol0_para - pol0_perp)/(pol0_para + pol0_para)\n",
    "                AR_pol90 = (pol90_para - pol90_perp)/(pol90_para + pol90_para)\n",
    "\n",
    "                AR1 = (AR_pol0 + AR_pol90)/2\n",
    "\n",
    "                scale = np.nanmax([abs(np.nanpercentile(AR1.sel(energy = slice(start_en, end_en)), 1)), abs(np.nanpercentile(AR1.sel(energy = slice(start_en, end_en)), 99))])\n",
    "                # scale = np.nanmax([abs(np.nanmin(AR.sel(energy = slice(280, 295)))), abs(np.nanmax(AR.sel(energy = slice(280, 295))))])\n",
    "\n",
    "                im_cbar = AR1.plot(\n",
    "                    x = 'q', y = 'energy',\n",
    "                    ax = axes[1, 0],\n",
    "                    vmin = -scale,\n",
    "                    vmax = scale,\n",
    "                    cmap = 'RdBu',\n",
    "                    add_colorbar = False)\n",
    "\n",
    "                cbar_ax = fig.add_axes([0.5, 0.14, 0.02, 0.31])\n",
    "                cbar = fig.colorbar(im_cbar, cax = cbar_ax)\n",
    "\n",
    "                # plt.setp(cbar.ax.get_yticklabels(), ha = \"right\")\n",
    "                # cbar.ax.tick_params(pad = 50)\n",
    "                # cbar.ax.set_ylabel('Anisotropy')\n",
    "\n",
    "                axes[1, 0].set_title('')\n",
    "                # axes[1, 0].text(0.975, 0.05, f'AR', color = 'black', fontsize = 15, transform = axes[1, 0].transAxes, horizontalalignment = 'right')\n",
    "\n",
    "                axes[1, 0].set_xlabel('$\\it{q}$ (Å$^{-1}$)')\n",
    "                axes[1, 0].set_ylabel('Energy (eV)')\n",
    "                axes[1, 0].set_xlim([start_q, end_q])\n",
    "                # axes[1, 0].set_xticks([0.01, 0.03, 0.05, 0.07, 0.09], ['', '0.03', '0.05', '0.07', '0.09'])\n",
    "                axes[1, 0].set_ylim([AR1.energy[0], (AR1.energy[1] - AR1.energy[0])/2 + AR1.energy[0]])\n",
    "                # axes[1, 0].set_ylim([np.nanmin(AR1.energy), np.nanmax(AR1.energy)])\n",
    "                axes[1, 0].yaxis.set_major_formatter(ticker.FormatStrFormatter('%.2f'))\n",
    "\n",
    "                # ticks = [tick for tick in axes[1, 0].get_yticks()]\n",
    "                # del ticks[-1]\n",
    "                # axes[1, 0].set_yticks(ticks)\n",
    "\n",
    "                axes[1, 0].xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'inout', top = 'on')\n",
    "                axes[1, 0].xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'inout', top = 'on')\n",
    "                axes[1, 0].yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "                axes[1, 0].yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "            else:\n",
    "\n",
    "                AR1.plot(\n",
    "                    x = 'q', y = 'energy',\n",
    "                    ax = axes[1, 0],\n",
    "                    vmin = -scale,\n",
    "                    vmax = scale,\n",
    "                    cmap = 'RdBu',\n",
    "                    add_colorbar = False,\n",
    "                    alpha = 0.2)\n",
    "\n",
    "                im_cbar = AR.plot(\n",
    "                    x = 'q', y = 'energy',\n",
    "                    ax = axes[1, 0],\n",
    "                    vmin = -scale,\n",
    "                    vmax = scale,\n",
    "                    cmap = 'RdBu',\n",
    "                    add_colorbar = False)\n",
    "\n",
    "                cbar_ax = fig.add_axes([0.5, 0.14, 0.02, 0.31])\n",
    "                cbar = fig.colorbar(im_cbar, cax = cbar_ax)\n",
    "\n",
    "                # plt.setp(cbar.ax.get_yticklabels(), ha = \"right\")\n",
    "                # cbar.ax.tick_params(pad = 50)\n",
    "                cbar.ax.set_ylabel('Anisotropy')\n",
    "\n",
    "                axes[1, 0].set_title('')\n",
    "                # axes[1, 0].text(0.975, 0.05, f'AR', color = 'black', fontsize = 15, transform = axes[1, 0].transAxes, horizontalalignment = 'right')\n",
    "\n",
    "                axes[1, 0].set_xlabel('$\\it{q}$ (Å$^{-1}$)')\n",
    "                axes[1, 0].set_ylabel('Energy (eV)')\n",
    "                axes[1, 0].set_xlim([start_q, end_q])\n",
    "                # axes[1, 0].set_xticks([0.01, 0.03, 0.05, 0.07, 0.09], ['', '0.03', '0.05', '0.07', '0.09'])\n",
    "\n",
    "                # ticks = [tick for tick in axes[1, 0].get_yticks()]\n",
    "                # del ticks[-1]\n",
    "                # axes[1, 0].set_yticks(ticks)\n",
    "                axes[1, 0].set_ylim([np.nanmin(AR.energy), np.nanmax(AR.energy)])\n",
    "                axes[1, 0].yaxis.set_major_formatter(ticker.FormatStrFormatter('%.2f'))\n",
    "\n",
    "                axes[1, 0].xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'inout', top = 'on')\n",
    "                axes[1, 0].xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'inout', top = 'on')\n",
    "                axes[1, 0].yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "                axes[1, 0].yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "            clear_output(wait = True)\n",
    "\n",
    "            # plt.tight_layout()\n",
    "\n",
    "            plt.show()\n",
    "            \n",
    "            save_fig(fig, current_scan_pol0, current_energy, rel_linecuts_save_dir, 'linecuts', 450)\n",
    "            \n",
    "            print(f'{len(scans)} scans remaining')\n",
    "            \n",
    "        clear_output(wait = True)    \n",
    "        \n",
    "        print(f'{len(scans)} scans remaining')\n",
    "        \n",
    "        print('Making gif of rel detector images')    \n",
    "        make_gif(current_scan_pol0, rel_detImage_save_dir, reduced_save_dir, 'rel_detImage')\n",
    "        \n",
    "        print('Making gif of rel linecuts images')    \n",
    "        make_gif(current_scan_pol0, rel_linecuts_save_dir, reduced_save_dir, 'rel_linecuts')\n",
    "        \n",
    "        print('rel gifs done, purging temporary images')  \n",
    "        shutil.rmtree(rel_detImage_save_dir)\n",
    "        shutil.rmtree(rel_linecuts_save_dir)\n",
    " \n",
    "        \n",
    "               \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "\n",
    "\n",
    "        ### Loop for iterating through energies\n",
    "        for current_energy in np.sort(integrated_data_pol0_unstacked.sel(energy = slice(start_en, end_en)).energy):\n",
    "\n",
    "            # ## Make plots of detector image            \n",
    "            # current_scan_pol0_unstacked = current_scan_pol0.unstack('system')\n",
    "            # current_scan_pol90_unstacked = current_scan_pol90.unstack('system')\n",
    "\n",
    "            # current_vmax_detImage.append(np.min([np.nanpercentile(current_scan_pol0_unstacked.sel(energy = current_energy)*total_intensity_pol0_corr, 97.5), \n",
    "            #                             np.nanpercentile(current_scan_pol90_unstacked.sel(energy = current_energy)*total_intensity_pol90_corr, 97.5)]))\n",
    "\n",
    "            current_vmax_cake = current_vmax_detImage\n",
    "\n",
    "            fig, axes = plt.subplots(2, 2, figsize = (7, 7))\n",
    "\n",
    "            # im_pol0 = axes[0, 0].imshow(current_scan_pol0_unstacked.sel(energy = current_energy, polarization = 0).where(~np.array(integ.mask)),\n",
    "            im_pol0 = axes[0, 0].imshow(current_scan_pol0_unstacked.sel(energy = current_energy, polarization = 0)*total_intensity_pol0_corr,\n",
    "                            # norm = LogNorm(1, np.nanmax(current_vmax_detImage)),\n",
    "                            vmin = 0,\n",
    "                            vmax = np.nanmax(current_vmax_detImage),\n",
    "                            cmap = cmap,\n",
    "                            aspect = 1)\n",
    "\n",
    "            axes[0, 0].set_xticks([])\n",
    "            axes[0, 0].set_yticks([])\n",
    "\n",
    "            # im_pol90 = axes[0, 1].imshow(current_scan_pol90_unstacked.sel(energy = current_energy, polarization = 90).where(~np.array(integ.mask)),\n",
    "            im_pol90 = axes[0, 1].imshow(current_scan_pol90_unstacked.sel(energy = current_energy, polarization = 90)*total_intensity_pol90_corr,\n",
    "                            # norm = LogNorm(1, np.nanmax(current_vmax_detImage)),\n",
    "                            vmin = 0,\n",
    "                            vmax = np.nanmax(current_vmax_detImage),\n",
    "                            cmap = cmap,\n",
    "                            aspect = 1)\n",
    "\n",
    "            axes[0, 1].set_xticks([])\n",
    "            axes[0, 1].set_yticks([])\n",
    "\n",
    "            cake_pol0 = integrated_data_pol0_unstacked.sel(energy = current_energy, q = slice(start_q, end_q))*total_intensity_pol0_corr\n",
    "            \n",
    "            im = cake_pol0.plot(\n",
    "                x = 'chi', y = 'q',\n",
    "                ax = axes[1, 0],\n",
    "                # norm = LogNorm(1, current_vmax_detImage[-1]),\n",
    "                vmin = 0,\n",
    "                vmax = np.nanmax(current_vmax_detImage),\n",
    "                cmap = cmap,\n",
    "                add_colorbar = False)\n",
    "\n",
    "            axes[1, 0].set_ylabel('$\\it{q}$ (Å$^{-1}$)')\n",
    "            axes[1, 0].set_xlabel('Azimuth, $\\it{χ}$')\n",
    "            # axes[1, 0].set_yticks([0.01, 0.03, 0.05, 0.07, 0.09])\n",
    "            axes[1, 0].set_xticks([-180, -90, 0, 90, 180], ['', '-90°', '0°', '90°', ''])\n",
    "            # axes[1, 0].set_xticks([-180, -90, 0, 90, 180], ['-π', '', '0', '', 'π'])\n",
    "\n",
    "            axes[1, 0].xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', top = 'on')\n",
    "            axes[1, 0].xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', top = 'on')\n",
    "            axes[1, 0].yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'inout', right = 'on')\n",
    "            axes[1, 0].yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'inout', right = 'on')\n",
    "            \n",
    "            \n",
    "            cake_pol90 = integrated_data_pol90_unstacked.sel(energy = current_energy, q = slice(start_q, end_q))*total_intensity_pol90_corr\n",
    "\n",
    "            im = cake_pol90.plot(\n",
    "                x = 'chi', y = 'q',\n",
    "                ax = axes[1, 1],\n",
    "                # norm = LogNorm(1, current_vmax_detImage[-1]),\n",
    "                vmin = 0,\n",
    "                vmax = np.nanmax(current_vmax_detImage),\n",
    "                cmap = cmap,\n",
    "                add_colorbar = False)\n",
    "\n",
    "            axes[1, 1].set_xlabel('Azimuth, $\\it{χ}$')\n",
    "            axes[1, 1].set_ylabel('')\n",
    "            axes[1, 1].set_xticks([-180, -90, 0, 90, 180], ['', '-90°', '0°', '90°', '180°'])\n",
    "            # axes[1, 1].set_xticks([-180, -90, 0, 90, 180], ['-π', '', '0', '', 'π'])\n",
    "            # axes[1, 1].set_yticks([0.01, 0.03, 0.05, 0.07, 0.09])\n",
    "            axes[1, 1].set_yticklabels([])\n",
    "\n",
    "            axes[1, 1].xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', top = 'on')\n",
    "            axes[1, 1].xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', top = 'on')\n",
    "            axes[1, 1].yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'inout', right = 'on')\n",
    "            axes[1, 1].yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'inout', right = 'on')\n",
    "\n",
    "            cbar_ax = fig.add_axes([0.975, 0.145, 0.02, 0.725])\n",
    "            cbar = fig.colorbar(im_pol90, cax = cbar_ax)\n",
    "            # cbar.ax.set_ylabel('Intensity (a.u.)')\n",
    "\n",
    "            fig.suptitle(f'{current_scan_pol0.sample_name}, {current_edge} edge', fontsize = 17.5, y = 0.95, x = 0.55)\n",
    "            axes[0, 0].set_title('')\n",
    "            axes[0, 1].set_title('')\n",
    "            axes[1, 0].set_title('')\n",
    "            axes[1, 1].set_title('')\n",
    "\n",
    "            axes[0, 0].text(0.975, 0.91, f'{current_energy:.2f} eV, 0°', color = 'white', fontsize = 15, transform = axes[0, 0].transAxes, horizontalalignment = 'right')\n",
    "            axes[0, 1].text(0.975, 0.91, f'{current_energy:.2f} eV, 90°', color = 'white', fontsize = 15, transform = axes[0, 1].transAxes, horizontalalignment = 'right')\n",
    "            axes[1, 0].text(0.975, 0.91, f'{current_energy:.2f} eV, 0°', color = 'white', fontsize = 15, transform = axes[1, 0].transAxes, horizontalalignment = 'right')\n",
    "            axes[1, 1].text(0.975, 0.91, f'{current_energy:.2f} eV, 90°', color = 'white', fontsize = 15, transform = axes[1, 1].transAxes, horizontalalignment = 'right')\n",
    "\n",
    "            plt.tight_layout()\n",
    "\n",
    "            # plt.show()\n",
    "\n",
    "            save_fig(fig, current_scan_pol0, current_energy, abs_detImage_save_dir, 'detImage', 450)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "            [fig, axes] = plt.subplots(2, 2, figsize = (7, 7))\n",
    "            fig.suptitle(f'{current_scan_pol0.sample_name}, {current_edge} edge', fontsize = 20, y = 0.975, x = 0.5)\n",
    "            axes[1, 1].axis('off')\n",
    "\n",
    "            pol0_para = integrated_data_pol0_unstacked.rsoxs.slice_chi(90, chi_width = 22.5).sel(energy = current_energy, q = slice(start_q, end_q,), polarization = 0)*total_intensity_pol0_corr\n",
    "            pol0_perp = integrated_data_pol0_unstacked.rsoxs.slice_chi(180, chi_width = 22.5).sel(energy = current_energy, q = slice(start_q, end_q), polarization = 0)*total_intensity_pol0_corr\n",
    "\n",
    "            pol90_para = integrated_data_pol90_unstacked.rsoxs.slice_chi(180, chi_width = 22.5).sel(energy = current_energy, q = slice(start_q, end_q), polarization = 90)*total_intensity_pol90_corr\n",
    "            pol90_perp = integrated_data_pol90_unstacked.rsoxs.slice_chi(90, chi_width = 22.5).sel(energy = current_energy, q = slice(start_q, end_q), polarization = 90)*total_intensity_pol90_corr\n",
    "\n",
    "            para = pol0_para + pol90_para\n",
    "            perp = pol0_perp + pol90_perp\n",
    "\n",
    "            axes[0, 0].errorbar(np.asarray(para.q).flatten(), np.asarray(para).flatten(), np.sqrt(para.values).flatten(),\n",
    "                        linestyle ='none', color = (57/255, 115/255, 180/255), alpha=0.25)\n",
    "\n",
    "            axes[0, 0].plot(para.q, para, '-', color = (57/255, 115/255, 180/255), label = 'Parallel')\n",
    "\n",
    "            axes[0, 0].errorbar(np.asarray(perp.q).flatten(), np.asarray(perp).flatten(), np.sqrt(perp.values).flatten(),\n",
    "                        linestyle ='none', color = (170/255, 21/255, 40/255), alpha=0.25)\n",
    "\n",
    "            axes[0, 0].plot(perp.q, perp, '-', color = (170/255, 21/255, 40/255), label = 'Perpendicular')\n",
    "\n",
    "            axes[0, 0].legend(loc = 'upper right', frameon = False, fontsize = 15)\n",
    "\n",
    "            axes[0, 0].set_title('')\n",
    "            axes[0, 0].text(0.395, 0.6, f'{current_energy:.2f} eV', color = 'black', fontsize = 15, transform = axes[0, 0].transAxes, horizontalalignment = 'left')\n",
    "\n",
    "            axes[0, 0].set_ylabel('Intensity (a.u.)')\n",
    "            axes[0, 0].set_xlabel('')\n",
    "            axes[0, 0].set_xlim([start_q, end_q])\n",
    "            axes[0, 0].set_ylim([np.nanmin(current_linecuts_ylim_min), \n",
    "                                 50*(np.nanmax(current_linecuts_ylim_max) - np.nanmin(current_linecuts_ylim_min)) + np.nanmin(current_linecuts_ylim_min)])\n",
    "            axes[0, 0].set_yticklabels([])\n",
    "            axes[0, 0].set_xticklabels([])\n",
    "            # axes[0, 0].set_xticks([0.01, 0.03, 0.05, 0.07, 0.09])\n",
    "\n",
    "            axes[0, 0].yaxis.set_major_formatter(ticker.FuncFormatter(fmt))\n",
    "            axes[0, 0].set_yscale('log')\n",
    "\n",
    "            axes[0, 0].xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'inout', top = 'on')\n",
    "            axes[0, 0].xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'inout', top = 'on')\n",
    "            axes[0, 0].yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "            axes[0, 0].yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "            # current_linecuts_ylim_min.append(np.nanmin(axes[0, 0].get_ylim()))\n",
    "            # current_linecuts_ylim_max.append(np.nanmax(axes[0, 0].get_ylim()))\n",
    "\n",
    "\n",
    "            FY_NEXAFS.plot(\n",
    "                ax = axes[0, 1],\n",
    "                color = 'k',\n",
    "                alpha = 0.2)\n",
    "\n",
    "            integrated_data_unstacked.mean('chi').sel(q = slice(0.055, 0.075), energy = slice(start_en, current_energy)).mean('q').plot(\n",
    "                ax = axes[0, 1],\n",
    "                color = 'k',\n",
    "                alpha = 1)\n",
    "\n",
    "            axes[0, 1].plot(current_energy, integrated_data_unstacked.mean('chi').sel(q = slice(0.055, 0.075)).mean('q').sel(energy = current_energy), 'o', color = 'black')\n",
    "\n",
    "            axes[0, 1].set_title('Fluorescence NEXAFS', color = 'black', fontsize = 18)\n",
    "            # axes[0, 1].text(0.975, 0.05, f'FY NEXAFS', color = 'black', fontsize = 15, transform = axes[0, 1].transAxes, horizontalalignment = 'right')\n",
    "\n",
    "            # axes[0, 1].set_ylabel('Intensity (a.u.)')\n",
    "            axes[0, 1].set_xlabel('Energy (eV)')\n",
    "            axes[0, 1].set_xlim([start_en, end_en])\n",
    "            axes[0, 1].set_ylim([0.9*np.nanmin(FY_NEXAFS), \n",
    "                                 1.1*np.nanmax(FY_NEXAFS)])\n",
    "            axes[0, 1].set_yticklabels([])\n",
    "            # axes[0, 1].set_yscale('log')\n",
    "\n",
    "            axes[0, 1].xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', top = 'on')\n",
    "            axes[0, 1].xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', top = 'on')\n",
    "            axes[0, 1].yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "            axes[0, 1].yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "\n",
    "            pol0_para = integrated_data_pol0_unstacked.rsoxs.slice_chi(90, chi_width = 22.5).sel(energy = slice(start_en, current_energy), q = slice(start_q, end_q,), polarization = 0)*total_intensity_pol0_corr\n",
    "            pol0_perp = integrated_data_pol0_unstacked.rsoxs.slice_chi(180, chi_width = 22.5).sel(energy = slice(start_en, current_energy), q = slice(start_q, end_q), polarization = 0)*total_intensity_pol0_corr\n",
    "\n",
    "            pol90_para = integrated_data_pol90_unstacked.rsoxs.slice_chi(180, chi_width = 22.5).sel(energy = slice(start_en, current_energy), q = slice(start_q, end_q), polarization = 90)*total_intensity_pol90_corr\n",
    "            pol90_perp = integrated_data_pol90_unstacked.rsoxs.slice_chi(90, chi_width = 22.5).sel(energy = slice(start_en, current_energy), q = slice(start_q, end_q), polarization = 90)*total_intensity_pol90_corr\n",
    "\n",
    "            AR_pol0 = (pol0_para - pol0_perp)/(pol0_para + pol0_para)\n",
    "            AR_pol90 = (pol90_para - pol90_perp)/(pol90_para + pol90_para)\n",
    "\n",
    "            AR = (AR_pol0 + AR_pol90)/2\n",
    "            \n",
    "            if len(AR.energy) < 2:\n",
    "                pol0_para = integrated_data_pol0_unstacked.rsoxs.slice_chi(90, chi_width = 22.5).sel(energy = slice(start_en, end_en), q = slice(start_q, end_q), polarization = 0)*total_intensity_pol0_corr\n",
    "                pol0_perp = integrated_data_pol0_unstacked.rsoxs.slice_chi(180, chi_width = 22.5).sel(energy = slice(start_en, end_en), q = slice(start_q, end_q), polarization = 0)*total_intensity_pol0_corr\n",
    "\n",
    "                pol90_para = integrated_data_pol90_unstacked.rsoxs.slice_chi(180, chi_width = 22.5).sel(energy = slice(start_en, end_en), q = slice(start_q, end_q), polarization = 90)*total_intensity_pol90_corr\n",
    "                pol90_perp = integrated_data_pol90_unstacked.rsoxs.slice_chi(90, chi_width = 22.5).sel(energy = slice(start_en, end_en), q = slice(start_q, end_q), polarization = 90)*total_intensity_pol90_corr\n",
    "\n",
    "                AR_pol0 = (pol0_para - pol0_perp)/(pol0_para + pol0_para)\n",
    "                AR_pol90 = (pol90_para - pol90_perp)/(pol90_para + pol90_para)\n",
    "\n",
    "                AR1 = (AR_pol0 + AR_pol90)/2\n",
    "\n",
    "                scale = np.nanmax([abs(np.nanpercentile(AR1.sel(energy = slice(start_en, end_en)), 1)), abs(np.nanpercentile(AR1.sel(energy = slice(start_en, end_en)), 99))])\n",
    "                # scale = np.nanmax([abs(np.nanmin(AR.sel(energy = slice(280, 295)))), abs(np.nanmax(AR.sel(energy = slice(280, 295))))])\n",
    "\n",
    "                im_cbar = AR1.plot(\n",
    "                    x = 'q', y = 'energy',\n",
    "                    ax = axes[1, 0],\n",
    "                    vmin = -scale,\n",
    "                    vmax = scale,\n",
    "                    cmap = 'RdBu',\n",
    "                    add_colorbar = False)\n",
    "\n",
    "                cbar_ax = fig.add_axes([0.5, 0.14, 0.02, 0.31])\n",
    "                cbar = fig.colorbar(im_cbar, cax = cbar_ax)\n",
    "\n",
    "                # plt.setp(cbar.ax.get_yticklabels(), ha = \"right\")\n",
    "                # cbar.ax.tick_params(pad = 50)\n",
    "                # cbar.ax.set_ylabel('Anisotropy')\n",
    "\n",
    "                axes[1, 0].set_title('')\n",
    "                # axes[1, 0].text(0.975, 0.05, f'AR', color = 'black', fontsize = 15, transform = axes[1, 0].transAxes, horizontalalignment = 'right')\n",
    "\n",
    "                axes[1, 0].set_xlabel('$\\it{q}$ (Å$^{-1}$)')\n",
    "                axes[1, 0].set_ylabel('Energy (eV)')\n",
    "                axes[1, 0].set_xlim([start_q, end_q])\n",
    "                # axes[1, 0].set_xticks([0.01, 0.03, 0.05, 0.07, 0.09], ['', '0.03', '0.05', '0.07', '0.09'])\n",
    "                axes[1, 0].set_ylim([AR1.energy[0], (AR1.energy[1] - AR1.energy[0])/2 + AR1.energy[0]])\n",
    "                # axes[1, 0].set_ylim([np.nanmin(AR1.energy), np.nanmax(AR1.energy)])\n",
    "                axes[1, 0].yaxis.set_major_formatter(ticker.FormatStrFormatter('%.2f'))\n",
    "\n",
    "                # ticks = [tick for tick in axes[1, 0].get_yticks()]\n",
    "                # del ticks[-1]\n",
    "                # axes[1, 0].set_yticks(ticks)\n",
    "\n",
    "                axes[1, 0].xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'inout', top = 'on')\n",
    "                axes[1, 0].xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'inout', top = 'on')\n",
    "                axes[1, 0].yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "                axes[1, 0].yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "            else:\n",
    "\n",
    "                AR1.plot(\n",
    "                    x = 'q', y = 'energy',\n",
    "                    ax = axes[1, 0],\n",
    "                    vmin = -scale,\n",
    "                    vmax = scale,\n",
    "                    cmap = 'RdBu',\n",
    "                    add_colorbar = False,\n",
    "                    alpha = 0.2)\n",
    "\n",
    "                im_cbar = AR.plot(\n",
    "                    x = 'q', y = 'energy',\n",
    "                    ax = axes[1, 0],\n",
    "                    vmin = -scale,\n",
    "                    vmax = scale,\n",
    "                    cmap = 'RdBu',\n",
    "                    add_colorbar = False)\n",
    "\n",
    "                cbar_ax = fig.add_axes([0.5, 0.14, 0.02, 0.31])\n",
    "                cbar = fig.colorbar(im_cbar, cax = cbar_ax)\n",
    "\n",
    "                # plt.setp(cbar.ax.get_yticklabels(), ha = \"right\")\n",
    "                # cbar.ax.tick_params(pad = 50)\n",
    "                cbar.ax.set_ylabel('Anisotropy')\n",
    "\n",
    "                axes[1, 0].set_title('')\n",
    "                # axes[1, 0].text(0.975, 0.05, f'AR', color = 'black', fontsize = 15, transform = axes[1, 0].transAxes, horizontalalignment = 'right')\n",
    "\n",
    "                axes[1, 0].set_xlabel('$\\it{q}$ (Å$^{-1}$)')\n",
    "                axes[1, 0].set_ylabel('Energy (eV)')\n",
    "                axes[1, 0].set_xlim([start_q, end_q])\n",
    "                # axes[1, 0].set_xticks([0.01, 0.03, 0.05, 0.07, 0.09], ['', '0.03', '0.05', '0.07', '0.09'])\n",
    "\n",
    "                # ticks = [tick for tick in axes[1, 0].get_yticks()]\n",
    "                # del ticks[-1]\n",
    "                # axes[1, 0].set_yticks(ticks)\n",
    "                axes[1, 0].set_ylim([np.nanmin(AR.energy), np.nanmax(AR.energy)])\n",
    "                axes[1, 0].yaxis.set_major_formatter(ticker.FormatStrFormatter('%.2f'))\n",
    "\n",
    "                axes[1, 0].xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'inout', top = 'on')\n",
    "                axes[1, 0].xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'inout', top = 'on')\n",
    "                axes[1, 0].yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "                axes[1, 0].yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "            clear_output(wait = True)\n",
    "\n",
    "            # plt.tight_layout()\n",
    "\n",
    "            plt.show()\n",
    "            \n",
    "            save_fig(fig, current_scan_pol0, current_energy, abs_linecuts_save_dir, 'linecuts', 450)\n",
    "            \n",
    "            print(f'{len(scans)} scans remaining')\n",
    "            \n",
    "        clear_output(wait = True)    \n",
    "        \n",
    "        print(f'{len(scans)} scans remaining')\n",
    "        \n",
    "        print('Making gif of abs det images')    \n",
    "        make_gif(current_scan_pol0, abs_detImage_save_dir, reduced_save_dir, 'abs_detImage')\n",
    "        \n",
    "        print('Making gif of abs linecuts images') \n",
    "        make_gif(current_scan_pol0, abs_linecuts_save_dir, reduced_save_dir, 'abs_linecuts')\n",
    "        \n",
    "        print('ans gifs done, purging temporary images')  \n",
    "        shutil.rmtree(abs_detImage_save_dir)\n",
    "        shutil.rmtree(abs_linecuts_save_dir)\n",
    "        \n",
    "        # save variables to pickle\n",
    "        with open('PyHyperScattering_Batch_SST1_JupyterHub.pkl', 'wb') as file:\n",
    "            pickle.dump([sample_names, edges, ARs, AR_scales, paras, para_scales, perps, perp_scales, circs, circ_scales], file, protocol = -1)\n",
    "\n",
    "    else:\n",
    "        print(f'{current_scan_pol0.sample_name} already processed, moving on')\n",
    "\n",
    "\n",
    "\n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "\n",
    "print(f'{len(scans)} scans leftover')"
   ]
  },
  {
   "cell_type": "raw",
   "id": "1d40771f-6f44-4d9d-97e3-cd3a16e7ff2a",
   "metadata": {},
   "source": [
    "with open('PyHyperScattering_Batch_SST1_JupyterHub_2022_10_11.pkl', 'rb') as file:\n",
    "    [sample_names, edges, ARs, AR_scales, paras, para_scales, perps, perp_scales, circs, circ_scales] = pickle.load(file)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "1e8e9d27-1e39-40c0-ae2e-be5fb633d3f9",
   "metadata": {},
   "source": [
    "with open('PyHyperScattering_Batch_SST1_JupyterHub_2022_10_11.pkl', 'wb') as file:\n",
    "    pickle.dump([sample_names, edges, ARs, AR_scales, paras, para_scales, perps, perp_scales, circs, circ_scales], file, protocol = -1))"
   ]
  },
  {
   "cell_type": "raw",
   "id": "72855202-5c19-4fa1-a9e4-c2ebad110e3e",
   "metadata": {},
   "source": [
    "dicts = dict(zip(sample_names, ARs))"
   ]
  },
  {
   "cell_type": "raw",
   "id": "e4dc0be0-93fd-46db-8cf2-157f9621749c",
   "metadata": {},
   "source": [
    "# create directories for output files\n",
    "mosaics_save_dir = os.path.join(base_save_dir, 'mosaics')\n",
    "os.makedirs(mosaics_save_dir, exist_ok = True)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "c4a11038-cf54-4f77-aa2b-8020e626f26e",
   "metadata": {},
   "source": [
    "[fig, axes] = plt.subplots(3, 6, figsize = (3.5*6, 3.5*3))\n",
    "\n",
    "y_ticks = axis.get_yticks()\n",
    "\n",
    "# format all axes in subplot\n",
    "for axis in axes.reshape(-1):\n",
    "    axis.set_xticklabels('')\n",
    "    axis.set_yticklabels('')\n",
    "    axis.set_xticks([0.01, 0.03, 0.05, 0.07, 0.09])\n",
    "    y_ticks = list(range(start_en, end_en + 1, 5))\n",
    "    axis.set_yticks(y_ticks)\n",
    "\n",
    "    axis.xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'inout', top = 'on')\n",
    "    axis.xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'inout', top = 'on')\n",
    "    axis.yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'inout', right = 'on')\n",
    "    axis.yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'inout', right = 'on')\n",
    "\n",
    "key = 'EPRA2'; axis = axes[2, 0]\n",
    "\n",
    "im_cbar = dicts[key].sel(energy = slice(start_en, end_en)).plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axis,\n",
    "            vmin = -np.nanmax(current_AR_scale),\n",
    "            vmax = np.nanmax(current_AR_scale),\n",
    "            cmap = 'RdBu',\n",
    "            add_colorbar = False)\n",
    "\n",
    "axis.text(0.965, 0.9, f'{key}', color = 'black', fontsize = 18, transform = axis.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "key = 'EPAN2'; axis = axes[0, 1]\n",
    "\n",
    "im_cbar = dicts[key].sel(energy = slice(start_en, end_en)).plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axis,\n",
    "            vmin = -np.nanmax(current_AR_scale),\n",
    "            vmax = np.nanmax(current_AR_scale),\n",
    "            cmap = 'RdBu',\n",
    "            add_colorbar = False)\n",
    "\n",
    "axis.text(0.965, 0.9, f'{key}', color = 'black', fontsize = 18, transform = axis.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "key = 'EPAV2'; axis = axes[1, 1]\n",
    "\n",
    "im_cbar = dicts[key].sel(energy = slice(start_en, end_en)).plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axis,\n",
    "            vmin = -np.nanmax(current_AR_scale),\n",
    "            vmax = np.nanmax(current_AR_scale),\n",
    "            cmap = 'RdBu',\n",
    "            add_colorbar = False)\n",
    "\n",
    "axis.text(0.965, 0.9, f'{key}', color = 'black', fontsize = 18, transform = axis.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "key = 'EPAA2'; axis = axes[2, 1]\n",
    "\n",
    "im_cbar = dicts[key].sel(energy = slice(start_en, end_en)).plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axis,\n",
    "            vmin = -np.nanmax(current_AR_scale),\n",
    "            vmax = np.nanmax(current_AR_scale),\n",
    "            cmap = 'RdBu',\n",
    "            add_colorbar = False)\n",
    "\n",
    "axis.text(0.965, 0.9, f'{key}', color = 'black', fontsize = 18, transform = axis.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "key = 'EPBN2'; axis = axes[0, 2]\n",
    "\n",
    "im_cbar = dicts[key].sel(energy = slice(start_en, end_en)).plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axis,\n",
    "            vmin = -np.nanmax(current_AR_scale),\n",
    "            vmax = np.nanmax(current_AR_scale),\n",
    "            cmap = 'RdBu',\n",
    "            add_colorbar = False)\n",
    "\n",
    "axis.text(0.965, 0.9, f'{key}', color = 'black', fontsize = 18, transform = axis.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "key = 'EPBV1'; axis = axes[1, 2]\n",
    "\n",
    "im_cbar = dicts[key].sel(energy = slice(start_en, end_en)).plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axis,\n",
    "            vmin = -np.nanmax(current_AR_scale),\n",
    "            vmax = np.nanmax(current_AR_scale),\n",
    "            cmap = 'RdBu',\n",
    "            add_colorbar = False)\n",
    "\n",
    "axis.text(0.965, 0.9, f'{key}', color = 'black', fontsize = 18, transform = axis.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "key = 'EPCN2'; axis = axes[0, 3]\n",
    "\n",
    "im_cbar = dicts[key].sel(energy = slice(start_en, end_en)).plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axis,\n",
    "            vmin = -np.nanmax(current_AR_scale),\n",
    "            vmax = np.nanmax(current_AR_scale),\n",
    "            cmap = 'RdBu',\n",
    "            add_colorbar = False)\n",
    "\n",
    "axis.text(0.965, 0.9, f'{key}', color = 'black', fontsize = 18, transform = axis.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "key = 'EPCV2'; axis = axes[1, 3]\n",
    "\n",
    "im_cbar = dicts[key].sel(energy = slice(start_en, end_en)).plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axis,\n",
    "            vmin = -np.nanmax(current_AR_scale),\n",
    "            vmax = np.nanmax(current_AR_scale),\n",
    "            cmap = 'RdBu',\n",
    "            add_colorbar = False)\n",
    "\n",
    "axis.text(0.965, 0.9, f'{key}', color = 'black', fontsize = 18, transform = axis.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "key = 'EPCA2'; axis = axes[2, 3]\n",
    "\n",
    "im_cbar = dicts[key].sel(energy = slice(start_en, end_en)).plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axis,\n",
    "            vmin = -np.nanmax(current_AR_scale),\n",
    "            vmax = np.nanmax(current_AR_scale),\n",
    "            cmap = 'RdBu',\n",
    "            add_colorbar = False)\n",
    "\n",
    "axis.text(0.965, 0.9, f'{key}', color = 'black', fontsize = 18, transform = axis.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "key = 'EPDN2'; axis = axes[0, 4]\n",
    "\n",
    "im_cbar = dicts[key].sel(energy = slice(start_en, end_en)).plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axis,\n",
    "            vmin = -np.nanmax(current_AR_scale),\n",
    "            vmax = np.nanmax(current_AR_scale),\n",
    "            cmap = 'RdBu',\n",
    "            add_colorbar = False)\n",
    "\n",
    "axis.text(0.965, 0.9, f'{key}', color = 'black', fontsize = 18, transform = axis.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "key = 'EPDV1'; axis = axes[1, 4]\n",
    "\n",
    "im_cbar = dicts[key].sel(energy = slice(start_en, end_en)).plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axis,\n",
    "            vmin = -np.nanmax(current_AR_scale),\n",
    "            vmax = np.nanmax(current_AR_scale),\n",
    "            cmap = 'RdBu',\n",
    "            add_colorbar = False)\n",
    "\n",
    "axis.text(0.965, 0.9, f'{key}', color = 'black', fontsize = 18, transform = axis.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "key = 'EPDA2'; axis = axes[2, 4]\n",
    "\n",
    "im_cbar = dicts[key].sel(energy = slice(start_en, end_en)).plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axis,\n",
    "            vmin = -np.nanmax(current_AR_scale),\n",
    "            vmax = np.nanmax(current_AR_scale),\n",
    "            cmap = 'RdBu',\n",
    "            add_colorbar = False)\n",
    "\n",
    "axis.text(0.965, 0.9, f'{key}', color = 'black', fontsize = 18, transform = axis.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "key = 'EPEN2'; axis = axes[0, 5]\n",
    "\n",
    "im_cbar = dicts[key].sel(energy = slice(start_en, end_en)).plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axis,\n",
    "            vmin = -np.nanmax(current_AR_scale),\n",
    "            vmax = np.nanmax(current_AR_scale),\n",
    "            cmap = 'RdBu',\n",
    "            add_colorbar = False)\n",
    "\n",
    "axis.text(0.965, 0.9, f'{key}', color = 'black', fontsize = 18, transform = axis.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "key = 'EPEV2'; axis = axes[1, 5]\n",
    "\n",
    "im_cbar = dicts[key].sel(energy = slice(start_en, end_en)).plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axis,\n",
    "            vmin = -np.nanmax(current_AR_scale),\n",
    "            vmax = np.nanmax(current_AR_scale),\n",
    "            cmap = 'RdBu',\n",
    "            add_colorbar = False)\n",
    "\n",
    "axis.text(0.965, 0.9, f'{key}', color = 'black', fontsize = 18, transform = axis.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "key = 'EPEA1'; axis = axes[2, 5]\n",
    "\n",
    "im_cbar = dicts[key].sel(energy = slice(start_en, end_en)).plot(\n",
    "            # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "            x = 'q', y = 'energy',\n",
    "            ax = axis,\n",
    "            vmin = -np.nanmax(current_AR_scale),\n",
    "            vmax = np.nanmax(current_AR_scale),\n",
    "            cmap = 'RdBu',\n",
    "            add_colorbar = False)\n",
    "\n",
    "axis.text(0.965, 0.9, f'{key}', color = 'black', fontsize = 18, transform = axis.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "# format empty axis\n",
    "\n",
    "key = 'EPRN'; axis = axes[0, 0]\n",
    "\n",
    "axis.text(0.965, 0.9, f'{key}', color = 'black', fontsize = 18, transform = axis.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "key = 'EPRV'; axis = axes[1, 0]\n",
    "\n",
    "axis.text(0.965, 0.9, f'{key}', color = 'black', fontsize = 18, transform = axis.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "key = 'EPBA'; axis = axes[2, 2]\n",
    "\n",
    "axis.text(0.965, 0.9, f'{key}', color = 'black', fontsize = 18, transform = axis.transAxes, horizontalalignment = 'right')\n",
    "\n",
    "# remove axes labels from xarray's plot function, reset limits\n",
    "for axis in axes.reshape(-1):\n",
    "    axis.set_ylabel('')\n",
    "    axis.set_xlabel('')\n",
    "    axis.set_ylim([start_en, end_en])\n",
    "    axis.set_xlim([start_q, end_q])\n",
    "\n",
    "# format bottom axis along bottom row\n",
    "\n",
    "for i in range(6):\n",
    "    axis = axes[2, i]\n",
    "    \n",
    "    axis.set_xlabel('$\\it{q}$ (Å$^{-1}$)')\n",
    "    axis.set_xticks([0.01, 0.03, 0.05, 0.07, 0.09], ['', '0.03', '0.05', '0.07', '0.09'])\n",
    "    \n",
    "# format left axis along left column\n",
    "\n",
    "for i in range(3):\n",
    "    axis = axes[i, 0]\n",
    "    \n",
    "    axis.set_ylabel('Energy (eV)')\n",
    "    \n",
    "    # y_ticks[0] = ''\n",
    "    \n",
    "    axis.set_yticklabels(y_ticks)\n",
    "\n",
    "fig.suptitle('C 1s edge', fontsize = 25, y = 1, x = 0.5)\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "cbar_ax = fig.add_axes([1, 0.11, 0.01, 0.79])\n",
    "cbar = fig.colorbar(im_cbar, cax = cbar_ax)\n",
    "cbar.ax.set_ylabel('Anisotropy')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "save_path = os.path.join(mosaics_save_dir, f'mosaic_{current_edge_folder}.png')\n",
    "fig.savefig(save_path, bbox_inches = 'tight', pad_inches = 0.2, dpi = 300, transparent = False)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "7be36bdc",
   "metadata": {
    "tags": []
   },
   "source": [
    "base_path_tracker = 0;\n",
    "for base_path in base_paths:\n",
    "    base_path_tracker += 1\n",
    "    scans = pd.Series(get_immediate_subdirectories(base_path))\n",
    "    for current_scan in np.sort(scans):\n",
    "        if ((current_scan >= '44065') and (current_scan != '44068') and (current_scan != '44109')):\n",
    "            clear_output(wait=True)\n",
    "            print(f'Currently pulling data from \\'{base_path}\\' ({base_path_tracker}/{len(base_paths)})')\n",
    "            print(f'Currently processing scan {current_scan}/{np.nanmax(scans)}')\n",
    "            scan = file_loader.loadFileSeries(base_path + current_scan,['energy','polarization'])\n",
    "            if scan.rsoxs_config == 'waxs':\n",
    "                maskmethod = 'nika'\n",
    "                mask = 'SST1_Sept2022_WAXS_mask.hdf'\n",
    "                # import mask\n",
    "                integ.loadNikaMask(filetoload=mask,rotate_image=False)\n",
    "\n",
    "                # rotate mask to match experimental orientation\n",
    "                integ.mask = np.flipud(np.rot90(integ.mask))\n",
    "\n",
    "                # compare experimental and mask orientation\n",
    "                # f, axarr = plt.subplots(1,2)\n",
    "                # axarr[0].imshow(scan_unstacked.sel(energy=np.max(scan.energy)),\n",
    "                #                 # norm = LogNorm(1,np.nanpercentile(scan_unstacked.sel(energy=current_energy),99)\n",
    "                #                 vmin = 0,\n",
    "                #                 vmax = np.nanpercentile(scan_unstacked.sel(energy=np.max(scan.energy)), 90),\n",
    "                #                 cmap = 'viridis',\n",
    "                #                 aspect='1')\n",
    "                # axarr[1].imshow(integ.mask,# norm = LogNorm(1,np.nanpercentile(scan_unstacked.sel(energy=np.max(scan.energy)),90)))\n",
    "\n",
    "            elif scan.rsoxs_config == 'saxs':\n",
    "                maskmethod = 'none'\n",
    "            else:\n",
    "                maskmethod = 'none'\n",
    "\n",
    "            scan_unstacked = scan.unstack('system')\n",
    "\n",
    "            if np.min(scan.energy) == 270:\n",
    "                current_edge = 'C 1s'\n",
    "                current_edge_folder = 'C1s'\n",
    "            elif np.min(scan.energy) == 385:\n",
    "                current_edge = 'N 1s'\n",
    "                current_edge_folder = 'N1s'\n",
    "            elif np.min(scan.energy) == 670:\n",
    "                current_edge = 'F 1s'\n",
    "                current_edge_folder = 'F1s'\n",
    "            else:\n",
    "                print('error: unknown edge')\n",
    "\n",
    "            # Create sample directory for processed data\n",
    "            if scan_unstacked.polarization == 0:\n",
    "                pol = 0\n",
    "            elif scan_unstacked.polarization == 90:\n",
    "                pol = 90\n",
    "            else:\n",
    "                print('Error in finding polarization')\n",
    "\n",
    "            sample_scan = f'{scan_unstacked.sample_name}_{current_edge_folder}_pol{pol}'\n",
    "\n",
    "            # create directories for output files\n",
    "            save_dir = os.path.join(base_save_dir, sample_scan)\n",
    "            os.makedirs(save_dir, exist_ok=True)\n",
    "\n",
    "            sample_scan_reduced = f'reduced_pol{pol}_{current_edge_folder}_{scan_unstacked.sample_name}'\n",
    "\n",
    "            # create directories for reduced output files\n",
    "            save_dir_reduced = os.path.join(base_save_dir, sample_scan_reduced)\n",
    "            os.makedirs(save_dir_reduced, exist_ok=True)\n",
    "\n",
    "            save_path = os.path.join(save_dir_reduced, sample_scan + '_AR.png')\n",
    "            tail = os.path.split(save_path)[1]\n",
    "\n",
    "            if tail not in os.listdir(save_dir_reduced):\n",
    "            # if True:\n",
    "                # Remove existing files in sample directory\n",
    "                for filename in os.listdir(save_dir):\n",
    "                    file_path = os.path.join(save_dir, filename)\n",
    "                    try:\n",
    "                        if os.path.isfile(file_path) or os.path.islink(file_path):\n",
    "                            os.unlink(file_path)\n",
    "                        elif os.path.isdir(file_path):\n",
    "                            shutil.rmtree(file_path)\n",
    "                    except Exception as e:\n",
    "                        print('Failed to delete %s. Reason: %s' % (file_path, e))\n",
    "\n",
    "                # Remove existing files in sample directory\n",
    "                for filename in os.listdir(save_dir_reduced):\n",
    "                    file_path = os.path.join(save_dir_reduced, filename)\n",
    "                    try:\n",
    "                        if os.path.isfile(file_path) or os.path.islink(file_path):\n",
    "                            os.unlink(file_path)\n",
    "                        elif os.path.isdir(file_path):\n",
    "                            shutil.rmtree(file_path)\n",
    "                    except Exception as e:\n",
    "                        print('Failed to delete %s. Reason: %s' % (file_path, e))\n",
    "\n",
    "                integrated_data = integ.integrateImageStack(scan)\n",
    "                integrated_data = integrated_data.unstack('system')\n",
    "\n",
    "                ## Make sacrificial image\n",
    "                print('Making sacrificial FY NEXAFS image')\n",
    "\n",
    "                plt.rcParams.update({\"font.size\": 18})\n",
    "                plt.rcParams.update({\"axes.linewidth\": 2})\n",
    "\n",
    "                fig = plt.figure(figsize = (3, 3))\n",
    "                ax = fig.add_axes([0, 0, 1, 1])\n",
    "\n",
    "                integrated_data.mean('chi').sel(q=slice(0.06,1)).mean('q').plot()\n",
    "                plt.gca().get_lines()[0].set_color('k')\n",
    "\n",
    "                plt.title(f'{integrated_data.sample_name} {current_edge} FY NEXAFS',fontsize = 17.5, pad = 15)\n",
    "                plt.ylabel('Intensity (a.u.)')\n",
    "                plt.xlabel('Energy (eV)')\n",
    "                plt.xlim([np.nanmin(integrated_data.energy), np.nanmax(integrated_data.energy)])\n",
    "\n",
    "                # Edit the major and minor ticks of the x and y axes\n",
    "                ax.xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', top = 'on')\n",
    "                ax.xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', top = 'on')\n",
    "                ax.yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "                ax.yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "                # plt.show()\n",
    "\n",
    "                current_vmax_detImage = []\n",
    "                current_vmax_cake = []\n",
    "                current_linecuts_ylim_min = []\n",
    "                current_linecuts_ylim_max = []\n",
    "\n",
    "                ### Loop for iterating through energies\n",
    "                for current_energy in np.sort(integrated_data.energy):\n",
    "                    clear_output(wait=True)\n",
    "                    print(f'Currently pulling data from \\'{base_path}\\' ({base_path_tracker}/{len(base_paths)})')\n",
    "                    print(f'Currently processing {scan_unstacked.sample_name} {current_edge} scan ({pol}° X-ray polarization) (scan {current_scan}/{np.nanmax(scans)}) (on scaled contrast images (1/2)) ({current_energy:.2f} eV/{np.nanmax(integrated_data.energy):.2f} eV)')\n",
    "\n",
    "                    ## Make plots of detector image\n",
    "\n",
    "                    plt.rcParams.update({\"axes.linewidth\": 2})\n",
    "                    plt.rcParams.update({\"font.size\": 18})\n",
    "\n",
    "                    fig = plt.figure(figsize = (3, 3))\n",
    "                    ax = fig.add_axes([0, 0, 1, 1])\n",
    "\n",
    "                    current_vmax_detImage.append(np.nanpercentile(scan_unstacked.sel(energy=current_energy), 97.5))\n",
    "\n",
    "                    im = ax.imshow(scan_unstacked.sel(energy=current_energy),\n",
    "                                    # norm = LogNorm(1,current_vmax_detImage[-1]),\n",
    "                                    vmin = 0,\n",
    "                                    vmax = current_vmax_detImage[-1],\n",
    "                                    cmap = 'viridis',\n",
    "                                    aspect = 1)\n",
    "\n",
    "                    plt.colorbar(im,\n",
    "                                 # format = ticker.FuncFormatter(fmt),\n",
    "                                 fraction = 0.046, pad = 0.04)\n",
    "                    plt.title(f'{integrated_data.sample_name} {current_edge}, {current_energy:.2f} eV, {pol:.0f}°', fontsize = 17.5, pad = 15)\n",
    "                    # plt.xlabel('$\\it{q_{x}}$ (Å$^{-1}$)')\n",
    "                    # plt.ylabel('$\\it{q_{y}}$ (Å$^{-1}$)')\n",
    "                    plt.xticks([])\n",
    "                    plt.yticks([])\n",
    "\n",
    "                    # plt.show()\n",
    "\n",
    "                    fname = f'{integrated_data.sample_name}_{current_energy}eV_{current_edge_folder}_{pol}pol_detImage_scaled'\n",
    "                    save_path = os.path.join(save_dir, fname + '.png')\n",
    "                    fig.savefig(save_path, bbox_inches = 'tight', pad_inches = 0.2, dpi = 150, transparent = False)\n",
    "\n",
    "                    ## Make linecuts plot for current energy\n",
    "\n",
    "                    plt.rcParams.update({\"font.size\": 18})\n",
    "                    plt.rcParams.update({\"axes.linewidth\": 2})\n",
    "\n",
    "                    fig = plt.figure(figsize = (3, 3))\n",
    "                    ax = fig.add_axes([0, 0, 1, 1])\n",
    "\n",
    "                    # integrated_data.rsoxs.slice_chi(-90,chi_width=10).sel(energy=current_energy,q=slice(0.01,1)).plot(label = '-90°')\n",
    "                    # integrated_data.rsoxs.slice_chi(0,chi_width=10).sel(energy=current_enepara_polrgy,q=slice(0.01,1)).plot(label = '0°')\n",
    "                    # integrated_data.rsoxs.slice_chi(90,chi_width=10).sel(energy=current_energy,q=slice(0.01,1)).plot(label = '90°')\n",
    "                    # integrated_data.rsoxs.slice_chi(180,chi_width=10).sel(energy=current_energy,q=slice(0.01,1)).plot(label = '±180°')\n",
    "\n",
    "                    # para_pol = integrated_data.rsoxs.slice_chi([90, -90],chi_width=10).sel(energy=current_energy,q=slice(0.01,1)) + integrated_data.rsoxs.slice_chi(90,chi_width=10).sel(energy=current_energy,q=slice(0.01,1))\n",
    "                    # perp_pol = integrated_data.rsoxs.slice_chi([0, 180],chi_width=10).sel(energy=current_energy,q=slice(0.01,1)) + integrated_data.rsoxs.slice_chi(180,chi_width=10).sel(energy=current_energy,q=slice(0.01,1))\n",
    "\n",
    "                    if pol == 0:\n",
    "                        para_pol = integrated_data.rsoxs.slice_chi(-90,chi_width=10).sel(energy=current_energy,q=slice(0.01,1))\n",
    "                        perp_pol = integrated_data.rsoxs.slice_chi(0,chi_width=10).sel(energy=current_energy,q=slice(0.01,1))\n",
    "                    elif pol == 90:\n",
    "                        perp_pol = integrated_data.rsoxs.slice_chi(-90,chi_width=10).sel(energy=current_energy,q=slice(0.01,1))\n",
    "                        para_pol = integrated_data.rsoxs.slice_chi(0,chi_width=10).sel(energy=current_energy,q=slice(0.01,1))   \n",
    "                    else:\n",
    "                        print('Error in finding polarization')\n",
    "\n",
    "                    plt.errorbar(np.asarray(para_pol.q).flatten(), np.asarray(para_pol).flatten(), np.sqrt(para_pol.values).flatten(),\n",
    "                                linestyle ='none', color ='k', alpha=0.25)\n",
    "                    plt.errorbar(np.asarray(perp_pol.q).flatten(), np.asarray(perp_pol).flatten(), np.sqrt(perp_pol.values).flatten(),\n",
    "                                linestyle ='none', color ='r', alpha=0.25)\n",
    "\n",
    "                    plt.plot(para_pol.q, para_pol, 'k-', label = 'Parallel')\n",
    "                    plt.plot(perp_pol.q, perp_pol, 'r-', label = 'Perpendicular')\n",
    "\n",
    "                    ax.legend(bbox_to_anchor=(1, 1), loc=1, frameon=False, fontsize=15)\n",
    "\n",
    "                    plt.title(f'{integrated_data.sample_name} {current_edge}, {current_energy:.2f} eV, {pol:.0f}°', fontsize = 17.5, pad = 15)\n",
    "                    plt.ylabel('Intensity (a.u.)')\n",
    "                    plt.xlabel('$\\it{q}$ (Å$^{-1}$)')\n",
    "                    plt.xlim([start_q, end_q])\n",
    "                    # ax.set_yscale('log')\n",
    "\n",
    "                    plt.xticks([0.01, 0.05, 0.1])\n",
    "                    plt.gca().yaxis.set_major_formatter(ticker.FuncFormatter(fmt))\n",
    "\n",
    "                    # Edit the major and minor ticks of the x and y axes\n",
    "                    ax.xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', top = 'on')\n",
    "                    ax.xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', top = 'on')\n",
    "                    ax.yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "                    ax.yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "                    # plt.show()\n",
    "\n",
    "                    current_linecuts_ylim_min.append(ax.get_ylim()[0])\n",
    "                    current_linecuts_ylim_max.append(ax.get_ylim()[1])\n",
    "\n",
    "                    fname = f'{integrated_data.sample_name}_{current_energy}eV_{current_edge_folder}_{pol}pol_linecuts_scaled'\n",
    "                    save_path = os.path.join(save_dir, fname + '.png')\n",
    "                    fig.savefig(save_path, bbox_inches = 'tight', pad_inches = 0.2, dpi = 150, transparent = False)\n",
    "\n",
    "                    ## Make cake plot for current energy\n",
    "\n",
    "                    plt.rcParams.update({\"font.size\": 18})\n",
    "                    plt.rcParams.update({\"axes.linewidth\": 2})\n",
    "\n",
    "                    fig = plt.figure(figsize = (3, 3))\n",
    "                    ax = fig.add_axes([0, 0, 1, 1])\n",
    "\n",
    "                    current_vmax_cake.append(np.nanpercentile(integrated_data.sel(energy=current_energy,q=slice(0.01,1)), 90))\n",
    "\n",
    "                    integrated_data.sel(energy=current_energy,q=slice(0.01,1)).plot(\n",
    "                        # norm = LogNorm(1,current_vmax_cake[-1]),\n",
    "                        vmin = 0,\n",
    "                        vmax = current_vmax_cake[-1],\n",
    "                        cmap = 'viridis')\n",
    "\n",
    "                    plt.title(f'{integrated_data.sample_name} {current_edge}, {current_energy:.2f} eV, {pol:.0f}°', fontsize = 17.5, pad = 15)\n",
    "                    plt.ylabel('Azimuth, $\\it{χ}$')\n",
    "                    plt.xlabel('$\\it{q}$ (Å$^{-1}$)')\n",
    "\n",
    "                    plt.yticks([-180, -90, 0, 90, 180],['-180°', '-90°', '0°', '90°', '180°'])\n",
    "                    plt.xticks([0.01, 0.05, 0.1])\n",
    "\n",
    "                    # Edit the major and minor ticks of the x and y axes\n",
    "                    ax.xaxis.set_tick_params(which='major', size=5, width=2, direction='inout', top='on')\n",
    "                    ax.xaxis.set_tick_params(which='minor', size=2.5, width=2, direction='inout', top='on')\n",
    "                    ax.yaxis.set_tick_params(which='major', size=5, width=2, direction='inout', right='on')\n",
    "                    ax.yaxis.set_tick_params(which='minor', size=2.5, width=2, direction='inout', right='on')\n",
    "\n",
    "                    # plt.show()\n",
    "\n",
    "                    fname = f'{integrated_data.sample_name}_{current_energy}eV_{current_edge_folder}_{pol}pol_cake_scaled'\n",
    "                    save_path = os.path.join(save_dir, fname + '.png')\n",
    "                    fig.savefig(save_path, bbox_inches = 'tight', pad_inches = 0.2, dpi = 150, transparent = False)\n",
    "\n",
    "                    ## Make NEXAFS with dot for current energy\n",
    "\n",
    "                    plt.rcParams.update({\"font.size\": 18})\n",
    "                    plt.rcParams.update({\"axes.linewidth\": 2})\n",
    "\n",
    "                    fig = plt.figure(figsize = (3, 3))\n",
    "                    ax = fig.add_axes([0, 0, 1, 1])\n",
    "\n",
    "                    integrated_data.mean('chi').sel(q=slice(0.06,1)).mean('q').plot()\n",
    "                    plt.gca().get_lines()[0].set_color('k')\n",
    "                    plt.plot(current_energy, integrated_data.mean('chi').sel(q=slice(0.06,1)).mean('q').sel(energy = current_energy), 'ro')\n",
    "\n",
    "                    plt.title(f'{integrated_data.sample_name} {current_edge} FY NEXAFS, {pol:.0f}°',fontsize = 17.5, pad = 15)\n",
    "                    plt.ylabel('Intensity (a.u.)')\n",
    "                    plt.xlabel('Energy (eV)')\n",
    "                    plt.xlim([np.nanmin(integrated_data.energy), np.nanmax(integrated_data.energy)])\n",
    "\n",
    "                    # Edit the major and minor ticks of the x and y axes\n",
    "                    ax.xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', top = 'on')\n",
    "                    ax.xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', top = 'on')\n",
    "                    ax.yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "                    ax.yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "                    # plt.show()\n",
    "\n",
    "                    fname = f'{integrated_data.sample_name}_{current_energy}eV_{current_edge_folder}_{pol}pol_FY_NEXAFS_en'\n",
    "                    save_path = os.path.join(save_dir, fname + '.png')\n",
    "                    fig.savefig(save_path, bbox_inches = 'tight', pad_inches = 0.2, dpi = 150, transparent = False)\n",
    "\n",
    "                ### Second loop for iterating through energies\n",
    "\n",
    "                vmax_detImage = np.nanmax(current_vmax_detImage)\n",
    "                vmax_cake = np.nanmax(current_vmax_cake)\n",
    "                linecuts_ylim_min = np.nanmin(current_linecuts_ylim_min)\n",
    "                linecuts_ylim_max = np.nanmax(current_linecuts_ylim_max)\n",
    "\n",
    "                for current_energy in np.sort(integrated_data.energy):\n",
    "                    clear_output(wait=True)\n",
    "                    print(f'Currently pulling data from \\'{base_path}\\' ({base_path_tracker}/{len(base_paths)})')\n",
    "                    print(f'Currently processing {scan_unstacked.sample_name} {current_edge} scan ({pol}° X-ray polarization) (scan {current_scan}/{np.nanmax(scans)}) (on fixed contrast images (2/2)) ({current_energy:.2f} eV/{np.nanmax(integrated_data.energy):.2f} eV)')\n",
    "\n",
    "                    ## Make plots of detector image\n",
    "\n",
    "                    plt.rcParams.update({\"axes.linewidth\": 2})\n",
    "                    plt.rcParams.update({\"font.size\": 18})\n",
    "\n",
    "                    fig = plt.figure(figsize = (3, 3))\n",
    "                    ax = fig.add_axes([0, 0, 1, 1])\n",
    "\n",
    "                    im = ax.imshow(scan_unstacked.sel(energy=current_energy),\n",
    "                                    # norm = LogNorm(1,vmax_detImage),\n",
    "                                    vmin = 0,\n",
    "                                    vmax = vmax_detImage,\n",
    "                                    cmap = 'viridis',\n",
    "                                    aspect = 1)\n",
    "\n",
    "                    plt.colorbar(im,\n",
    "                                 # format = ticker.FuncFormatter(fmt),\n",
    "                                 fraction = 0.046, pad = 0.04)\n",
    "                    plt.title(f'{integrated_data.sample_name} {current_edge}, {current_energy:.2f} eV, {pol:.0f}°', fontsize = 17.5, pad = 15)\n",
    "                    plt.xticks([])\n",
    "                    plt.yticks([])\n",
    "\n",
    "                    # plt.show()\n",
    "\n",
    "                    fname = f'{integrated_data.sample_name}_{current_energy}eV_{current_edge_folder}_{pol}pol_detImage_abs'\n",
    "                    save_path = os.path.join(save_dir, fname + '.png')\n",
    "                    fig.savefig(save_path, bbox_inches = 'tight', pad_inches = 0.2, dpi = 150, transparent = False)\n",
    "\n",
    "                    ## Make linecuts plot for current energy\n",
    "\n",
    "                    plt.rcParams.update({\"font.size\": 18})\n",
    "                    plt.rcParams.update({\"axes.linewidth\": 2})\n",
    "\n",
    "                    fig = plt.figure(figsize = (3, 3))\n",
    "                    ax = fig.add_axes([0, 0, 1, 1])\n",
    "\n",
    "                    if pol == 0:\n",
    "                        para_pol = integrated_data.rsoxs.slice_chi(-90,chi_width=10).sel(energy=current_energy,q=slice(0.01,1))\n",
    "                        perp_pol = integrated_data.rsoxs.slice_chi(0,chi_width=10).sel(energy=current_energy,q=slice(0.01,1))\n",
    "                    elif pol == 90:\n",
    "                        perp_pol = integrated_data.rsoxs.slice_chi(-90,chi_width=10).sel(energy=current_energy,q=slice(0.01,1))\n",
    "                        para_pol = integrated_data.rsoxs.slice_chi(0,chi_width=10).sel(energy=current_energy,q=slice(0.01,1))   \n",
    "                    else:\n",
    "                        print('Error in finding polarization')\n",
    "\n",
    "                    plt.errorbar(np.asarray(para_pol.q).flatten(), np.asarray(para_pol).flatten(), np.sqrt(para_pol.values).flatten(),\n",
    "                                linestyle ='none', color ='k', alpha=0.25)\n",
    "                    plt.errorbar(np.asarray(perp_pol.q).flatten(), np.asarray(perp_pol).flatten(), np.sqrt(perp_pol.values).flatten(),\n",
    "                                linestyle ='none', color ='r', alpha=0.25)\n",
    "\n",
    "                    plt.plot(para_pol.q, para_pol, 'k-', label = 'Parallel')\n",
    "                    plt.plot(perp_pol.q, perp_pol, 'r-', label = 'Perpendicular')\n",
    "\n",
    "                    ax.legend(bbox_to_anchor=(1, 1), loc=1, frameon=False, fontsize=15)\n",
    "\n",
    "                    plt.title(f'{integrated_data.sample_name} {current_edge}, {current_energy:.2f} eV, {pol:.0f}°', fontsize = 17.5, pad = 15)\n",
    "                    plt.ylabel('Intensity (a.u.)')\n",
    "                    plt.xlabel('$\\it{q}$ (Å$^{-1}$)')\n",
    "                    plt.xlim([start_q, end_q])\n",
    "                    plt.ylim([linecuts_ylim_min, linecuts_ylim_max])\n",
    "                    # ax.set_yscale('log')\n",
    "                    plt.xticks([0.01, 0.05, 0.1])\n",
    "                    plt.gca().yaxis.set_major_formatter(ticker.FuncFormatter(fmt))\n",
    "\n",
    "                    # Edit the major and minor ticks of the x and y axes\n",
    "                    ax.xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', top = 'on')\n",
    "                    ax.xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', top = 'on')\n",
    "                    ax.yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "                    ax.yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "                    # plt.show()\n",
    "\n",
    "                    fname = f'{integrated_data.sample_name}_{current_energy}eV_{current_edge_folder}_{pol}pol_linecuts_abs'\n",
    "                    save_path = os.path.join(save_dir, fname + '.png')\n",
    "                    fig.savefig(save_path, bbox_inches = 'tight', pad_inches = 0.2, dpi = 150, transparent = False)\n",
    "\n",
    "                    ## Make cake plot for current energy\n",
    "\n",
    "                    plt.rcParams.update({\"font.size\": 18})\n",
    "                    plt.rcParams.update({\"axes.linewidth\": 2})\n",
    "\n",
    "                    fig = plt.figure(figsize = (3, 3))\n",
    "                    ax = fig.add_axes([0, 0, 1, 1])\n",
    "\n",
    "                    integrated_data.sel(energy=current_energy,q=slice(0.01,1)).plot(\n",
    "                        # norm = LogNorm(1,vmax_cake),\n",
    "                        vmin = 0,\n",
    "                        vmax = vmax_cake,\n",
    "                        cmap = 'viridis')\n",
    "\n",
    "                    plt.title(f'{integrated_data.sample_name} {current_edge}, {current_energy:.2f} eV, {pol:.0f}°', fontsize = 17.5, pad = 15)\n",
    "                    plt.ylabel('Azimuth, $\\it{χ}$')\n",
    "                    plt.xlabel('$\\it{q}$ (Å$^{-1}$)')\n",
    "\n",
    "                    plt.yticks([-180, -90, 0, 90, 180],['-180°', '-90°', '0°', '90°', '180°'])\n",
    "                    plt.xticks([0.01, 0.05, 0.1])\n",
    "\n",
    "                    # Edit the major and minor ticks of the x and y axes\n",
    "                    ax.xaxis.set_tick_params(which='major', size=5, width=2, direction='inout', top='on')\n",
    "                    ax.xaxis.set_tick_params(which='minor', size=2.5, width=2, direction='inout', top='on')\n",
    "                    ax.yaxis.set_tick_params(which='major', size=5, width=2, direction='inout', right='on')\n",
    "                    ax.yaxis.set_tick_params(which='minor', size=2.5, width=2, direction='inout', right='on')\n",
    "\n",
    "                    # plt.show()\n",
    "\n",
    "                    fname = f'{integrated_data.sample_name}_{current_energy}eV_{current_edge_folder}_{pol}pol_cake_abs'\n",
    "                    save_path = os.path.join(save_dir, fname + '.png')\n",
    "                    fig.savefig(save_path, bbox_inches = 'tight', pad_inches = 0.2, dpi = 150, transparent = False)\n",
    "\n",
    "                # GIF making part\n",
    "                image_folder = os.fsencode(save_dir)\n",
    "\n",
    "                ## make gif of fluorescence yield NEXAFS plots\n",
    "                print('Making FY NEXAFS gif (1/7)')\n",
    "                filenames = []\n",
    "\n",
    "                for file in os.listdir(image_folder):\n",
    "                    filename = os.fsdecode(file)\n",
    "                    if filename.endswith( ('_FY_NEXAFS_en.png') ):\n",
    "                        filenames.append(filename)\n",
    "\n",
    "                filenames.sort() # this iteration technique has no built in order, so sort the frames\n",
    "\n",
    "                # Make it pause at the end so that the viewers can ponder\n",
    "                for _ in range(10):\n",
    "                    filenames.append(filenames[-1])\n",
    "\n",
    "                images = list(map(lambda filename: imageio.imread(os.path.join(save_dir, filename)), filenames))\n",
    "\n",
    "                imageio.mimsave(os.path.join(save_dir_reduced, sample_scan + '_FY_NEXAFS_en_movie.gif'), images, duration = 0.3) # modify duration as needed\n",
    "\n",
    "                ## make gif of scaled cake images\n",
    "                print('Making scaled cake gif (2/7)')\n",
    "                filenames = []\n",
    "\n",
    "                for file in os.listdir(save_dir):\n",
    "                    filename = os.fsdecode(file)\n",
    "                    if filename.endswith( ('_cake_scaled.png') ):\n",
    "                        filenames.append(filename)\n",
    "\n",
    "                filenames.sort() # this iteration technique has no built in order, so sort the frames\n",
    "\n",
    "                # Make it pause at the end so that the viewers can ponder\n",
    "                for _ in range(10):\n",
    "                    filenames.append(filenames[-1])\n",
    "\n",
    "                images = list(map(lambda filename: imageio.imread(os.path.join(save_dir, filename)), filenames))\n",
    "\n",
    "                imageio.mimsave(os.path.join(save_dir_reduced, sample_scan + '_cake_scaled_movie.gif'), images, duration = 0.3) # modify duration as needed\n",
    "\n",
    "                ## make gif of abs cake images\n",
    "                print('Making abs cake gif (3/7)')\n",
    "                filenames = []\n",
    "\n",
    "                for file in os.listdir(save_dir):\n",
    "                    filename = os.fsdecode(file)\n",
    "                    if filename.endswith( ('_cake_abs.png') ):\n",
    "                        filenames.append(filename)\n",
    "\n",
    "                filenames.sort() # this iteration technique has no built in order, so sort the frames\n",
    "\n",
    "                # Make it pause at the end so that the viewers can ponder\n",
    "                for _ in range(10):\n",
    "                    filenames.append(filenames[-1])\n",
    "\n",
    "                images = list(map(lambda filename: imageio.imread(os.path.join(save_dir, filename)), filenames))\n",
    "\n",
    "                imageio.mimsave(os.path.join(save_dir_reduced, sample_scan + '_cake_abs_movie.gif'), images, duration = 0.3) # modify duration as needed\n",
    "\n",
    "                ## make gif of scaled detector image\n",
    "                print('Making scaled detector image gif (4/7)')\n",
    "                filenames = []\n",
    "\n",
    "                for file in os.listdir(image_folder):\n",
    "                    filename = os.fsdecode(file)\n",
    "                    if filename.endswith( ('detImage_scaled.png') ):\n",
    "                        filenames.append(filename)\n",
    "\n",
    "                filenames.sort() # this iteration technique has no built in order, so sort the frames\n",
    "\n",
    "                # Make it pause at the end so that the viewers can ponder\n",
    "                for _ in range(10):\n",
    "                    filenames.append(filenames[-1])\n",
    "\n",
    "                images = list(map(lambda filename: imageio.imread(os.path.join(save_dir, filename)), filenames))\n",
    "\n",
    "                imageio.mimsave(os.path.join(save_dir_reduced, sample_scan + '_detImage_scaled_movie.gif'), images, duration = 0.3) # modify duration as needed\n",
    "\n",
    "                ## make gif of abs detector image\n",
    "                print('Making abs detector image gif (5/7)')\n",
    "                filenames = []\n",
    "\n",
    "                for file in os.listdir(image_folder):\n",
    "                    filename = os.fsdecode(file)\n",
    "                    if filename.endswith( ('detImage_abs.png') ):\n",
    "                        filenames.append(filename)\n",
    "\n",
    "                filenames.sort() # this iteration technique has no built in order, so sort the frames\n",
    "\n",
    "                # Make it pause at the end so that the viewers can ponder\n",
    "                for _ in range(10):\n",
    "                    filenames.append(filenames[-1])\n",
    "\n",
    "                images = list(map(lambda filename: imageio.imread(os.path.join(save_dir, filename)), filenames))\n",
    "\n",
    "                imageio.mimsave(os.path.join(save_dir_reduced, sample_scan + '_detImage_abs_movie.gif'), images, duration = 0.3) # modify duration as needed\n",
    "\n",
    "                ## make gif of scaled linecuts plots\n",
    "                print('Making scaled linecuts gif (6/7)')\n",
    "                filenames = []\n",
    "\n",
    "                for file in os.listdir(image_folder):\n",
    "                    filename = os.fsdecode(file)\n",
    "                    if filename.endswith( ('_linecuts_scaled.png') ):\n",
    "                        filenames.append(filename)\n",
    "\n",
    "                filenames.sort() # this iteration technique has no built in order, so sort the frames\n",
    "\n",
    "                # Make it pause at the end so that the viewers can ponder\n",
    "                for _ in range(10):\n",
    "                    filenames.append(filenames[-1])\n",
    "\n",
    "                images = list(map(lambda filename: imageio.imread(os.path.join(save_dir, filename)), filenames))\n",
    "\n",
    "                imageio.mimsave(os.path.join(save_dir_reduced, sample_scan + '_linecuts_scaled_movie.gif'), images, duration = 0.3) # modify duration as needed\n",
    "\n",
    "                ## make gif of abs linecuts plots\n",
    "                print('Making abs linecuts gif (7/7)')\n",
    "                filenames = []\n",
    "\n",
    "                for file in os.listdir(image_folder):\n",
    "                    filename = os.fsdecode(file)\n",
    "                    if filename.endswith( ('_linecuts_abs.png') ):\n",
    "                        filenames.append(filename)\n",
    "\n",
    "                filenames.sort() # this iteration technique has no built in order, so sort the frames\n",
    "\n",
    "                # Make it pause at the end so that the viewers can ponder\n",
    "                for _ in range(10):\n",
    "                    filenames.append(filenames[-1])\n",
    "\n",
    "                images = list(map(lambda filename: imageio.imread(os.path.join(save_dir, filename)), filenames))\n",
    "\n",
    "                imageio.mimsave(os.path.join(save_dir_reduced, sample_scan + '_linecuts_abs_movie.gif'), images, duration = 0.3) # modify duration as needed\n",
    "\n",
    "                ## End of gif making part\n",
    "                ## Make fluoresence yield NEXAFS for scan across all energies\n",
    "                print('Making FY NEXAFS image (1/3)')\n",
    "\n",
    "                plt.rcParams.update({\"font.size\": 18})\n",
    "                plt.rcParams.update({\"axes.linewidth\": 2})\n",
    "\n",
    "                fig = plt.figure(figsize = (3, 3))\n",
    "                ax = fig.add_axes([0, 0, 1, 1])\n",
    "\n",
    "                integrated_data.mean('chi').sel(q=slice(0.06,1)).mean('q').plot()\n",
    "                plt.gca().get_lines()[0].set_color('k')\n",
    "\n",
    "                plt.title(f'{integrated_data.sample_name} {current_edge} FY NEXAFS',fontsize = 17.5, pad = 15)\n",
    "                plt.ylabel('Intensity (a.u.)')\n",
    "                plt.xlabel('Energy (eV)')\n",
    "                plt.xlim([np.nanmin(integrated_data.energy), np.nanmax(integrated_data.energy)])\n",
    "\n",
    "                # Edit the major and minor ticks of the x and y axes\n",
    "                ax.xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', top = 'on')\n",
    "                ax.xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', top = 'on')\n",
    "                ax.yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "                ax.yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "                # plt.show()\n",
    "\n",
    "                fname = f'{integrated_data.sample_name}_{current_edge_folder}_pol{pol}_FY_NEXAFS'\n",
    "                save_path = os.path.join(save_dir_reduced, fname + '.png')\n",
    "                fig.savefig(save_path, bbox_inches = 'tight', pad_inches = 0.2, dpi=300, transparent=False)\n",
    "\n",
    "                ## Make intensity vs q for scan across all energies\n",
    "                print('Making intensity vs q image (2/3)')\n",
    "\n",
    "                plt.rcParams.update({\"font.size\": 18})\n",
    "                plt.rcParams.update({\"axes.linewidth\": 2})\n",
    "\n",
    "                fig = plt.figure(figsize = (3, 3))\n",
    "                ax = fig.add_axes([0, 0, 1, 1])\n",
    "\n",
    "                integrated_data.mean('chi').sel(q=slice(0.01,1)).plot(\n",
    "                    norm = LogNorm(1, np.nanpercentile(integrated_data.mean('chi'), 97.5)),\n",
    "                    # vmin = 0,\n",
    "                    # vmax = np.nanpercentile(integrated_data.mean('chi'), 97.5),\n",
    "                    cmap = 'viridis')\n",
    "\n",
    "                plt.title(f'{integrated_data.sample_name} {current_edge}',fontsize = 17.5, pad = 15)\n",
    "                plt.ylabel('$\\it{q}$ (Å$^{-1}$)')\n",
    "                plt.xlabel('Energy (eV)')\n",
    "                # plt.yticks([0.01, 0.02, 0.04, 0.06, 0.08, 0.1])\n",
    "\n",
    "                # Edit the major and minor ticks of the x and y axes\n",
    "                ax.xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', top = 'on')\n",
    "                ax.xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', top = 'on')\n",
    "                ax.yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "                ax.yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "                # plt.show()\n",
    "\n",
    "                fname = f'{integrated_data.sample_name}_{current_edge_folder}_pol{pol}_IqE'\n",
    "                save_path = os.path.join(save_dir_reduced, fname + '.png')\n",
    "                fig.savefig(save_path, bbox_inches = 'tight', pad_inches = 0.2, dpi=300, transparent=False)\n",
    "\n",
    "                ## Make anisotropy ratio plot for scan across all energies\n",
    "                print('Making anisotropy ratio image (3/3)')\n",
    "\n",
    "                plt.rcParams.update({\"font.size\": 18})\n",
    "                plt.rcParams.update({\"axes.linewidth\": 2})\n",
    "\n",
    "                fig = plt.figure(figsize = (3, 3))\n",
    "                ax = fig.add_axes([0, 0, 1, 1])\n",
    "\n",
    "                integrated_data.rsoxs.AR(chi_width=10).sel(q=slice(0.01,1)).plot(\n",
    "                    vmin = -1,\n",
    "                    vmax = 1,\n",
    "                    cmap = 'RdBu')\n",
    "\n",
    "                plt.title(f'{integrated_data.sample_name} {current_edge}, {pol:.0f}° AR',fontsize = 17.5, pad = 15)\n",
    "                plt.xlabel('Energy (eV)')\n",
    "                plt.ylabel('$\\it{q}$ (Å$^{-1}$)')\n",
    "                # plt.yticks([0.01, 0.02, 0.04, 0.06, 0.08, 0.1])\n",
    "\n",
    "                # Edit the major and minor ticks of the x and y axes\n",
    "                ax.xaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', top = 'on')\n",
    "                ax.xaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', top = 'on')\n",
    "                ax.yaxis.set_tick_params(which = 'major', size = 5, width = 2, direction = 'in', right = 'on')\n",
    "                ax.yaxis.set_tick_params(which = 'minor', size = 2.5, width = 2, direction = 'in', right = 'on')\n",
    "\n",
    "                # plt.show()\n",
    "\n",
    "                fname = f'reduced_pol{pol}_{current_edge_folder}_{scan_unstacked.sample_name}_AR'\n",
    "                save_path = os.path.join(save_dir_reduced, fname + '.png')\n",
    "                fig.savefig(save_path, bbox_inches = 'tight', pad_inches = 0.2, dpi=300, transparent=False)\n",
    "            else:\n",
    "                print(f'Skipping scan {current_scan} - already fully processed')\n",
    "        else:\n",
    "            print(f'Skipping scan {current_scan} - ID not in range')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}