MuonLifetimeMeaurement_school.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f161080",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import ROOT\n",
    "from iminuit import Minuit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00c4f996",
   "metadata": {},
   "source": [
    "Part1: Create exponential random numbers for a given value of tau and plot it ;\n",
    "Exponential pdf = 1/τ*e^(-t/τ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6446b36-1ed0-4d93-b142-130455003e08",
   "metadata": {},
   "outputs": [],
   "source": [
    "##transforming U(0,1) to exponential random number t\n",
    "def generate_time(n1,tau):\n",
    "    root_time_measurements = []\n",
    "    #TODO:generate exponential random numbers->store in root_time_measurements\n",
    "    return root_time_measurements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78c94f53-4ca7-4496-9511-82e71c4ce332",
   "metadata": {},
   "outputs": [],
   "source": [
    "##generating exponential random number t using numpy\n",
    "def generate_time_numpy(n1,tau):\n",
    "    #TODO:generate random number using numpy->store in numpy_time_measurements\n",
    "    return numpy_time_measurements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35314008",
   "metadata": {},
   "outputs": [],
   "source": [
    "##filling tree with the time measuremnets\n",
    "def fill_tree_with_time(outfile, time_measurements):\n",
    "    tree = ROOT.TTree(\"time_measurements\", \"Time Measurement Data\")\n",
    "    #TODO: Fill the tree branch with time measurements from outfile\n",
    "    tree.Write()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01c5c6d2-089d-4fa7-9a08-670ff50a5fbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "##plot the time measurements\n",
    "def plot_time_measurements(outfile):\n",
    "    tree = outfile.Get(\"time_measurements\")\n",
    "    x_data = np.array([entry.x for entry in tree])\n",
    "    plt.figure(figsize=(8, 6))\n",
    "    plt.hist(x_data, bins=50, density=True, alpha=0.7, label=\"Generated Time Measurement Data\")\n",
    "    plt.xlabel(\"time\")\n",
    "    plt.ylabel(\"\")\n",
    "    plt.legend()\n",
    "    plt.title(\"Generated Exponential Distribution\")\n",
    "    plt.show()\n",
    "    plt.close('all')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02597319",
   "metadata": {},
   "outputs": [],
   "source": [
    "outfile = ROOT.TFile(\"DecayTimeMeasurements.root\", \"recreate\")\n",
    "n1 = 1000 #number of time measurements in a particular experiment\n",
    "tau = 5 #mean muon lifetime\n",
    "time_measurements = generate_time(n1,tau) \n",
    "#time_measurements = generate_time_numpy(n1,tau) #generating exponential random number t using numpy\n",
    "fill_tree_with_time(outfile, time_measurements) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1dfc80b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_time_measurements(outfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c92d6940",
   "metadata": {},
   "source": [
    "Part2: Create exponential random numbers for a given value of tau and plot it ;\n",
    "Exponential pdf = 1/τ*e^(-t/τ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6381eda2",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Unbinned NLL fit to the generated time_measurements data using Minuit\n",
    "def fit_time_measurements(outfile):\n",
    "    tree = outfile.Get(\"time_measurements\")\n",
    "    def fcn(tau):                                                                                                                  \n",
    "        time = [entry.x for entry in tree]                                                                                         \n",
    "        lnL = 0.0                                                                                                                  \n",
    "        for x in time:                                                                                                             \n",
    "          #TODO: construct the log likelihood                                                                                      \n",
    "        return -2.0*lnL     \n",
    "    minuit = Minuit(fcn, tau=2)\n",
    "    minuit.limits[\"tau\"] = (2, 10)\n",
    "    minuit.migrad()\n",
    "    minuit.hesse()\n",
    "    tau_expected = minuit.values[\"tau\"]\n",
    "    error_tau_expected = minuit.errors[\"tau\"]\n",
    "    print(\"Expected tau:\", tau_expected, \"+/-\", error_tau_expected)\n",
    "    return tau_expected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc50a735",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot the fit and the data\n",
    "def plot_time_measurement_fit(outfile, tau_expected):\n",
    "    tree = outfile.Get(\"time_measurements\")\n",
    "    x_data = np.array([entry.x for entry in tree])\n",
    "    x_fit = np.linspace(0, max(x_data), 1000)\n",
    "    y_fit = np.exp(-x_fit/tau_expected)/tau_expected\n",
    "    plt.figure(figsize=(8, 6))\n",
    "    plt.hist(x_data, bins=50, density=True, alpha=0.7, label=\"Generated Time Measurement Data\")\n",
    "    plt.plot(x_fit, y_fit, label=\"Fitted Exponential Distribution\", color=\"red\", linewidth=2)\n",
    "    plt.xlabel(\"time\")\n",
    "    plt.ylabel(\"\")\n",
    "    plt.legend()\n",
    "    plt.text(\n",
    "        0.6,\n",
    "        0.8,\n",
    "        f'Fitted tau: {tau_expected:.4f}',\n",
    "        transform=plt.gca().transAxes,\n",
    "        bbox=dict(facecolor='white', alpha=0.8),\n",
    "    )\n",
    "    plt.title(\"Fitted Exponential Distribution\")\n",
    "    plt.show()\n",
    "    plt.close('all')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6f2ecce",
   "metadata": {},
   "outputs": [],
   "source": [
    " # Fit distribution and then plot the fit along with data\n",
    "tau_expected = fit_time_measurements(outfile)\n",
    "plot_time_measurement_fit(outfile, tau_expected)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f19f71f7",
   "metadata": {},
   "source": [
    "Part 3: Do n2 such experiments and plot the fitted taus ; Please have a look at the mean of the gaussian \n",
    "The Central Limit Theorem:\n",
    "For samples of size  30 or more, the sample mean is approximately normally distributed, with mean μ  and standard deviation = σ/√n; where σ = std of the sample=tau, n  is the sample size. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37f4fa23",
   "metadata": {},
   "outputs": [],
   "source": [
    "n2=1000 #total number of experiments\n",
    "tau_expected_samples = []\n",
    "print(\"Fitting the samples\")\n",
    "#TODO: Do n2 such experiments: for each experiment fit the n1 random number's  distribution->calculate the expected value of tau->store the expected values in tau_expected_sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60ab5005",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_tau_expected(tau_expected_samples):\n",
    "    plt.figure(figsize=(12, 6))\n",
    "    plt.hist(tau_expected_samples, bins=25, density=True, alpha=0.7, label=\"tau expected\")\n",
    "    plt.xlabel(\"tau expected\")\n",
    "    plt.ylabel(\" \")\n",
    "    plt.title(\"Histogram of expected tau in many experiments\")\n",
    "    plt.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0593edeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_tau_expected(tau_expected_samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71e7b6c9",
   "metadata": {},
   "source": [
    "Part 4: From the distribution of expected tau, calculate the std and match with the expected (σ/√n1) one from CLT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d5ab115",
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_error(tau_expected_samples,n1,n2,tau):\n",
    "    #TODO:Find the Error on tau expected from tau_expected_samples -> tau_expected_std\n",
    "    time_measurements_std = tau/(np.sqrt(n1))\n",
    "    print(f\"Error on the tau expected: {tau_expected_std:.4f}\")\n",
    "    print(f\"Expected std : {time_measurements_std:.4f}\")\n",
    "    print(f\"Ratio: {tau_expected_std / time_measurements_std:.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6f7b1ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "check_error(tau_expected_samples, n1, n2, tau)"
   ]
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}