spec_software.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "%matplotlib qt\n",
    "\n",
    "try:\n",
    "    from tkinter import Tk\n",
    "    from tkFileDialog import askopenfilenames\n",
    "except:\n",
    "    from tkinter import Tk\n",
    "    from tkinter import filedialog"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gaussian(x, N, mu, sigma):\n",
    "    return N*np.exp(-0.5*((x-mu)/sigma)**2)\n",
    "\n",
    "def doubleGaussian(x, N1, mu1, sigma1, N2, mu2, sigma2):\n",
    "    return N1*np.exp(-0.5*((x-mu1)/sigma1)**2) + N2*np.exp(-0.5*((x-mu2)/sigma2)**2) \n",
    "\n",
    "def evaluate_fwhm(x, y, threshold = 0.05):\n",
    "    y_norm = y / np.max(y)\n",
    "\n",
    "    try:\n",
    "        i_min = np.where(y_norm >= 0.5*(1.-threshold))[0][0]\n",
    "        i_max = np.where(y_norm >= 0.5*(1.-threshold))[0][-1] \n",
    "        if i_min == i_max:\n",
    "            i_min -= 1\n",
    "            i_max += 1\n",
    "        x_min = x[i_min]\n",
    "        x_max = x[i_max]\n",
    "            \n",
    "        fwhm = x_max - x_min\n",
    "    except:\n",
    "        fwhm = np.inf\n",
    "\n",
    "    return fwhm, np.nan\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set-up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "use_calibrated = False #default = True\n",
    "fit_type = 'dgauss' #Use 'dgauss' to fit a double gaussian (default), 'gauss' to fit a normal gaussian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calibration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "root = Tk()\n",
    "root.withdraw() \n",
    "root.wm_attributes('-topmost', 1) \n",
    "calib_file_name = filedialog.askopenfilenames(parent = root) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "chanel_cal, energy_cal = np.loadtxt(calib_file_name[0], unpack = 'True')\n",
    "\n",
    "#NB: Ordered from n-th order to 0-th order\n",
    "calib_coeff = np.polyfit(chanel_cal, energy_cal, deg = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "mock_channel = np.linspace(0,2000,2000)\n",
    "cal_line = mock_channel*calib_coeff[0] + calib_coeff[1]\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot()\n",
    "\n",
    "ax.set_title('Calibration Plot')\n",
    "ax.plot(mock_channel, cal_line, color = 'grey', label = 'Cal. Line')\n",
    "ax.scatter(chanel_cal, energy_cal, marker = '+', color ='red', label = 'Cal. Points')\n",
    "ax.set_xlabel('Channel')\n",
    "ax.set_ylabel('Energy [keV]')\n",
    "ax.legend()\n",
    "ax.grid()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Data Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "root = Tk()\n",
    "root.withdraw() \n",
    "root.wm_attributes('-topmost', 1) \n",
    "data_file_name = filedialog.askopenfilenames(parent = root) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "chanel, counts = np.loadtxt(data_file_name[0], unpack = True, skiprows=10, delimiter =';', usecols=(0,1))\n",
    "if (use_calibrated):\n",
    "    calib_energies =  chanel*calib_coeff[0] + calib_coeff[1]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "if use_calibrated:\n",
    "    fig, ax = plt.subplots(1,2, figsize = (15,7))\n",
    "    plt.suptitle(data_file_name[0])\n",
    "\n",
    "    ax[0].step(chanel, counts, color = 'gray')\n",
    "    ax[0].set_xlabel('Channel')\n",
    "    ax[0].set_ylabel('Counts')\n",
    "    ax[0].grid()\n",
    "\n",
    "    ax[1].step(calib_energies, counts, color = 'gray')\n",
    "    ax[1].set_xlabel('Energy [keV]')\n",
    "    ax[1].set_ylabel('Counts')\n",
    "    ax[1].grid()\n",
    "else:\n",
    "    fig, ax = plt.subplots(1,1, figsize = (7,7))\n",
    "    plt.suptitle(data_file_name[0])\n",
    "\n",
    "    ax.step(chanel, counts, color = 'gray')\n",
    "    ax.set_xlabel('Channel')\n",
    "    ax.set_ylabel('Counts')\n",
    "    ax.grid()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data File:  C:/Users/lisaf/Desktop/due2lab/test_Am\n",
      "Calibration:  C:/Users/lisaf/Desktop/due2lab/calib.txt\n",
      "\n",
      "\n",
      "Fitting range: 80 - 120 ch\n",
      "Fitting model: double gaussian\n",
      "\n",
      "\n",
      "FIT PARAMETERS:\n",
      "N1:     300.40 ± 27.00 counts\n",
      "Mu1:    94.05 ± 1.41 ch\n",
      "Sigma1: 12.08 ± 0.65 ch\n",
      "N2:     1521.90 ± 37.65 counts\n",
      "Mu2:    103.91 ± 0.07 ch\n",
      "Sigma2: 5.57 ± 0.10 ch\n",
      "\n",
      "\n",
      "FWHM (main peak): 13.11 ± 0.24 ch\n"
     ]
    }
   ],
   "source": [
    "#CAREFUL! If use_calibrated = True, those two values are energies, otherwhise they represent channel numbers\n",
    "min_val = 80\n",
    "max_val = 120\n",
    "\n",
    "if use_calibrated:\n",
    "      mask = np.logical_and(calib_energies >= min_val, calib_energies <= max_val)\n",
    "      energies_to_fit = calib_energies[mask]\n",
    "      x = energies_to_fit\n",
    "else:\n",
    "      mask = np.logical_and(chanel >= min_val, chanel <= max_val)\n",
    "      chanel_to_fit = chanel[mask]\n",
    "      x = chanel_to_fit\n",
    "\n",
    "counts_to_fit = counts[mask]\n",
    "\n",
    "x_fit = np.linspace(np.min(x), np.max(x),1000)\n",
    "\n",
    "if fit_type == 'dgauss': \n",
    "\n",
    "      p0 = [np.max(counts_to_fit), np.min(x), 1, \n",
    "            np.max(counts_to_fit), np.min(x), 1]\n",
    "      bounds = [(0,min(x),0,0,min(x),0),\n",
    "            (np.inf, max(x),np.inf, np.inf, max(x), np.inf)]\n",
    "\n",
    "      \n",
    "      par, cov = curve_fit(doubleGaussian, xdata = x, ydata = counts_to_fit, p0 = p0, bounds = bounds)      \n",
    "      y_fit = doubleGaussian(x_fit, *par)\n",
    "\n",
    "      #fwhm, err_fwhm = evaluate_fwhm(x_fit, y_fit)\n",
    "      fwhm = par[5]*2.355\n",
    "      #Must be corrected to take into account covariance\n",
    "      errs = np.sqrt(np.diag(cov))\n",
    "      err_fwhm = errs[5]*2.355\n",
    "\n",
    "elif fit_type == 'gauss':\n",
    "      p0 = [np.max(x), np.min(x), 1]\n",
    "      bounds = [(0,min(x),0), (np.inf, max(x),np.inf)]\n",
    "\n",
    "      par, cov = curve_fit(gaussian, xdata = x, ydata = counts_to_fit, p0 = p0, bounds = bounds)\n",
    "      y_fit = gaussian(x_fit, *par)\n",
    "\n",
    "      fwhm = par[2]*2.355\n",
    "      #Must be corrected to take into account covariance\n",
    "      errs = np.sqrt(np.diag(cov))\n",
    "      err_fwhm = errs[2]*2.355\n",
    "\n",
    "else:\n",
    "      raise Exception(\"Select a valid fit function: either double gaussian (\\'dgauss\\'), or gaussian (\\'gauss\\')\")\n",
    "\n",
    "fig, ax = plt.subplots(1,1, figsize = (9,7))\n",
    "ax.step(x, counts_to_fit, color = 'black', label = 'Data')\n",
    "if use_calibrated:\n",
    "      ax.set_xlabel('Energy [keV]')\n",
    "      um = 'keV'\n",
    "else:\n",
    "      ax.set_xlabel('Channel')\n",
    "      um = 'ch'\n",
    "ax.set_ylabel('Counts')\n",
    "ax.grid()\n",
    "\n",
    "ax.plot(x_fit, y_fit, color = 'red', label = 'Fit')\n",
    "ax.legend()\n",
    "\n",
    "print('Data File: ', data_file_name[0])\n",
    "print('Calibration: ', calib_file_name[0])\n",
    "print('\\n')\n",
    "print('Fitting range: {0} - {1} {2}'.format(min_val, max_val, um))\n",
    "print('Fitting model: double gaussian')\n",
    "print('\\n')\n",
    "print('FIT PARAMETERS:')\n",
    "print('N1:     {0:.2f} ± {1:.2f} counts'.format(par[0],errs[0]))\n",
    "print('Mu1:    {0:.2f} ± {1:.2f} {2}'.format(par[1],errs[1], um))\n",
    "print('Sigma1: {0:.2f} ± {1:.2f} {2}'.format(par[2],errs[2], um))\n",
    "if fit_type == 'dgauss':\n",
    "      print('N2:     {0:.2f} ± {1:.2f} counts'.format(par[3],errs[3]))\n",
    "      print('Mu2:    {0:.2f} ± {1:.2f} {2}'.format(par[4],errs[4], um))\n",
    "      print('Sigma2: {0:.2f} ± {1:.2f} {2}'.format(par[5],errs[5], um))\n",
    "print('\\n')\n",
    "\n",
    "print('FWHM (main peak): {0:.2f} ± {1:.2f} {2}'.format(fwhm, err_fwhm, um))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "astro",
   "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.12.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}