{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Window function and amplitude correction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from radiocalibrationtoolkit import *\n",
    "from scipy.signal import butter, lfilter\n",
    "from scipy.signal import blackmanharris, boxcar, hann\n",
    "from scipy.fft import rfft"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some global plot settings\n",
    "plt.rcParams[\"axes.labelweight\"] = \"bold\"\n",
    "plt.rcParams[\"font.weight\"] = \"bold\"\n",
    "plt.rcParams[\"font.size\"] = 16\n",
    "plt.rcParams[\"legend.fontsize\"] = 12\n",
    "\n",
    "plt.rcParams[\"xtick.major.width\"] = 2\n",
    "plt.rcParams[\"ytick.major.width\"] = 2\n",
    "\n",
    "plt.rcParams[\"xtick.major.size\"] = 5\n",
    "plt.rcParams[\"ytick.major.size\"] = 5\n",
    "\n",
    "plt.rcParams[\"xtick.labelsize\"] = 14\n",
    "plt.rcParams[\"ytick.labelsize\"] = 14"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Examine effects on a single trace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 2048   # samples\n",
    "fs_Hz = 250e6 # sampling frequency\n",
    "signal_freq_Hz = 55.5e6 # signal\n",
    "win = blackmanharris(N) # window function type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tx = np.arange(N) / fs_Hz\n",
    "def make_simple_time_trace(signal_amplitude=0.7, signal_frequency=55.5e+6, noise_amplitude=1, N=2048, apply_filter=True):\n",
    "    time_trace = signal_amplitude * np.sin(2*np.pi*signal_frequency*tx) + noise_amplitude * np.random.normal(size=N)\n",
    "    b, a = butter(5, [30e6, 80e6], fs=250e6, btype='band')\n",
    "    if apply_filter:\n",
    "        time_trace = lfilter(b, a, time_trace)\n",
    "    return np.round(time_trace * N) / N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tt = make_simple_time_trace(signal_amplitude=1)\n",
    "fx = np.fft.rfftfreq(N, 1/fs_Hz)*1e-6\n",
    "\n",
    "fig, ax = plt.subplots(2,1, figsize=(6,6))\n",
    "labels = ['no windows', 'hann', 'blackman-harris']\n",
    "for i, w in enumerate([1, hann(N), blackmanharris(N)]):\n",
    "    ax[0].plot(tx*1e+6, tt*w)\n",
    "    ax[1].plot(fx, 10*np.log10(np.abs(np.fft.rfft(tt*w))), label=labels[i])\n",
    "\n",
    "fig.subplots_adjust(hspace=0.5)\n",
    "ax[1].legend()\n",
    "\n",
    "ax[0].set_xlabel(\"$\\mu$s\")\n",
    "ax[1].set_xlabel(\"frequency [MHz]\")\n",
    "ax[0].set_ylabel(\"amplitude [ADC]\")\n",
    "ax[1].set_ylabel(\"amplitude [ADC]\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate time trace\n",
    "time_trace = make_simple_time_trace()\n",
    "\n",
    "# calculate two sided spectrum\n",
    "spectrum = np.abs(fft(time_trace))\n",
    "# calculate one sided spectrum (not corrected for the one side)\n",
    "rspectrum = np.abs(np.fft.rfft(time_trace))\n",
    "\n",
    "# calculate one sided spectrum (corrected for the one side)\n",
    "r2spectrum = correct_energy_of_one_sided_spectrum(rspectrum)\n",
    "\n",
    "# calculate amplitude window function correction\n",
    "Aw = N/np.sum(win) /np.sqrt(2)\n",
    "\n",
    "# calculate two sided spectrum using a window function (not corrected for the one side)\n",
    "spectrum_w = np.abs(fft(time_trace*win))\n",
    "# calculate one sided spectrum using a window function (not corrected for the one side)\n",
    "rspectrum_w = np.abs(np.fft.rfft(time_trace*win))\n",
    "\n",
    "# calculate one sided spectrum (corrected for the one side)\n",
    "r2spectrum = correct_energy_of_one_sided_spectrum(rspectrum)\n",
    "r2spectrum_w = correct_energy_of_one_sided_spectrum(rspectrum_w)\n",
    "\n",
    "# define X-axis values\n",
    "fx = np.fft.rfftfreq(N, 1/fs_Hz)/1e+6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate energy from time trace, one and two sided spectrum with and without window\n",
    "# when the window is used, the amplitudes are corrected\n",
    "print(\n",
    "    f\"Energy calculated from:\\n\"\n",
    "    f\"time trace = {get_energy_from_time_trace(time_trace):.2f},\\n\"\n",
    "    f\"two sided spectrum = {get_energy_from_two_sided_spectrum(spectrum):.2f},\\n\"\n",
    "    f\"one sided spectrum = {get_energy_from_one_sided_spectrum(rspectrum):.2f},\\n\"\n",
    "    f\"windowed two sided spectrum = {get_energy_from_two_sided_spectrum(spectrum_w):.2f},\\n\"\n",
    "    f\"one sided spectrum corrected for being the one sided spectrum = {get_energy_from_one_sided_spectrum_corrected4one_side(r2spectrum):.2f},\\n\"\n",
    "    f\"windowed two sided spectrum corrected by window function loss = {get_energy_from_two_sided_spectrum(spectrum_w * Aw):.2f},\\n\"\n",
    "    f\"windowed one sided spectrum corrected by window function loss = {get_energy_from_one_sided_spectrum(rspectrum_w * Aw):.2f},\\n\"\n",
    "    f\"one sided spectrum corrected for being the one sided spectrum and for window function loss = {get_energy_from_one_sided_spectrum_corrected4one_side(r2spectrum_w*Aw):.2f},\\n\"\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trace with a broad band pulse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# create a broad band noise trace\n",
    "bb_spec=np.zeros(65)\n",
    "bb_spec[26:33] = 200\n",
    "\n",
    "std_dev = 10  # Standard deviation\n",
    "\n",
    "# Randomize the array values using a normal distribution\n",
    "bb_spec = bb_spec + std_dev * np.random.randn(len(bb_spec))\n",
    "bb_spec = np.abs(bb_spec)\n",
    "\n",
    "bb_tt = np.real(ifft(one_sided_2_complex_two_sided(bb_spec)))\n",
    "\n",
    "fx_bb = np.arange(0, bb_spec.size)/(2*(bb_spec.size-1)/250)\n",
    "tx_bb = np.arange(bb_tt.size)/250\n",
    "\n",
    "fig, ax = plt.subplots(2,1, figsize=(6,6))\n",
    "\n",
    "ax[0].set_xlabel(\"frequency [MHz]\")\n",
    "ax[0].set_ylabel(\"amplitude [ADC]\")\n",
    "ax[1].set_ylabel(\"amplitude [ADC]\")\n",
    "ax[1].set_xlabel(\"$\\mu$s\")\n",
    "\n",
    "ax[0].plot(fx_bb, bb_spec)\n",
    "ax[1].plot(tx_bb, bb_tt)\n",
    "\n",
    "fig.subplots_adjust(hspace=0.35)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(np.abs(rfft(bb_tt)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add broad band pulse to the time trace\n",
    "time_trace_e = time_trace.copy()\n",
    "time_trace_e[np.arange(1024, 1024+bb_tt.size)] = time_trace_e[np.arange(1024, 1024+bb_tt.size)] + bb_tt\n",
    "\n",
    "rspectrum_e = np.abs(np.fft.rfft(time_trace_e))\n",
    "rspectrum_w_e = np.abs(np.fft.rfft(time_trace_e*win))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_plots(time_trace, show_window_func=True):\n",
    "    rspectrum = np.abs(np.fft.rfft(time_trace))\n",
    "    rspectrum_w = np.abs(np.fft.rfft(time_trace * win))\n",
    "    Aw = N / np.sum(win) / np.sqrt(2)\n",
    "\n",
    "    e_tt = get_energy_from_time_trace(time_trace)\n",
    "    e_spec = get_energy_from_one_sided_spectrum(rspectrum)\n",
    "    e_spec_w = get_energy_from_one_sided_spectrum(rspectrum_w * Aw)\n",
    "\n",
    "    fig, ax = plt.subplots(2, 1, figsize=(6, 6))\n",
    "\n",
    "    mu = 1e6\n",
    "    ax[0].plot(\n",
    "        tx * mu, time_trace, label=\"no window: E={:.1f} a.u.\".format(e_tt), alpha=0.7\n",
    "    )\n",
    "    ax[0].plot(tx * mu, time_trace * win, label=\"windowed\", alpha=0.7, c='r')\n",
    "    \n",
    "    if show_window_func:\n",
    "        ax[0].plot(tx * mu, win * 10, label=\"window func. x10\", alpha=1, c='g', lw=2)\n",
    "\n",
    "    ax[1].plot(\n",
    "        fx , rspectrum, label=\"no window:         E={:.1f} a.u.\".format(e_spec), alpha=0.7\n",
    "    )\n",
    "\n",
    "    ax[1].plot(\n",
    "        fx ,\n",
    "        rspectrum_w * Aw,\n",
    "        label=\"windowed & A$_w$: E={:.1f} a.u.\".format(e_spec_w),\n",
    "        alpha=0.7,\n",
    "        c='r'\n",
    "    )\n",
    "\n",
    "    display(\n",
    "        get_energy_from_time_trace(time_trace),\n",
    "        get_energy_from_one_sided_spectrum(rspectrum),\n",
    "        get_energy_from_one_sided_spectrum(rspectrum_w * Aw),\n",
    "    )\n",
    "\n",
    "    ax[0].set_xlabel(\"$\\mu$s\")\n",
    "    ax[1].set_xlabel(\"frequency [MHz]\")\n",
    "    ax[0].set_ylabel(\"amplitude [ADC]\")\n",
    "    ax[1].set_ylabel(\"amplitude [ADC]\")\n",
    "    ax[0].legend()\n",
    "    ax[1].legend()\n",
    "    ylims = ax[0].get_ylim()\n",
    "    \n",
    "    ax[0].set_ylim(-abs(np.max(ylims)*2), abs(np.max(ylims))*1.5)\n",
    "    ylims = ax[1].get_ylim()\n",
    "    ax[1].set_ylim(ylims[0], ylims[1]*1.4)\n",
    "    ax[1].set_xlim(20, 90)\n",
    "    \n",
    "    fig.subplots_adjust(hspace=0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# trace without the broad band pulse\n",
    "show_plots(time_trace, show_window_func=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# trace with the broad band pulse\n",
    "time_trace_e = time_trace.copy()\n",
    "time_trace_e[np.arange(1024, 1024+bb_tt.size)] = time_trace_e[np.arange(1024, 1024+bb_tt.size)] + bb_tt\n",
    "show_plots(time_trace_e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rolled trace to shift the broad band pulse\n",
    "show_plots(np.roll(time_trace_e, 500))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The learning here is that a amplitudes of a trace with broad band pulse cannot be safely corrected\n",
    "after the window function because the window function is symmetric and the broad band pulse appearing \n",
    "in random parts of the trace will be each time differently supressed by the window function.\n",
    "\n",
    "Note that the energy of the trace is after the rolling of the trace still the same when no window is used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Distributions of the energy ratios of not windowed spectra and windowed spectra with amplitude correction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_averaged_spectra_and_diffs(arr):\n",
    "    diffs = np.array([])\n",
    "    n, N = arr.shape\n",
    "    avr_rspectrum = 0\n",
    "    avr_rspectrum_w = 0\n",
    "    win = blackmanharris(N)\n",
    "    Aw = N / np.sum(win) / np.sqrt(2)\n",
    "    for i in range(n):\n",
    "        time_trace = arr[i, :]\n",
    "        rspectrum = np.abs(np.fft.rfft(time_trace))\n",
    "        rspectrum_w = np.abs(np.fft.rfft(time_trace * win))\n",
    "        diffs = np.append(\n",
    "            diffs,\n",
    "            get_energy_from_time_trace(time_trace)\n",
    "            / get_energy_from_one_sided_spectrum(rspectrum_w * Aw),\n",
    "        )\n",
    "        avr_rspectrum += rspectrum\n",
    "        avr_rspectrum_w += rspectrum_w\n",
    "\n",
    "    avr_rspectrum /= n\n",
    "    avr_rspectrum_w /= n\n",
    "    return avr_rspectrum, avr_rspectrum_w, diffs\n",
    "\n",
    "\n",
    "def show_results(\n",
    "    avr_rspectrum,\n",
    "    avr_rspectrum_w,\n",
    "    diffs,\n",
    "    histo_edge=2,\n",
    "    N=2048,\n",
    "    bins=None,\n",
    "    xax_min=None,\n",
    "    xax_max=None,\n",
    "    xlim= [None, None]\n",
    "):\n",
    "    fig, ax = plt.subplots()\n",
    "    fx = np.fft.rfftfreq(N, 1 / fs_Hz) / 1e6\n",
    "\n",
    "    if bins == None:\n",
    "        bins = linspace_with_middle_value(\n",
    "            np.mean(diffs), histo_edge * np.std(diffs), 20\n",
    "        )\n",
    "    ax.hist(diffs, bins=bins)\n",
    "    # Calculate mean and standard deviation\n",
    "    mean_diff = np.mean(diffs)\n",
    "    std_diff = np.std(diffs)\n",
    "\n",
    "    print(mean_diff)\n",
    "    print(std_diff)\n",
    "\n",
    "    # Add text box with mean and standard deviation\n",
    "    text_box = f\"$\\mu$: {mean_diff:.2f}\\n$\\sigma$: {std_diff:.2f}\"\n",
    "    ax.text(\n",
    "        0.95,\n",
    "        0.95,\n",
    "        text_box,\n",
    "        transform=ax.transAxes,\n",
    "        verticalalignment=\"top\",\n",
    "        horizontalalignment=\"right\",\n",
    "        bbox=dict(boxstyle=\"round\", facecolor=\"white\", alpha=0.5),\n",
    "    )\n",
    "\n",
    "    ax.set_xlim(xax_min, xax_max)\n",
    "    ax.set_xlabel(\n",
    "        r\"$\\frac{\\mathrm{energy: \\ spectrum \\ with \\ no\\ window}}{\\mathrm{energy:spectrum \\ with\\ window,\\ amplitudes\\ corrected}}$\"\n",
    "    )\n",
    "    ax.set_ylabel(\"entries\")\n",
    "\n",
    "    fig, ax = plt.subplots()\n",
    "\n",
    "    ax.plot(fx, voltageAmp2dB(avr_rspectrum), alpha=0.7, label=\"no window\")\n",
    "    ax.plot(fx, voltageAmp2dB(avr_rspectrum_w), alpha=0.7, label=\"with window\", c='g')\n",
    "    ax.plot(\n",
    "        fx,\n",
    "        voltageAmp2dB(avr_rspectrum_w * Aw),\n",
    "        alpha=0.7,\n",
    "        label=\"with window \" + \"\\n\" + \"and corrected\", c='r'\n",
    "    )\n",
    "    \n",
    "    \n",
    "    # ax.axes.axvspan(30,80, color='b', alpha=0.1)\n",
    "    ax.set_xlim(*xlim)\n",
    "    autoscale_y(ax)     \n",
    "    \n",
    "    ax.set_xlabel(\"frequency [MHz]\")\n",
    "    ax.set_ylabel(\"dB\")\n",
    "    ax.legend()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simple traces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a set of 1000 time traces\n",
    "n = 1000\n",
    "avr_rspectrum = 0\n",
    "avr_rspectrum_w = 0\n",
    "time_traces = make_simple_time_trace()\n",
    "for i in range(n-1):\n",
    "    time_traces = np.vstack((time_traces, make_simple_time_trace()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate energy ratios of unwindowed spectra and windowed with amplitude correction\n",
    "# and average spectra\n",
    "avr_rspectrum, avr_rspectrum_w, diffs = get_averaged_spectra_and_diffs(time_traces)\n",
    "diffs1 = diffs\n",
    "label4final_histo = ['simple traces']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_results(avr_rspectrum, avr_rspectrum_w, diffs, histo_edge=3, xlim=[45, 65])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a set of 1000 time traces with bb pulses\n",
    "n = 1000\n",
    "avr_rspectrum = 0\n",
    "avr_rspectrum_w = 0\n",
    "time_traces = make_simple_time_trace()\n",
    "for i in range(n-1):\n",
    "    time_trace = make_simple_time_trace()\n",
    "    time_trace[0:bb_tt.size] = time_trace[0:bb_tt.size] + bb_tt\n",
    "    time_trace = np.roll(time_trace_e, random.randint(0, 2048 - bb_tt.size) )\n",
    "    time_traces = np.vstack((time_traces, time_trace))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "avr_rspectrum, avr_rspectrum_w, diffs = get_averaged_spectra_and_diffs(time_traces)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_results(avr_rspectrum, avr_rspectrum_w, diffs, histo_edge=3, bins='auto')\n",
    "# ,xax_min=-10, xax_max=80)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mock traces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read HW response\n",
    "hw_file_path = \"./antenna_setup_files/HardwareProfileList_realistic.xml\"\n",
    "hw_dict = read_hw_file(hw_file_path, interp_args={\"fill_value\": \"extrapolate\"})\n",
    "\n",
    "hw_reponse_1 = hw_dict[\"RResponse\"][\"LNA\"]\n",
    "hw_reponse_2 = hw_dict[\"RResponse\"][\"digitizer\"]\n",
    "hw_reponse_3 = hw_dict[\"RResponse\"][\"cable_fromLNA2digitizer\"]\n",
    "hw_reponse_4 = hw_dict[\"RResponse\"][\"impedance_matching_EW\"]\n",
    "\n",
    "\n",
    "# merge all hw responses to one function\n",
    "def hw_response_func(x):\n",
    "    return dB2PowerAmp(\n",
    "        hw_reponse_1(x) + hw_reponse_2(x) + hw_reponse_3(x) + hw_reponse_4(x)\n",
    "    )\n",
    "\n",
    "\n",
    "# impedance function\n",
    "impedance_func = hw_dict[\"IImpedance\"][\n",
    "    \"antenna_EW\"\n",
    "]\n",
    "\n",
    "# read sidereal voltage square spectral density\n",
    "sidereal_voltage2_density_DF = pd.read_csv(\n",
    "    \"./voltage2_density/voltage2_density_Salla_EW_GSM16.csv\",\n",
    "    index_col=0,\n",
    ")\n",
    "sidereal_voltage2_density_DF.columns = sidereal_voltage2_density_DF.columns.astype(\n",
    "    float\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mock_trace_generator = Mock_trace_generator(\n",
    "    sidereal_voltage2_density_DF=sidereal_voltage2_density_DF,\n",
    "    hw_response_func=hw_response_func,\n",
    "    impedance_func=impedance_func,\n",
    "    voltage2ADC=2048,\n",
    "    time_trace_size=2048,\n",
    "    sampling_frequency_MHz=250,\n",
    ")\n",
    "freq_MHz_bins = mock_trace_generator.get_frequency_bins()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "piko = 1e-12\n",
    "additional_noise = 5e-4*piko\n",
    "debug_spectra_dict = mock_trace_generator.generate_mock_trace(\n",
    "    1,\n",
    "    lst=15,\n",
    "    temp_celsius=30,\n",
    "    nbi={\"67.2\": 1},\n",
    "    nbi_err=0.2,\n",
    "    return_debug_dict=True,\n",
    "    additional_noise=additional_noise,\n",
    ")[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "number_of_traces = 1000\n",
    "mock_traces_DF = mock_trace_generator.generate_mock_trace(\n",
    "    number_of_traces,\n",
    "    temp_celsius=[-10,30],\n",
    "    additional_noise=additional_noise,\n",
    "    nbi={\"67.25\": 1},\n",
    "    nbi_err=0.3,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "avr_rspectrum, avr_rspectrum_w, diffs = get_averaged_spectra_and_diffs(mock_traces_DF.iloc[:,2:].values)\n",
    "diffs2 = diffs\n",
    "label4final_histo.append('mock traces')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_results(avr_rspectrum, avr_rspectrum_w, diffs, histo_edge=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mock_traces_with_BB_DF = mock_traces_DF.copy(deep=True)\n",
    "for i in range(mock_traces_with_BB_DF.index.size):\n",
    "    start_bb_index = np.random.randint(0, 2048-40)\n",
    "    time_trace = mock_traces_with_BB_DF.iloc[i,2:].values\n",
    "    time_trace[0:bb_tt.size] = time_trace[0:bb_tt.size] + bb_tt\n",
    "    time_trace = np.roll(time_trace_e, random.randint(0, 2048 - bb_tt.size) )\n",
    "    mock_traces_with_BB_DF.iloc[i,2:] = time_trace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "avr_rspectrum, avr_rspectrum_w, diffs = get_averaged_spectra_and_diffs(mock_traces_with_BB_DF.iloc[:,2:].values)\n",
    "diffs3 = diffs\n",
    "label4final_histo.append('mock traces with BB pulse')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_results(avr_rspectrum, avr_rspectrum_w, diffs, histo_edge=3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "fx = np.fft.rfftfreq(N, 1 / fs_Hz) / 1e6\n",
    "\n",
    "diffs_list = np.asarray([diffs1, diffs2, diffs3], dtype=\"object\") - 1\n",
    "# diffs_list = [diffs1, diffs2]\n",
    "\n",
    "bins = np.histogram_bin_edges(np.concatenate(diffs_list), bins=20)\n",
    "\n",
    "\n",
    "# if bins == None:\n",
    "bins = linspace_with_middle_value(\n",
    "    np.mean(np.concatenate(diffs_list)), 2 * np.std(np.concatenate(diffs_list)), 20\n",
    ")\n",
    "\n",
    "bins = np.linspace(0.5, 1.5, 100) - 1\n",
    "\n",
    "for i, diffs in enumerate(diffs_list):\n",
    "    mean_diff = np.mean(diffs)\n",
    "    std_diff = np.std(diffs)\n",
    "    mean_diff_trun, std_diff_trun = calculate_truncated_stats(diffs, 1, 99)\n",
    "    label = (\n",
    "        text_box\n",
    "    ) = f\"$\\mu_{{trunc}}$: {round(mean_diff_trun,2)+0:.2f}\\n$\\sigma_{{trunc}}$: {round(std_diff_trun,2)+0:.2f}\"\n",
    "    ax.hist(diffs, bins=bins, alpha=1, density=True, histtype=\"step\", lw=3, label=label)\n",
    "    # Calculate mean and standard deviation\n",
    "    print(label4final_histo[i])\n",
    "    print('Mean and STD:{:.4f}, {:.4f}'.format(mean_diff, std_diff))\n",
    "    print('Truncated Mean and STD: {:.4f}, {:.4f}'.format(mean_diff, std_diff))\n",
    "    print(\"******\")\n",
    "\n",
    "# ax.set_xlim(xax_min, xax_max)\n",
    "ax.set_xlabel(\n",
    "    r\"$\\frac{\\mathrm{energy: \\ spectrum \\ with \\ no\\ window}}{\\mathrm{energy:spectrum \\ with\\ window,\\ amplitudes\\ corrected}} -1$\"\n",
    ")\n",
    "ax.set_ylabel(\"entries\")\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('The colors corespond to')\n",
    "print(label4final_histo)\n",
    "print(', respectively.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CONCLUSION: The broad band pulse in the trace spoils the energy recovery."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}