Fit_lifetime_2022.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3.5",
      "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.5.3"
    },
    "colab": {
      "provenance": [],
      "include_colab_link": true
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/mbarbetti/unifi-physics-lab3/blob/main/Fit_lifetime_2022.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ROGuWpjT97Ed"
      },
      "source": [
        "## Misura della vita media del muone\n",
        "Lo scopo di questo programma è misurare la vita media del muone facendo un fit hai dati raccolti in laboratorio\n",
        "Il primo blocco qui sotto serve solo a importare (e installare quelli che non sono disponibili) alcuni pacchetti che ci sono necessari e in particolare Minuit, il programma di minimizzazione, che si trova nella libreria iminuit, che deve essere installata con pip perché non è presente tra le librerie disponibili di default in Colab.\n",
        "L'ultima riga serve per ridefinire le dimensioni dei plots."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BqM5F5EACp6G"
      },
      "source": [
        "import pandas as pd\n",
        "from urllib.request import urlopen\n",
        "import sys\n",
        "import os\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.optimize import curve_fit\n",
        "from scipy.stats import poisson\n",
        "from scipy.stats import chi2\n",
        "!pip install -q iminuit\n",
        "from iminuit import Minuit\n",
        "\n",
        "\n",
        "plt.rcParams['figure.figsize'] = [16, 8]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xer78ISteVAI"
      },
      "source": [
        "Prima di tutto dovete **modificare il codice inserendo il nome del gruppo**, in modo che il programma trovi il file con i dati. Le istruzioni che seguono stampano l'header del file con alcune informazioni che ci serviranno in seguito. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CDXHdN4Qk9Ia"
      },
      "source": [
        "La prima operazione necessaria è importare i dati nel file scritto da LabView in formato pandas. Dopo questa operazione i nostri dati si trovano nell'oggetto chiamato **data** che è di tipo DataFrame (una tabella).\n",
        "Alla fine si moltiplicano i valori TDC per 1e09 per trasformarli in ns. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4emZiP_vtn_z"
      },
      "source": [
        "## INSERIRE NUMERO DEL GRUPPO\n",
        "group = \"M2\"   # esempio\n",
        "\n",
        "### SET THE VALUE OF A AND B\n",
        "a = 0.97693 #ns/a.u.\n",
        "b = -90.14805 #ns\n",
        "    \n",
        "header = urlopen (f\"https://raw.githubusercontent.com/mbarbetti/unifi-physics-lab3/main/data/2022/data_{group}.txt\")\n",
        "for i, line in enumerate(header):\n",
        "  print (line.decode('utf-8'))\n",
        "  if i == 5: break\n",
        "\n",
        "print ('Calibration factors:')\n",
        "print ('    a = {:.5f} ns/a.u.' . format (a))\n",
        "print ('    b = {:.5f} ns' . format (b))\n",
        "print ('\\nt_corr = ({:.5f}) * t_mis + ({:.5f}) ns' . format (a, b))\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Xwkpl9B5Cp6p"
      },
      "source": [
        "file_path = f\"https://raw.githubusercontent.com/mbarbetti/unifi-physics-lab3/main/data/2022/data_{group}.txt\"\n",
        "data = pd.read_csv (file_path, header = 4, delim_whitespace = True).drop([0])\n",
        "data['QDC1'] = -data['QDC1'] * 1e9\n",
        "data['QDC2'] = -data['QDC2'] * 1e9\n",
        "data['TDC'] = data['TDC'] * 1e9\n",
        "print(data)   # dataframe\n",
        "print(min(data['QDC1']),\" \",min(data['QDC1']))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xJW4VN7zq3Su"
      },
      "source": [
        "Qui si riempiono degli istogrammi bi- e uni-dimensionali con i valori del QDC1 e QDC2 (rispettivamente muone ed **elettrone**). Le unità sono arbitrarie e sono moltiplicate per 1e09 per comodità di visualizzazione. Queste distribuzioni servono per selezionare gli eventi di fisica rispetto a quelli di calibrazione: guardando allo scatter plot **dovete decidere dove mettere i tagli per selezionare gli eventi corrispondenti ai decadimenti dei muoni**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "g10znHKvoDBJ"
      },
      "source": [
        "### Plot qdc values\n",
        "### SET THE QDC RANGES FOR MUON DECAY SIGNAL \n",
        "qdc1_min = 3.5\n",
        "qdc1_max = 10 \n",
        "qdc2_min = 3.4\n",
        "qdc2_max = 10\n",
        "\n",
        "data_muon = data.query (f\"(QDC1 > {qdc1_min}) & (QDC1 < {qdc1_max}) & (QDC2 > {qdc2_min}) & (QDC2 < {qdc2_max})\")\n",
        "qdc1 = data_muon['QDC1']\n",
        "qdc2 = data_muon['QDC2']\n",
        "qdc1_all = data['QDC1']\n",
        "qdc2_all = data['QDC2']\n",
        "\n",
        "plt.xlabel (\"QDC1\")\n",
        "plt.ylabel (\"QDC2\")\n",
        "plt.plot ([0,10], [qdc2_min, qdc2_min], color = \"red\", lw = 2)\n",
        "plt.plot ([qdc1_min, qdc1_min], [0,10], color = \"red\", lw = 2)\n",
        "plt.scatter (qdc1_all, qdc2_all, s=1)\n",
        "plt.axis ([2,10,2,8])\n",
        "plt.show()\n",
        "\n",
        "nbin = 100\n",
        "\n",
        "plt.xlabel (\"QDC1\")\n",
        "plt.ylabel (\"Entries\")\n",
        "plt.hist (qdc1_all, bins=nbin, range = (2,10), alpha = 0.75, label='QDC1')\n",
        "plt.hist (qdc1, bins=nbin, range = (2,10), alpha = 0.75, label='QDC1')\n",
        "plt.yscale(\"log\")\n",
        "plt.show()\n",
        "\n",
        "plt.xlabel (\"QDC2\")\n",
        "plt.ylabel (\"Entries\")\n",
        "plt.hist (qdc2_all, bins=nbin, range = (2,10), alpha = 0.75, label='QDC2')\n",
        "plt.hist (qdc2, bins=nbin, range = (2,10), alpha = 0.75, label='QDC2')\n",
        "plt.yscale(\"log\")\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DnCHcpB0lMQq"
      },
      "source": [
        "Si riempie adesso l'istogramma con i valori dei tempi. Questo istogramma è poi usato per il fit, quindi si deve prestare attenzione a definire **tmin** e **tmax** in modo da eliminare eventuali valori non fisici e le regioni vicine ai valori limite di accettanza della nostra misura.  \n",
        "Scegliete il numero di bin **nbins** e di conseguenza la larghezza del bin **delta_t**. La larghezza del bin deve essere scelta in modo che la densità di probabilità non cambi significativamente nell'intervallo del bin, ma anche che sia abbastanza grande perché si possa applicare il test del $\\chi^2$ di Pearson (ovvero i conteggi devono essere > 10). \n",
        "Per ingrandire una zona dell'istogramma potete modificate gli estremi dell'asse x con **plt.xlim(right,left)** e la scala verticale può essere impostata in modo logaritmico con **plt.yscale(\"log\")**\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dq0rT9ImCp7T"
      },
      "source": [
        "### Prepare and fill histogram of time intervals\n",
        "### SET TMIN, TMAX AND NBINS\n",
        "tmin = 100\n",
        "tmax = 10000\n",
        "nbins = 200\n",
        "\n",
        "delta_t = (tmax - tmin) / nbins\n",
        "\n",
        "### Apply calibration factors\n",
        "\n",
        "tdc = a*data_muon['TDC'].to_numpy() + b\n",
        "\n",
        "histo_times = np.histogram(tdc, bins = nbins, range = (tmin, tmax))\n",
        "#print(histo_times)\n",
        "\n",
        "plt.xlabel('Time interval [ns]')\n",
        "plt.ylabel('Number of events')\n",
        "plt.hist (tdc, bins = nbins, range = (tmin, tmax), label='Measured time intervals')\n",
        "plt.yscale(\"linear\") #linear or log\n",
        "plt.show()\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "V8KBR51FrQMJ"
      },
      "source": [
        "Si riempiono adesso i vettori con i valori di x e y che saranno usati per il fit"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Oi4j0ZGyCp7d"
      },
      "source": [
        "### Prepare data for the fit\n",
        "# histo_times #print the bin values and edges \n",
        "### Super cool way to get bin centers\n",
        "x_values = (histo_times[1][1:] + histo_times[1][:-1]) / 2\n",
        "y_values = histo_times[0]\n",
        "#bin_edges = histo_times[1].tolist()\n",
        "print(x_values)\n",
        "print(y_values)\n",
        "print (delta_t)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rK9LY23trcw4"
      },
      "source": [
        "Si costruisce la negative log-likelihood (NLL) considerando il numero di eventi in ogni bin come il risultato di una distribuzione di poisson con valore medio dato da $ \\nu_i = \\nu_{tot}\\int_{x_i}^{x_{i+1}} \\frac{1}{\\tau}\\exp(-\\frac{x}{\\tau}) + A\\Delta t \\approx \\nu_{tot}\\frac{1}{\\tau}\\exp(-\\frac{x}{\\tau})\\Delta t + A\\Delta t $\n",
        "\n",
        "L'ultimo termine $A\\Delta t$ rappresenta il numero di coincidenze casuali in ogni bin, con il parametro $ A = R^2\\cdot T_M \\cdot 10^{-9} \\, [ns^{-1}]$, dove $ R $ è il ritmo di arrivo dei muoni in *Hz* e $ T_M $ il tempo totale di misura in *s*. Il fattore $10^{-9}$ tiene conto che $\\Delta t$ è espresso in *ns*. \n",
        "\n",
        "Il fattore $A$ è posto inizialmente a 0 ed è fissato nel fit dall'opzione di Minuit *fix_A = True*. Successivamente potete provare a lasciarlo libero nel fit o a fissarlo al valore calcolato tramite la formula qui sopra. \n",
        "Al posto della semplice esponenziale si può usare la funzione *fullExp()* che tiene conto dell'assorbimento dei muoni negativi.  \n",
        "Le varie likelihood che possono essere usate nel fit sono: \n",
        "\n",
        "1) *NLLbin()*: binned extended maximum likelihood con la funzione *simpleExp()* per i decadimenti, che non tiene conto dell'assorbimento dei muoni\n",
        "\n",
        "2) *NLLass()*: come sopra ma con *fullExp()* per i decadimenti, che tiene conto dell'assordimento\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ldHQ1ThRCp7j"
      },
      "source": [
        "### Build Negative Log Likelihood (NLL) to minimize\n",
        "\n",
        "def simpleExp(t, nu, tau):\n",
        "\n",
        "    return nu*(np.exp(-t/tau)/tau)\n",
        "\n",
        "def fullExp(t, nu, tau):\n",
        "\n",
        "    r = 1.27 \n",
        "    lass = 38.8e-6\n",
        "\n",
        "    return (nu/tau)/(1+r)*(r+np.exp(-lass*t))*np.exp(-t/tau)\n",
        "\n",
        "def NLLbin(nu, tau, A):\n",
        "    \n",
        "    nlogL = 0\n",
        "\n",
        "    for i, count in enumerate(y_values):\n",
        "        nlogL = nlogL - poisson.logpmf(count, (simpleExp(x_values[i], nu, tau) + A)*delta_t )\n",
        "    \n",
        "    return nlogL\n",
        "\n",
        "\n",
        "def NLLass(nu, tau, A):\n",
        "    \n",
        "    nlogL = 0\n",
        "\n",
        "    for i, count in enumerate(y_values):\n",
        "        nlogL = nlogL - poisson.logpmf(count, (fullExp(x_values[i], nu, tau) + A)*delta_t)\n",
        "    \n",
        "    return nlogL\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_rQ85MZWtOB7"
      },
      "source": [
        "Adesso si inizializza *Minuit* passando la funzione da minimizzare e il valore dei vari parametri. Successivamente si chiama *migrad()* per minimizzare NLL e poi *hesse()* per ri-calcolare più precisamente la matrice di covarianza e infine *minos()* per calcolare gli errori. Per chi vuole approfondire, la documentazione di *Minuit* è disponibile al link https://iminuit.readthedocs.io/en/stable/\n",
        "\n",
        "**Il primo argomento è la funzione di fit:** si possono indicare *NLLbin*, *NLLass*. Fate riferimento al blocco precedente per la descrizione.\n",
        "\n",
        "**I parametri possono essere lasciati liberi o fissati nel fit** impostando m.fixed['nome del parametro'] rispettivamente su *False* o *True*. Per esempio il parametro *A* è inizialmente fissato.  "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OSMejb5mCp7p"
      },
      "source": [
        "### Call Migrad to find minimum of the NLL\n",
        "### SET NLL FUNCTION TYPE AS NLLbin or NLLass\n",
        "NLL = NLLass\n",
        "m = Minuit(NLL, nu=4000, tau=2200., A = 0.0)\n",
        "m.errordef = Minuit.LIKELIHOOD     # == 0.5\n",
        "m.fixed['A'] = True\n",
        "m.migrad()\n",
        "m.hesse()\n",
        "m.minos()\n",
        "m.migrad()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "U_aUZlCzyM_8"
      },
      "source": [
        "Stampiamo i parametri ottenuti dal fit in una forma più leggibile"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mJ7uFP9YCp75"
      },
      "source": [
        "### Get fit results and errors from Minuit\n",
        "nu_fit = m.values['nu']\n",
        "tau_fit = m.values['tau']\n",
        "A_fit = m.values['A']\n",
        "\n",
        "nu_fit_err = m.errors['nu']\n",
        "tau_fit_err = m.errors['tau']\n",
        "A_fit_err = m.errors['A']\n",
        "\n",
        "### Print the results\n",
        "if m.fixed['A']:\n",
        "  print(f'Parameters from fit:\\n\\tnu_tot = {nu_fit:.0f} +- {nu_fit_err:.0f} \\n\\tTau = {tau_fit:.2f} +- {tau_fit_err:.2f} ns\\n\\tA = {A_fit:.3f} 1/ns  (fixed)')\n",
        "else:\n",
        "  print(f'Parameters from fit:\\n\\tnu_tot = {nu_fit:.0f} +- {nu_fit_err:.0f} ns^-1\\n\\tTau = {tau_fit:.2f} +- {tau_fit_err:.2f} ns\\n\\tA = {A_fit:.3f} +- {A_fit_err:.3f} 1/ns')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "t63Zex1B7TDX"
      },
      "source": [
        "Confrontiamo il risultato ottenuto con la distribuzione nei dati. Possiamo anche calcolare il $\\chi^2$. \n",
        "\n",
        "**Si noti che il numero di gradi di libertà *ndf* è definito tenendo conto del numero di parametri liberi nel fit che si ottiene con m.nfit**\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mkQT9WdDCp8E"
      },
      "source": [
        "### Plot both the histograms as well as the fitted distribution\n",
        "\n",
        "if NLL != NLLass:\n",
        "  pdf = simpleExp\n",
        "else:\n",
        "  pdf = fullExp\n",
        "\n",
        "t_plot = np.arange (0, 10000, 10)\n",
        "y_plot = (pdf(t_plot, nu_fit, tau_fit) + A_fit)*delta_t\n",
        "plt.plot(t_plot, y_plot, label='Fitted distribution', color = \"red\", zorder = 1, lw = 2)\n",
        "\n",
        "y_errors = np.sqrt(y_values)\n",
        "plt.errorbar(x_values, y_values, yerr=y_errors, elinewidth=2, linewidth=0, color='blue', markerfacecolor='blue', markersize=4, markeredgewidth=3, marker='v', label='Data points', zorder = 0)\n",
        "plt.legend()\n",
        "plt.show()\n",
        "\n",
        "### Calculate chi2 and print fit results\n",
        "y_fit = (pdf(x_values, nu_fit, tau_fit) + A_fit)*delta_t\n",
        "residuals = y_values - y_fit\n",
        "squares = np.square(residuals)/y_fit\n",
        "chi2fit = squares.sum()   ## chi2 value\n",
        "ndf = len(squares) - m.nfit ## number of degree of freedom\n",
        "p_val = 1 - chi2.cdf (chi2fit, ndf)\n",
        "\n",
        "print(f'Chi square of the fitted distribution:\\n\\tchi2 = {chi2fit:.4f}\\n\\tndf = {ndf:.4f}\\n\\tp-value = {p_val:.4f}')\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ndSTvUvKZwHN"
      },
      "source": [
        "Usiamo la funzione *draw_mnprofile()* per disegnare il profilo della negative log-likelihood (NLL) in un intorno del minimo al variare di $\\tau$. \n",
        "\n",
        "Si noti che per ogni valore di tau vengono ricalcolati i valori degli altri parametri che minimizzano la NLL: questo equivale a considerare la densità di probabilità marginale per il parametro per cui si disegna il profilo.\n",
        "\n",
        "La banda mostrata corrisponde ad una variazione di $\\pm 1 \\sigma$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ilHXyMbPZnyd"
      },
      "source": [
        "m.draw_mnprofile('tau')\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "la8e-iIYzlpK"
      },
      "source": [
        "Se il parametro *A* non è fissato, usiamo la funzione *draw_mncontour()* per disegnare le linee di livello corrispondenti ad un livello di confidenza pari al 39% e al 68%. Per curiosità potete provare a vedere come viene il contorno anche per *nu* verso *tau* o verso *A*. \n",
        "\n",
        "**Sapreste dire a quale variazione del valore della negative log-likelihood corrispondono queste linee di livello?**"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "B01kBW_mCp8K"
      },
      "source": [
        "### Ask Minuit for the 2D contour at CL = 39% and CL = 68% if A is not fixed\n",
        "if not m.fixed['A']:  \n",
        "  m.draw_mncontour('tau','A', cl=[0.39,0.68])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BXR4nKU-DDRn"
      },
      "source": [
        "**Passi successivi**\n",
        "\n",
        "1) Adesso provate a lasciare libero il parametro A. Che valori ottenete? Come cambia il chi2? E l'errore su *tau*? Provate a guardare l'ellisse di correlazione tra *A* e *tau*: sapete spiegarvi il risultato? \n",
        "\n",
        "2) Provate adesso a tenere conto dell'effetto dell'assorbimento dei muoni negativi. Il valore trovato è compatibile con quello noto in letteratura $\\tau = 2197$ ns? \n",
        "\n",
        "3) Provate infine a calcolare A utilizzando il valore di $R$ misurato in laboratorio. \n",
        "\n",
        "Utilizzate il tempo $T_M$ della misura che è scritto nell'header del file con i dati e che abbiamo stampato all'inizio. Ricordate di sottrarvi un numero di secondi pari numero di eventi raccolti per tener conto del tempo morto dovuto all'acquiaizione. Questo tempo morto può essere considerato come l'errore sistematico du $T_M$.\n",
        "\n",
        "Dopo tutte queste considerazioni qual è l'errore sulla vostra stima di A? \n",
        "Il risultato ottenuto è compatibile con quello ottenuto dal fit al punto precente (con NLLass)?\n",
        "\n",
        "Per tener conto di questa incertezza ripetete il fit (sempre con NLLass) modificando A di più o meno l'errore. La variazione nel risultato del fit per $\\tau$ è l'errore sistematico dovuto ad A. Quanto è importante se confrontato con l'errore statistico del fit? \n",
        "\n",
        "Quanto vale l'errore totale sulla misura? Confrontatelo con quello che avete ottenuto al punto 2) e commentate la differenza. Per una discussione finale e per tirare le vostre conclusioni tenete anche conto dei valori del chi2 (e dei relativi p-valori) che avete ottenuto\n"
      ]
    }
  ]
}