c6-data-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting code using Gaussian modulation patterns on simulated measurement data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import curve_fit\n",
    "import os\n",
    "import re\n",
    "import os.path\n",
    "from os import path\n",
    "from definitions import *\n",
    "from instrument import *\n",
    "from util import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_P(I_up,I_down, dI_up, dI_down):\n",
    "    P = (I_up - I_down) / (I_up + I_down)\n",
    "    dP = 2 * np.sqrt((I_down * dI_up)**2 + (I_up * dI_down)**2) / (I_up + I_down) ** 2 \n",
    "    return P, dP\n",
    "\n",
    "def compute_I_diff(I_up,I_down, dI_up, dI_down):\n",
    "    I_diff = (I_up - I_down)\n",
    "    dI_diff = np.sqrt(dI_up ** 2 + dI_down ** 2)\n",
    "    return I_diff, dI_diff\n",
    "\n",
    "def compute_rms(x, dx):\n",
    "    rms = np.sqrt(np.mean(np.square(x)))\n",
    "    squared_errors = np.square(dx)\n",
    "    d_rms = np.sqrt(np.mean(squared_errors) / len(x))\n",
    "    return rms, d_rms \n",
    "\n",
    "def get_measurement_point(i, folder='data'):\n",
    "    d1 = np.genfromtxt(f'{folder}/up/{i}/det.dat', delimiter=' ', usecols=(0,1,2,3), unpack=True)\n",
    "    d2 = np.genfromtxt(f'{folder}/down/{i}/det.dat', delimiter=' ', usecols=(0,1,2,3), unpack=True)\n",
    "    d3 = np.genfromtxt(f'{folder}/empty_up/{i}/det.dat', delimiter=' ', usecols=(0,1,2,3), unpack=True)\n",
    "    d4 = np.genfromtxt(f'{folder}/empty_down/{i}/det.dat', delimiter=' ', usecols=(0,1,2,3), unpack=True)\n",
    "    up = DetectorReading(*d1)\n",
    "    down = DetectorReading(*d2)\n",
    "    empty_up = DetectorReading(*d3)\n",
    "    empty_down = DetectorReading(*d4)\n",
    "    parameters = extract_parameters(f'{folder}/up/{i}/det.dat')\n",
    "    mp = MeasurementPoint(up, down, empty_up, empty_down, parameters)\n",
    "    return mp\n",
    "    \n",
    "\n",
    "class DetectorReading:\n",
    "    def __init__(self, y, I, I_err, N):\n",
    "        self.y = y * 1e-2\n",
    "        self.I = I\n",
    "        self.dI = I_err\n",
    "        self.N = N\n",
    "\n",
    "class MeasurementPoint:\n",
    "    def __init__(self, up, down, empty_up, empty_down, parameters):\n",
    "        self.up = up\n",
    "        self.down = down\n",
    "        self.empty_up = empty_up\n",
    "        self.empty_down = empty_down\n",
    "        self.parameters = parameters\n",
    "    \n",
    "    def y(self):\n",
    "        return self.up.y\n",
    "    \n",
    "    # Computes just I_up - I_down\n",
    "    def compute_I_diffs(self):\n",
    "        I_diff_s, dI_diff_s = compute_I_diff(self.up.I, self.down.I, self.up.dI, self.down.dI)\n",
    "        I_diff_b, dI_diff_b = compute_I_diff(self.empty_up.I, self.empty_down.I, self.empty_up.dI, self.empty_down.dI)\n",
    "        return self.y(), I_diff_b, dI_diff_b, I_diff_s, dI_diff_s\n",
    "\n",
    "    def compute_P(self):\n",
    "        P_s, dP_s = compute_P(self.up.I, self.down.I, self.up.dI, self.down.dI)\n",
    "        P_b, dP_b = compute_P(self.empty_up.I, self.empty_down.I, self.empty_up.dI, self.empty_down.dI)\n",
    "        return self.y(), P_b, dP_b, P_s, dP_s\n",
    "    \n",
    "    def compute_rms(self, y_range = None):\n",
    "        y, P_b, dP_b, P_s, dP_s = self.compute_P()\n",
    "        # y*=1e-3\n",
    "        # y, P_b, dP_b, P_s, dP_s = self.compute_I_diffs()\n",
    "        if y_range == None:\n",
    "            rms_b, d_rms_b = compute_rms(P_b, dP_b)\n",
    "            rms_s, d_rms_s = compute_rms(P_s, dP_s)\n",
    "        else:\n",
    "            # range is assumed to be of the form (y_min, y_max)\n",
    "            indices = indices_within_range(y, *y_range)\n",
    "            rms_b, d_rms_b = compute_rms(P_b[indices], dP_b[indices])\n",
    "            rms_s, d_rms_s = compute_rms(P_s[indices], dP_s[indices])\n",
    "        return rms_b, d_rms_b, rms_s, d_rms_s\n",
    "\n",
    "    def rms_ratio(self, y_range = None):\n",
    "        rms_b, d_rms_b, rms_s, d_rms_s = self.compute_rms(y_range=y_range)\n",
    "        R = rms_s / rms_b\n",
    "        dR = np.sqrt((d_rms_s / rms_b)**2 + (rms_s * d_rms_b / rms_b**2)**2)\n",
    "        return R, dR\n",
    "    \n",
    "    def plot_intensities(self):\n",
    "        plt.plot(self.y() * 1e3, self.empty_up.I, label='Up')\n",
    "        plt.plot(self.y() * 1e3, self.empty_down.I, label = 'Down')\n",
    "        plt.title(r'Base $I_{up}, I_{down}$')\n",
    "        plt.grid()\n",
    "        plt.legend()\n",
    "        plt.show()    \n",
    "        plt.plot(self.y() * 1e3, self.up.I, label='Up')\n",
    "        plt.plot(self.y() * 1e3, self.down.I, label = 'Down')\n",
    "        plt.title(r'Sample $I_{up}, I_{down}$')\n",
    "        plt.grid()\n",
    "        plt.legend()\n",
    "        plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_gauss(y, I, dI, y_range, By, lambda_0, theta_0, wl_sigma):   \n",
    "    indices = indices_within_range(y, *y_range)\n",
    "    def I_fit(y, A_0, I_0):\n",
    "        return I_empty_analytical(y, lambda_0, By, theta_0, wl_sigma=wl_sigma) * A_0 + I_0\n",
    "    # Todo: use dI for better estimate?\n",
    "    popt, pcov = curve_fit(I_fit, y[indices], I[indices], absolute_sigma=True,sigma = dI[indices])\n",
    "    # print(popt, pcov)\n",
    "    A_0 = popt[0]\n",
    "    dA_0 = np.sqrt(pcov[0][0])\n",
    "    I_0 = popt[1]\n",
    "    dI_0 = np.sqrt(pcov[1][1])\n",
    "    # print(I_0, dI_0)\n",
    "    return A_0, dA_0, I_0, dI_0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gaussian envelope and RMS fitting "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from hankel import *\n",
    "# from definitions import form_factor\n",
    "radii = [500,3000, 20000] # Å\n",
    "\n",
    "def P_Q_corrected(Q_max, delta_min, delta_max, tau, R):\n",
    "\n",
    "    Z = 250\n",
    "    delta, Qn, dQ, N =setup_Q_delta_arrays(delta_min, delta_max, Z)\n",
    "    H_0, H_kernel = create_hankel_kernels(Qn, dQ, delta, Z)\n",
    "    I_Q = form_factor(Qn,R) * 0.5 * R ** 2\n",
    "    G_full, G_full_0 = compute_G_matrices(I_Q, H_0, H_kernel, Z)\n",
    "    G_delta_num, sigma_t_num = compute_total_G_sigma(G_full, G_full_0)\n",
    "    G_partials, sigma_partials, G_0_partial, P_partial = compute_partial_G_sigma(G_full, G_full_0, sigma_t_num, Z)\n",
    "    def get_Q_ix(Q):\n",
    "        return int(Q/dQ)\n",
    "    Q_ix = min([N-1,get_Q_ix(Q_max)])\n",
    "    print(Q_ix)\n",
    "    P_corrected = np.exp(sigma_partials[Q_ix] * tau * (G_0_partial[Q_ix,:] - 1))\n",
    "    return delta, P_corrected\n",
    "\n",
    "def analyse_instr_data(instr, folder = '', N = 30, dynamic_fit_scaling=False, logplot = True, show_corr=False):\n",
    "    for fit_method in ['gauss', 'rms']:\n",
    "        plt.figure(figsize=(10, 6), dpi=300)\n",
    "        print(f\"FITTING USING {fit_method.upper()}\")\n",
    "        for radius in radii:\n",
    "        # for radius in [3000]:\n",
    "            B = np.zeros(N)\n",
    "            delta = np.zeros(N)\n",
    "            P = np.zeros(N)\n",
    "            dP = np.zeros(N)\n",
    "            for i in range(N):\n",
    "                measurement = get_measurement_point(i, f'{folder}data_{instr.id}_{radius}')\n",
    "                # measurement.plot_intensities()\n",
    "                # y, P_b, dP_b, P_s, dP_s = measurement.compute_I_diffs()\n",
    "                # plt.plot(y,P_b)\n",
    "                # plt.show()\n",
    "                y, P_b, dP_b, P_s, dP_s = measurement.compute_P()\n",
    "                # plt.plot(y,P_b)\n",
    "                # plt.show()\n",
    "                By = float(measurement.parameters['By'])\n",
    "                lambda_0 = float(measurement.parameters['L0']) * 1e-10\n",
    "                lambda_sigma = float(measurement.parameters['DL']) * 1e-10\n",
    "                theta_0 = float(measurement.parameters['theta0'])\n",
    "                try: \n",
    "                    instr.L_s = float(measurement.parameters['L_s'])\n",
    "                except:\n",
    "                    pass\n",
    "                B[i] = By * abs(1 - instr.L_1 / instr.L_2)\n",
    "                p_0 = compute_p_0(B[i], instr.L0, instr.theta_0)\n",
    "                delta[i] = compute_z(B[i],instr.theta_0,instr.L0,instr.L_s)\n",
    "                if dynamic_fit_scaling:\n",
    "                    y_lim = np.min([np.max([p_0/2,3e-3]),instr.detector_size / 2 * 0.95]) # m\n",
    "                    print(y_lim, p_0/2, instr.detector_size / 2 * 0.95)\n",
    "                    y_range=(-y_lim,y_lim) # m\n",
    "                else:\n",
    "                    y_lim = 3e-3 # m\n",
    "                    y_range=(-y_lim,y_lim) # m\n",
    "                if fit_method == 'gauss':\n",
    "                    # Gaussian fit\n",
    "                    A_0_b, dA_0_b, I_0_b, dI_0_b = fit_gauss(y, P_b, dP_b, y_range, B[i], lambda_0, theta_0, lambda_sigma)\n",
    "                    A_0_s, dA_0_s, I_0_s, dI_0_s = fit_gauss(y, P_s, dP_s, y_range, B[i], lambda_0, theta_0, lambda_sigma)\n",
    "                    R = A_0_s / A_0_b\n",
    "                    dR = np.sqrt((dA_0_s / A_0_b)**2 + (A_0_s * dA_0_b / A_0_b**2)**2)\n",
    "                    P[i] = R\n",
    "                    dP[i] = dR\n",
    "                elif fit_method == 'rms':\n",
    "                    P[i], dP[i] = measurement.rms_ratio(y_range)\n",
    "            t = float(measurement.parameters['t'])\n",
    "            R = float(measurement.parameters['R']) * 1e-10\n",
    "            \n",
    "            wavelength = float(measurement.parameters['L0']) * 1e-10\n",
    "            delta_range = np.linspace(0, R * 3, 1000)\n",
    "            P_an = P_analytical(delta_range, R, t, wavelength)\n",
    "            plt.errorbar(delta * 1e9,P, ms=3,fmt='o', linestyle=' ', yerr=dP,label=rf'$R = {round(R*1e9)}$ nm, fit')\n",
    "            plt.plot(delta_range * 1e9, P_an, lw = 1.5, label=rf'$R = {round(R*1e9)}$ nm, analytical')\n",
    "            if show_corr:\n",
    "                tau = s_t(R,t,wavelength)\n",
    "                delta_cor, P_cor = P_Q_corrected(instr.Q_max(), R / 10, R * 5, tau, R)\n",
    "                plt.plot(delta_cor * 1e9, P_cor, lw = 1.5, linestyle='--', label=rf'$R = {round(R*1e9)}$ nm, corr')\n",
    "        name = instr.name\n",
    "        if logplot:\n",
    "            plt.xscale('log')\n",
    "        plt.xlim((10, R * 3 * 1e9))\n",
    "        # plt.ylim((0, 1))\n",
    "        \n",
    "        plt.xlabel(r'$\\delta$ [nm]')\n",
    "        plt.ylabel(r'$P(\\delta)$')\n",
    "        if logplot:\n",
    "            plt.legend(loc=3)\n",
    "        else:\n",
    "            plt.legend(loc=1)\n",
    "        \n",
    "        plt.grid()\n",
    "        plt.savefig(f\"docs/simulation-plot-{fit_method}-{name.replace(' ','-')}.eps\",bbox_inches=\"tight\", pad_inches=0, format='eps')\n",
    "        # plt.savefig(f\"docs/simulation-plot-{fit_method}-{name.replace(' ','-')}.eps\",bbox_inches=\"tight\", format='eps')\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Normal simulations "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "intrs = load_instruments('instruments.csv')\n",
    "for instr in intrs[3:-6]:\n",
    "    analyse_instr_data(instr, folder='working-data/', logplot=True, show_corr=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimized instrument"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "with open('optimized_instruments_2.pkl', 'rb') as file: \n",
    "    # Call load method to deserialze \n",
    "    instrs = pickle.load(file) \n",
    "\n",
    "    for (instr,_) in instrs:\n",
    "        prec_type = instr.prec_type\n",
    "        if instr.prec_type == 'wsp':\n",
    "            prec_type = 'iwsp'\n",
    "\n",
    "        if instr.prec_type != 'foil2':\n",
    "            continue\n",
    "        instr.id = instr.name\n",
    "        instr.id=instr.id.replace(' ','_')\n",
    "\n",
    "        instr.prec_type='foil'\n",
    "        analyse_instr_data(instr, '',2, logplot=False)\n",
    "        #analyse_instr_data(instr, '',50, logplot=False)\n",
    "        break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Raw data plots and demonstration\n",
    "It is important to also show some of the data you perform a fitting routine on etc. to illustrate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_lim = 3e-3 # m\n",
    "y_range=(-y_lim,y_lim) # m\n",
    "\n",
    "intrs = load_instruments('instruments.csv')\n",
    "instr = intrs[7]\n",
    "print(instr.name)\n",
    "\n",
    "radius = 3000\n",
    "i = 29\n",
    "measurement = get_measurement_point(i, f'working-data/data_{instr.id}_{radius}')\n",
    "y, P_b, dP_b, P_s, dP_s = measurement.compute_P()\n",
    "By = float(measurement.parameters['By'])\n",
    "lambda_0 = float(measurement.parameters['L0']) * 1e-10\n",
    "lambda_sigma = float(measurement.parameters['DL']) * 1e-10\n",
    "theta_0 = float(measurement.parameters['theta0'])\n",
    "try: \n",
    "    instr.L_s = float(measurement.parameters['L_s'])\n",
    "except:\n",
    "    pass\n",
    "delta = compute_z(By,instr.theta_0,instr.L0,instr.L_s)\n",
    "print(delta * 1e9, By)\n",
    "t = float(measurement.parameters['t'])\n",
    "R = float(measurement.parameters['R']) * 1e-10\n",
    "\n",
    "# Gauss\n",
    "A_0_b, dA_0_b, I_0_b, dI_0_b = fit_gauss(y, P_b, dP_b, y_range, By, lambda_0, theta_0, lambda_sigma)\n",
    "A_0_s, dA_0_s, I_0_s, dI_0_s = fit_gauss(y, P_s, dP_s, y_range, By, lambda_0, theta_0, lambda_sigma)\n",
    "R = A_0_s / A_0_b\n",
    "dR = np.sqrt((dA_0_s / A_0_b)**2 + (A_0_s * dA_0_b / A_0_b**2)**2)\n",
    "P_gauss = R\n",
    "dP_gauss = dR\n",
    "# RMS\n",
    "rms_b, d_rms_b, rms_s, d_rms_s = measurement.compute_rms(y_range=y_range)\n",
    "P_rms, dP_rms = measurement.rms_ratio(y_range)\n",
    "wavelength = float(measurement.parameters['L0']) * 1e-10\n",
    "\n",
    "P_an_val = P_analytical(delta, radius * 1e-10, t, wavelength)\n",
    "print(\"GAUSS\")\n",
    "print(A_0_b, dA_0_b)\n",
    "print(A_0_s, dA_0_s)\n",
    "print(\"RMS\")\n",
    "print(rms_b, d_rms_b)\n",
    "print(rms_s, d_rms_s)\n",
    "print(P_rms, dP_rms, P_an_val)\n",
    "print(P_gauss, dP_gauss, P_an_val)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I_empty_up = measurement.empty_up.I\n",
    "I_empty_down = measurement.empty_down.I\n",
    "I_up = measurement.up.I\n",
    "I_down = measurement.down.I\n",
    "dI_empty_up = measurement.empty_up.I\n",
    "dI_empty_down = measurement.empty_down.I\n",
    "dI_up = measurement.up.I\n",
    "I_down = measurement.down.I\n",
    "y = measurement.y()\n",
    "\n",
    "def plot_intensity(I, label=''):\n",
    "    plt.plot(y * 1e3, I ,'.',ms=3)\n",
    "    plt.xlabel(r'$y$ [mm]')\n",
    "    plt.ylabel(r'$I$ [a.u.]')\n",
    "    plt.grid()\n",
    "    plt.savefig(f\"docs/simulation-raw-intensity-{label}.eps\",bbox_inches=\"tight\", pad_inches=0, format='eps')\n",
    "    plt.show() \n",
    "\n",
    "\n",
    "for (label,I) in [('b-up',I_empty_up), ('b-down', I_empty_down), ('s-up', I_up), ('s-down',I_down)]:\n",
    "    plot_intensity(I,label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(y* 1e3, I_empty_up ,'.',ms=3, label=r'$I_{b,up}$')\n",
    "plt.plot(y* 1e3, I_up ,'.',ms=3, label=r'$I_{s,up}$')\n",
    "plt.xlabel(r'$y$ [mm]')\n",
    "plt.ylabel(r'$I$ [a.u.]')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.savefig(f\"docs/simulation-raw-intensity-up.eps\",bbox_inches=\"tight\", pad_inches=0, format='eps')\n",
    "plt.show() \n",
    "plt.plot(y* 1e3, I_empty_down ,'.',ms=3, label=r'$I_{b,down}$')\n",
    "plt.plot(y* 1e3, I_down ,'.',ms=3, label=r'$I_{s,down}$')\n",
    "plt.xlabel(r'$y$ [mm]')\n",
    "plt.ylabel(r'$I$ [a.u.]')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.savefig(f\"docs/simulation-raw-intensity-down.eps\",bbox_inches=\"tight\", pad_inches=0, format='eps')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def I_fit(y, A_0, I_0):\n",
    "    return I_empty_analytical(y, lambda_0, By, theta_0, wl_sigma=lambda_sigma) * A_0 + I_0\n",
    "y, P_b, dP_b, P_s, dP_s = measurement.compute_P()\n",
    "\n",
    "plt.plot(y * 1e3, I_empty_up + I_empty_down,'.',ms=3, label=r'$I_{b,up} + I_{b,down}$')\n",
    "plt.plot(y * 1e3, I_up + I_down ,'.',ms=3, label=r'$I_{s,up} + I_{s,down}$')\n",
    "plt.plot(y * 1e3, I_empty_up - I_empty_down,'.',ms=3, label=r'$I_{b,up} - I_{b,down}$')\n",
    "plt.plot(y * 1e3, I_up - I_down ,'.',ms=3, label=r'$I_{s,up} - I_{s,down}$')\n",
    "plt.xlabel(r'$y$ [mm]')\n",
    "plt.ylabel(r'$I$ [a.u.]')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.savefig(f\"docs/simulation-raw-intensity-differential.eps\",bbox_inches=\"tight\", pad_inches=0, format='eps')\n",
    "plt.show() \n",
    "plt.plot(y * 1e3, P_b ,'.',ms=3, label=r'$V_b$')\n",
    "plt.plot(y * 1e3, P_s ,'.',ms=3, label=r'$V_s$')\n",
    "\n",
    "fit_line_s = I_fit(y, A_0_s, I_0_s)\n",
    "fit_line_b = I_fit(y, A_0_b, I_0_b)\n",
    "plt.plot(y * 1e3, fit_line_b, lw = 1.5, label=r'$V_{fit,b}$')\n",
    "plt.plot(y * 1e3, fit_line_s, lw = 1.5, label=r'$V_{fit,s}$')\n",
    "plt.axvline(y_lim* 1e3, linestyle='--', color='red')\n",
    "plt.axvline(-y_lim* 1e3, linestyle='--', color='red')\n",
    "plt.xlabel(r'$y$ [mm]')\n",
    "plt.ylabel(r'$V$ ')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.savefig(f\"docs/simulation-raw-intensity-pol.eps\",bbox_inches=\"tight\", pad_inches=0, format='eps')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bep",
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}