xps_processing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b4fe2953-8c13-4df1-95ec-bb37484f9012",
   "metadata": {},
   "source": [
    "# Case Study - X-ray Photoelectron Spectroscopy (XPS)\n",
    "\n",
    "<div class=\"overview admonition\"> \n",
    "<p class=\"admonition-title\">Overview</p>\n",
    "\n",
    "Questions:\n",
    "\n",
    "* How can I use Python to fit analyze XPS data?\n",
    "\n",
    "Objectives:\n",
    "\n",
    "* Use pandas to read XPS data.\n",
    "* Use SciPy to find peaks.\n",
    "* Create interactive plots data.\n",
    "* Perform a background subtraction.\n",
    "* Find and  fit peaks with Gaussian functions.\n",
    "\n",
    "</div>\n",
    "\n",
    "\n",
    "XPS is a surface-sensitive quantitative technique that measures the elemental composition, empirical formula, chemical state, and electronic state of elements present in a material. \n",
    "It works on the principle of photoemission, where an X-ray beam ejects core-level electrons from the surface of the sample, and the kinetic energy of these electrons is measured.\n",
    "\n",
    "In this notebook, we will process [XPS spectra of silicon](https://pubs.aip.org/avs/sss/article-abstract/20/1/36/366432/Silicon-100-SiO2-by-XPS?redirectedFrom=fulltext).\n",
    "We will use many of the skills learned in previous notebooks including processing data using pandas and visualization using plotly.\n",
    "This notebook introduces a new library - SciPy. \n",
    "SciPy is an advanced Python library tailored for scientific computing and analysis.\n",
    "We will use SciPy for peak fitting and integration.\n",
    "\n",
    "Note that this notebook is **not** a best practices for XPS analysis - rather it shows you Python strategies that can be used for data processing and analysis.\n",
    "\n",
    "For this notebook, we will be using `pandas`, `plotly`, `scipy`, and `numpy`.\n",
    "\n",
    "We will follow these steps:\n",
    "\n",
    "1. Analyze wide scan spectrum - We will read in the wide scan and use SciPy `find_peaks` to find peaks.\n",
    "2. Analyze narrow spectra - For the oxygen and silicon peaks, we will do a more in-depth analysis.\n",
    "\n",
    "For the narrow spectra, we will:\n",
    "1. Load in the data.\n",
    "2. Perform a background correction using a polynomial fit to the background.\n",
    "3. Fit peaks in the spectrum using a Gaussian function.\n",
    "4. Integrate each peak."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d498013-ada8-4c60-accf-77b78a477692",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import plotly.express as px\n",
    "\n",
    "from scipy.signal import savgol_filter, find_peaks\n",
    "from scipy.optimize import curve_fit\n",
    "from scipy.integrate import simpson, quad"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c237ede4-d364-4408-a16c-c2a4e9be0c96",
   "metadata": {},
   "source": [
    "##  Reading in Data\n",
    "\n",
    "In this notebook, we will again use `pandas` to read in our data.\n",
    "For past notebooks, we have used the `read_csv` function. \n",
    "However, the data set we will use for this notebook is in a tab delimited file.\n",
    "\n",
    "We could use `read_csv` with this data, and specify that the delimiter is a tab character.\n",
    "However the [`read_table`](https://pandas.pydata.org/docs/reference/api/pandas.read_table.html) function is more general for this purpose.\n",
    "\n",
    "When you examine the data for this notebook by viewing one of the files, you will see something like the following:\n",
    "\n",
    "```\n",
    "1.3500000e+003\t6.5014286e+003\n",
    "1.3495000e+003\t6.4672997e+003\n",
    "1.3490000e+003\t6.4709445e+003\n",
    "1.3485000e+003\t6.4449666e+003\n",
    "```\n",
    "\n",
    "There are no headers for this data, so we will tell `pandas`that by putting in the argument `header=None`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb8a175a-7e35-45fe-877c-c300b1b563d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "wide_scan = pd.read_table(\"si_data/N0125701.asc\", header=None)\n",
    "wide_scan.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a9017da-2e01-4a4d-961e-175fb0e8a76c",
   "metadata": {},
   "source": [
    "Because there is no header, our columns do not yet have names.\n",
    "We will rename the columns to `energy` and `intensity`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2e02207-00d7-4377-92ab-c0e8cc861496",
   "metadata": {},
   "outputs": [],
   "source": [
    "wide_scan.columns = [\"energy\", \"intensity\"]\n",
    "wide_scan.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2cc35b0-a9c9-4004-ac14-44a81b7463ed",
   "metadata": {},
   "source": [
    "## Visualizing the Data\n",
    "\n",
    "Upon loading the data initially, you would probably wish to examine it.\n",
    "We will plot our data using plotly express and `px.line`.\n",
    "\n",
    "If you look at the order of the energy data in the `wide_scan` dataframe above, you will see that the energy\n",
    "is listed from highest to lowest. Typically XPS spectra are shown with the highest enregy on the left side of the axis.\n",
    "In order to have our graph appear this way, we can use `fig.update_xaxes(autorange=\"reversed\")`\n",
    "\n",
    "\n",
    "<div class=\"note admonition\"> \n",
    "<p class=\"admonition-title\">Why do we reverse the x axis?</p>\n",
    "\n",
    "We reverse the x axis on XPS spectra so that the order of the peaks represent how tightly bound electrons are.\n",
    "When the axis is reversed, more tightly bound electrons are on the left, and more loosely bound electrons are on the right.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b7d1785-5f05-4e68-93d0-20e1e983ceb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an interactive plot of the wide scan\n",
    "fig = px.line(wide_scan, x=\"energy\", y=\"intensity\") # line plot of data\n",
    "fig.update_xaxes(autorange=\"reversed\")  # reverse x axis\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "947f7bd7-a42d-4824-bc6a-e31da5c10319",
   "metadata": {},
   "source": [
    "## Data Smoothing\n",
    "\n",
    "Sometimes data may be noisy and require smoothing. \n",
    "The data set we are working with does not require smoothing, but an example is shown below in case you need to smooth data!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43b9da87-0a1c-4e29-9c95-fefb52709393",
   "metadata": {},
   "outputs": [],
   "source": [
    "smoothed = pd.DataFrame(savgol_filter(wide_scan, window_length=10, polyorder=3, axis=0))\n",
    "smoothed.columns = [\"energy\", \"intensity\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adf50fac-bc79-4b00-ac04-090d625f352d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an interactive plot of the wide scan\n",
    "fig = px.line(smoothed, x=\"energy\", y=\"intensity\")\n",
    "fig.update_xaxes(autorange=\"reversed\")  \n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b81125d-e837-4b97-8fd4-57daa2041260",
   "metadata": {},
   "source": [
    "## Peak Finding\n",
    "\n",
    "On our XPS spectrum, the position of peaks tells us what elements are present in the sample. \n",
    "In order to analyze our spectrum, we would want to find the peaks, then compare our peaks to peaks of known elements.\n",
    "\n",
    "SciPy provides a [function called `find_peaks`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html) that we will use to automatically pick peaks from our spectrum. \n",
    "A screenshot of the documentation `find_peaks` function is shown below. Note, that if you click the link above you will see a more in-depth explanation.\n",
    "\n",
    "![find_peaks](images/scipy_find_peaks.png)\n",
    "\n",
    "We will use the `prominence` keyword to define how much a peak must stand out from its neighbors. This will prevent us from picking peaks based on noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e50a87d1-6233-4723-ab56-e3ec287e4be0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit the peaks of the wide scan\n",
    "peak_index, properties = find_peaks(x=wide_scan[\"intensity\"], prominence=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da493d23-e7bb-40e9-ac64-102f365ebf38",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(peak_index)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08fc064e-5a5f-4cf5-8c0e-a7a80e85db94",
   "metadata": {},
   "source": [
    "SciPy has returned the *index numbers* for the peaks, rather than information about the peaks themselves.\n",
    "To understand what this means, consider the image below.\n",
    "If we were trying to find the peak for the list of numbers below, SciPy would return `2` to us - the index number of the maximum.\n",
    "\n",
    "![peak_index](images/peak_index.png)\n",
    "\n",
    "In our example above, if we wanted to get the actual value of the peak, we would need to index into our list.\n",
    "We will use `.iloc` along with the `peak_index` returned from SciPy to get the locations of our peaks from our dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "381252ce-0d03-4f72-b6d8-74237bc68513",
   "metadata": {},
   "outputs": [],
   "source": [
    "peaks = wide_scan.iloc[peak_index].copy()\n",
    "peaks.reset_index(inplace=True, drop=True)\n",
    "peaks.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9ab473b-f1fc-4b86-b0a0-55b549a96bfd",
   "metadata": {},
   "source": [
    "We can now visualize our peak locations using plotly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06ab9c57-933b-4335-af21-e3509fb61f56",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.line(wide_scan, x=\"energy\", y=\"intensity\")\n",
    "fig.add_scatter(x=peaks[\"energy\"], y=peaks[\"intensity\"], \n",
    "                mode=\"markers\", \n",
    "                name=\"peaks\")\n",
    "fig.update_xaxes(autorange=\"reversed\")   \n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c07f6bc7-a9e2-4658-8d7b-5360b8055c3d",
   "metadata": {},
   "source": [
    "<div class=\"exercise admonition\">\n",
    "<p class=\"admonition-title\">Exercise</p>\n",
    "\n",
    "As you can see, just using the `find_peaks` function returns many more peaks to us than we would choose to label. \n",
    "\n",
    "[Read the documentation for the `find_peaks` function.](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html)\n",
    "\n",
    "Try different arguments in order to find appropriate peaks in the spectrum.\n",
    "\n",
    "\n",
    "**To continue with the tutorial, change your `find_peaks` function above to**:\n",
    "\n",
    "```\n",
    "peak_index, properties = find_peaks(x=wide_scan[\"intensity\"], prominence=1000)\n",
    "```\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b70cb72f-4684-44f1-8de6-ec59d4e790e0",
   "metadata": {},
   "source": [
    "## Peak Labeling\n",
    "\n",
    "We might wish to add another column to our table that gives the peak identities.\n",
    "We can first create an empty column in our data frame for peak identities.\n",
    "We will start by creating a blank column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07df48a8-0368-4a7e-a52a-4ac7f96515a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a blanke column called \"labels\" with nothing in it.\n",
    "peaks[\"labels\"] = None\n",
    "\n",
    "# Preview the first 20 rows of the dataframe.\n",
    "peaks.head(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc7830ef-d890-49ea-90fb-6887af47bca9",
   "metadata": {},
   "source": [
    "Next, we will fill in identities for some of the peaks. \n",
    "For example, from looking up a [reference table](https://www.thermofisher.com/us/en/home/materials-science/learning-center/periodic-table.html), we would know that the peak at 100 is $Si$, and the peak at 104 is $Si$ in $SiO_2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "757a36ca-d58e-4225-bb21-259a6d0cbbe5",
   "metadata": {},
   "outputs": [],
   "source": [
    "peaks.loc[17, \"labels\"] = \"Si\"\n",
    "peaks.loc[16, \"labels\"] = \"Si (SiO2)\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "331c602a-d22e-4537-ba6e-7083e7976f10",
   "metadata": {},
   "outputs": [],
   "source": [
    "peaks.head(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e9d09e4-8282-42fa-bb91-13c38303e9fc",
   "metadata": {},
   "source": [
    "<div class=\"exercise admonition\">\n",
    "<p class=\"admonition-title\">Exercise</p>\n",
    "\n",
    "Use the following [reference table](https://www.thermofisher.com/us/en/home/materials-science/learning-center/periodic-table.html) to identify the oxygen peak.\n",
    "\n",
    "After you have identified it, add a label to the appropriate row.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edd39b5d-1e19-42b6-a230-9e6429f22514",
   "metadata": {},
   "source": [
    "We can add labels to our graphs by adding a `text` and `textposition` argument to our plot.\n",
    "\n",
    "Use the cell below to visualize your plot with labels.\n",
    "\n",
    "You may need to change the range of the y-axis once you have added the oxygen plot. \n",
    "You can do so with the command\n",
    "\n",
    "```python\n",
    "fig.update_yaxes(range=[minimum_y_value, maximum_y_value])\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "613b9852-2a64-4afe-b95c-7da4bda95a44",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = px.line(wide_scan, x=\"energy\", y=\"intensity\")\n",
    "fig.add_scatter(x=peaks[\"energy\"], y=peaks[\"intensity\"], \n",
    "                mode=\"markers+text\", \n",
    "                name=\"peaks\",\n",
    "                text=peaks[\"labels\"],  # this was added  to add labels\n",
    "                textposition=\"top center\",  # this was added to position the labels.\n",
    "               )\n",
    "fig.update_xaxes(autorange=\"reversed\")   \n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "361687d4-e0d6-4f23-af33-0831a4b65a9c",
   "metadata": {},
   "source": [
    "## Peak Analysis\n",
    "\n",
    "A common practice in XPS analysis is to analyze particular peaks in order to determine the chemical states, concentrations, and elemental composition of a sample.\n",
    "For the next section of this notebook, we will consider a close-up view of the Si peaks. \n",
    "We will fit and subtract the background signal, then we will fit Gaussian functions to the peaks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16c2e2b6-384c-4363-966e-a148ddc5d49d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# data loading\n",
    "si_peaks = pd.read_table(\"si_data/N0125702.asc\", header=None)\n",
    "\n",
    "# add column names\n",
    "si_peaks.columns = [\"energy\", \"intensity\"]\n",
    "\n",
    "# visualize\n",
    "si_fig = px.line(si_peaks, x=\"energy\", y=\"intensity\")\n",
    "si_fig.update_xaxes(autorange=\"reversed\")\n",
    "si_fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96c43097-5385-4e49-bf04-22b3d56518a9",
   "metadata": {},
   "source": [
    "### Background Subtraction\n",
    "\n",
    "Upon viewing the figure above, we can see the presence of some background signal, particularly on the left side of the graph.\n",
    "Using your mouse hovering over the graph to examine, you might decide that for the peak on the left the background corresponds to\n",
    "energies less than 98.5 eV or greater than 106.5 eV.\n",
    "\n",
    "In the cell below, we select parts of our signal that fit that criteria using\n",
    "\n",
    "```\n",
    "(98.5 > si_peaks[\"energy\"]) | (si_peaks[\"energy\"] > 106.5)\n",
    "```\n",
    "\n",
    "in this syntax, `|` represents \"or\". The phrase above says that the energy should be less than 98.5 or greater than 106.5. \n",
    "Similar to peak finding, this returns the index for where the statement is true.\n",
    "We get the actual values by using the index with the peaks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4d47523-ae78-4a2b-8a48-de475d1b984c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define what we want to fit for background subtraction\n",
    "background_index = (98.5 > si_peaks[\"energy\"]) | (si_peaks[\"energy\"] > 106.5)\n",
    "background = si_peaks[background_index]\n",
    "background.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a90382b-5b88-4592-85bc-e5d215d8b6ed",
   "metadata": {},
   "source": [
    "We can visualize what we've selected for the background by adding more data to the figure we created previously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68300d85-fc8a-40f5-b929-c74c9b896bf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "si_fig.add_scatter(x=background[\"energy\"], y=background[\"intensity\"], mode=\"markers\", name=\"background\")\n",
    "si_fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9226f432-3e4e-46cf-845d-87cd78e5399f",
   "metadata": {},
   "source": [
    "In order to subtract out this background, we will fit a polynomial function to our selected data then subtract it from the signal.\n",
    "This will adjust our baseline to 0.\n",
    "\n",
    "To do this in Python, we will use a [function from numpy called `polyfit`](https://numpy.org/doc/stable/reference/generated/numpy.polyfit.html). `polyfit` is specifically for fitting polynomials of any order.\n",
    "To use polyfit, you pass in your x values, your y values, and your desired polynomial degree.\n",
    "For this particular case, it looks as though a linear fit (polynomial degree=1) will be sufficient. \n",
    "However, you could increase this or change it if you needed a different polynomial degree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50b9daa0-371a-497a-b4d7-bc12723d0bbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit a polynomial of degree 1\n",
    "fit = np.polyfit(background[\"energy\"], background[\"intensity\"], 1)\n",
    "\n",
    "print(fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b2b6a23-9f61-45fb-b182-b84bc6ba28f1",
   "metadata": {},
   "source": [
    "What is returned to us in fit is the parameters for the fit. \n",
    "In this case, the first value is our slope and the second is the intercept.\n",
    "We can pass this into another function `np.polyval` to evaluate the value of the polynomial at given energies.\n",
    "We will simply evaluate the fit over the whole energy range for these peaks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55694bb7-4eac-452a-9aa6-1f518a2179b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get fit values\n",
    "vals = np.polyval(fit, si_peaks[\"energy\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f67339c4-9bc6-434d-8fe6-04c082fb2a9c",
   "metadata": {},
   "source": [
    "We can add this to our visualization to evaluate our fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b132b836-3362-4844-9caf-6c49a6be371b",
   "metadata": {},
   "outputs": [],
   "source": [
    "si_fig.add_scatter(x=si_peaks[\"energy\"], y=vals, name=\"background fit\")\n",
    "si_fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d2f9cea-eb97-4507-aea2-65ff68133f61",
   "metadata": {},
   "source": [
    "Now that we have determined the equation for our background, we can subtract it from our intensity to get our corrected peaks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa34c3e3-0bae-4436-9d40-1fe66273ed72",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Subtract the background\n",
    "si_peaks[\"corrected_intensity\"] = si_peaks[\"intensity\"] - vals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9033287e-1047-458b-a503-608cc347fa02",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig_corrected = px.line(si_peaks, x=\"energy\", y=\"corrected_intensity\")    \n",
    "fig_corrected.update_xaxes(autorange=\"reversed\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e68fbde1-b698-4cc4-b92d-5b8939027b37",
   "metadata": {},
   "source": [
    "<div class=\"exercise admonition\">\n",
    "<p class=\"admonition-title\">Exercise</p>\n",
    "\n",
    "Read in data for the oxygen peak and perform a background correction.\n",
    "\n",
    "The data for the oxygen peak is in the file `si_data/N0125704.asc`.\n",
    "\n",
    "Be sure to save oxygen data in appropriate variable names. \n",
    "You will use both O and Si data later.\n",
    "\n",
    "</div>\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "241e5c0b-65ae-490b-a2df-01a09c0367b3",
   "metadata": {},
   "source": [
    "### Numerical Integration of Peaks\n",
    "\n",
    "Although it is possible to fit a  function to each peak and perform  integration,\n",
    "we can also analyze our spectrum by doing numerical integration.\n",
    "\n",
    "SciPy offers a number of functions for [doing numerical integration](https://docs.scipy.org/doc/scipy/tutorial/integrate.html) (see \"Methods for Integrating Functions given fixed samples\"). \n",
    "We will somewhat arbitrarily use `simpson`.\n",
    "We imported this at the beginning of the notebook with `from scipy.integrate import simpson`.\n",
    "As a reminder, our peaks look like the following:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c868faf-4f9f-48dc-b287-70af659299a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig_corrected = px.line(si_peaks, x=\"energy\", y=\"corrected_intensity\")    \n",
    "fig_corrected.update_xaxes(autorange=\"reversed\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81ee9392-d1cb-4f49-ba7b-109c491237f8",
   "metadata": {},
   "source": [
    "To use [`simpson` for integration](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.simpson.html#scipy.integrate.simpson), you first pass in your `y` values, then your `x` values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07eb58c7-c191-4d63-99f9-198a7e6afc1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "si_peak_area = simpson(y=si_peaks[\"corrected_intensity\"], x=si_peaks[\"energy\"])\n",
    "print(si_peak_area)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76a9f2ec-e4f0-4832-9a70-aef6efd1f44a",
   "metadata": {},
   "source": [
    "Notice that our area is negative. \n",
    "This is because our energies are listed in descending order (highest to lowest).\n",
    "This matters when we are doing numerical integration.\n",
    "We can fix this by sorting our dataframe based on the energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ac6a4a1-d4d9-4d38-8687-c61c111aacc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "si_peaks.sort_values(by=\"energy\", inplace=True)\n",
    "si_peaks.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb635783-82fb-4d9f-a1f5-dac8fb30edc7",
   "metadata": {},
   "source": [
    "Now we can repeat the integration and get the correct sign."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb1bcdfe-4d5c-49bd-be57-5037f77a416b",
   "metadata": {},
   "outputs": [],
   "source": [
    "si_peak_area = simpson(y=si_peaks[\"corrected_intensity\"], x=si_peaks[\"energy\"])\n",
    "print(si_peak_area)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d83f1cc-81c9-4222-ba42-37ebbe10c166",
   "metadata": {},
   "source": [
    "<div class=\"exercise admonition\">\n",
    "<p class=\"admonition-title\">Exercise</p>\n",
    "\n",
    "Repeat the integration for the oxygen peak.\n",
    "\n",
    "After you have calculated the area of the oxygen peak, figure out the ratio of Si to O in you material by taking the ratio of the calculated areas.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d9007712-db67-4251-b675-dfba9a5dd677",
   "metadata": {},
   "source": [
    "### Fitting Peaks with Gaussian Functions\n",
    "\n",
    "The following is supplemental material. \n",
    "If you would like to be more rigorous with your peak fitting, you might choose to fit the peaks with a Gaussian instead of performing numerical integration.\n",
    "\n",
    "In order to integrate peaks in XPS, peaks are often fit with a Gaussian function.\n",
    "The Gaussian function is represented by the equation below:\n",
    "\n",
    "$$\n",
    "f(x) = A \\exp \\left( -\\frac{(x - B)^2}{2C^2} \\right)\n",
    "$$\n",
    "\n",
    "Where:\n",
    "- $A$ is the amplitude of the peak.\n",
    "- $B$ is the position of the center of the peak (mean or expected value).\n",
    "- $C$ determines the width of the Gaus\n",
    "In order to fit this data, we will use a function from `sclled `curve_fit`.led `c.\n",
    "`curve_fit` allows you to fit parameters for functions that you \n",
    "\n",
    "Sometimes, a peak may be fit with multiple Gaussian if one does not describe the peak well.define.\n",
    "\n",
    "In order to use `curve_fit`, we must first write aon for  functithe equation we would like to fit.\n",
    "The cell below defines a Gaussian function using the parameters $A$, $B$, a\n",
    "We also define a function called `double_gaussian` that contains a fit with two Gaussian functions added together.riance.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c3d6e9e-54d5-48e3-a16a-17c1244e15c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gaussian(x, A, B, C):\n",
    "    y = A*np.exp(-(x-B)**2/(2*C**2)) \n",
    "    return y\n",
    "\n",
    "def double_gaussian(x, a1, b1, c1, a2, b2, c2):\n",
    "    g1 = a1 * np.exp(-(x - b1)**2 / (2 * c1**2))\n",
    "    g2 = a2 * np.exp(-(x - b2)**2 / (2 * c2**2))\n",
    "    return g1 + g2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4c0aaf1-f830-408c-8994-ee00a96b0c0c",
   "metadata": {},
   "source": [
    "To peform our fit, we'll first need to pick a peak of interest.\n",
    "We will do this in a way that is very similar to what we did for the background. \n",
    "This time, we can use `between`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a67d9ce3-a7a6-4526-b7f0-9d4f10c92b7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pick a peak \n",
    "peak_1 = si_peaks[si_peaks[\"energy\"].between(99.25, 101.5)]\n",
    "peak_1.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e866574-2ecd-48ef-9cf9-742354001a5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualize the peak\n",
    "fig_narrow = px.line(peak_1, x=\"energy\", y=\"intensity\")\n",
    "fig_narrow.update_xaxes(autorange=\"reversed\")\n",
    "fig_narrow.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cc2d189-a871-4b5a-b36b-4c944b9fb255",
   "metadata": {},
   "source": [
    "We can see that this peak doesn't look symmetrical - that is, it may not be well-described by a Gaussian. \n",
    "We will try fitting it with our two Gaussian function.\n",
    "\n",
    "In the cell below,  we use the `curve_fit` function. \n",
    "The first argument is the function you want to fit, followed by the x values, then y values.\n",
    "Finally, we specify `p0`, this is an initial guess for the parameters.\n",
    "By examining the peak in the graph above, we might specify an initial guess of 30,000 for the peak height, 100 for the center of the peak, and 1 for the peak width for both Gaussians."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "841cc813-78af-468c-b558-247cf4fdd3aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "parameters_1, covariance_1 = curve_fit(double_gaussian, peak_1[\"energy\"], peak_1[\"intensity\"], \n",
    "                                       p0=[30000, 100, 0.1, 30000, 100, 0.1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b939762a-9dc4-4daf-9bad-1e8fb10237b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(parameters_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae14bf90-8716-490f-b3db-2d12c01efba3",
   "metadata": {},
   "source": [
    "We can compare our model to our actual peak by putting our x values into our gaussian function, along with our peak parameters.\n",
    "In the cell below, we use `*parameters_1`. \n",
    "When we use this with a list, it fills in each element into the function arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1aa8dff2-0cff-42c2-9df6-06464db061b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get fit values\n",
    "fit_y = double_gaussian(si_peaks[\"energy\"], *parameters_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "381da71c-d790-4d92-9fdd-74cbe5b6b944",
   "metadata": {},
   "source": [
    "Taking a wider view, we can examine our fit against the original Si peaks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae36b456-551f-4d1f-a7f1-5ba8243302ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig_fit = px.line(si_peaks, x=\"energy\", y=\"corrected_intensity\")\n",
    "fig_fit.add_scatter(x=si_peaks[\"energy\"], y=fit_y)\n",
    "fig_fit.update_xaxes(autorange=\"reversed\")\n",
    "fig_fit.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17046b5b-797c-4380-b8e3-2528b1863e67",
   "metadata": {},
   "source": [
    "### Peak Integration\n",
    "\n",
    "Now that we have our peak, we can integrate it using SciPy.\n",
    "We will use [the `quad` function](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) we imported at the top of the file.\n",
    "In the quad function, you first pass the function you would like to integrate, followed by the lower and upper integration bounds.\n",
    "In our case, we need to pass in our fit parameters as arguments to the `double_gaussian` function as well. That is done with the \n",
    "syntax `args=(*parameters_1,)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8acac0d0-a46f-4ae4-bd11-14e89ed6105d",
   "metadata": {},
   "outputs": [],
   "source": [
    "result, error = quad(double_gaussian, 98, 102, args=(*parameters_1,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6814bda-d3ed-4ecd-8dbc-56b070c5247b",
   "metadata": {},
   "outputs": [],
   "source": [
    "result"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93000dad-7bce-4f94-960d-bc800f49b345",
   "metadata": {},
   "source": [
    "<div class=\"exercise admonition\">\n",
    "<p class=\"admonition-title\">Exercise</p>\n",
    "\n",
    "Repeat the analysis for the second Si peak and for the O peak.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8574e317-64da-4dd7-a149-44679d3b1bbf",
   "metadata": {},
   "source": [
    "<div class=\"exercise admonition\">\n",
    "<p class=\"admonition-title\">Exercise</p>\n",
    "\n",
    "Repeat the full analysis starting from background subtraction for the oxygen peak.\n",
    "\n",
    "\n",
    "\n",
    "Once you have fit and integrated the oxygen peak, you can compare the ratio of Si to O in the material by\n",
    "taking the ratio of the area of the silicon peaks to the oxygen peak.\n",
    "</div>"
   ]
  }
 ],
 "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.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}