TS_DegFreedom_Sigma_Relation_FermiLAT.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "K6w0AljmdSpB"
      },
      "source": [
        "# Fermi-LAT: Convert TS + Degree of Freedom --> $\\sigma$  \n",
        "\n",
        "This short routine is intended as a tool to calculate the  statistical significance (\\sigma) of a $\\gamma$-ray signature, detected with the Fermi Large Area Telescope (LAT). The sigma value depends on the number of Degrees of Freedom (dof) involved in a Likelihood Analysis and the resulting Test Statistic (TS) value. \n",
        "\n",
        " \n",
        "**A standard analysis with Fermi-LAT** involves four degrees of freedom. Two associated to the uncertain position of the source (R.A.,Dec.), and another two associated with the shape of the $\\gamma$-ray spectrum that -in the simplest case- can be assumed to follow a power-law, $dN/dE = N_0 (E/E_0)^\\Gamma $. Here, $N_0$ is the pre-factor or normalization parameter [photons/cm^2/s/MeV], $\\Gamma$ is the photon index, and $E_0$ is the pivot-energy [MeV] where both $N_0$ and $\\Gamma$ are estimated. Usually, the $E_0$ value has to be fixed during the analysis, otherwise, the number of degrees of freedom brings far too much uncertainty to the $N_0$, $\\Gamma$ and $E_0$ values. Still, one should run the analysis testing over a reasoable range of $E_0$ values (from 500 MeV up to Tens of GeV) to optimize the $N_0$ and $\\Gamma$ estimates. This is the typical approach to Fermi-LAT data, if the analysis is solo based in $\\gamma$-ray information (blind to information collected from other bands).\n",
        "\n",
        "In this case, a signature with TS = 25, from a likelihood analysis with 4 degrees of freedom, corresponds to a 4.05$\\sigma$ detection. \n",
        "\n",
        "**Multifrequency Information:** However, it is possible to take into account multifrequency information, and lower the uncertainty of the analysis, by fixing the position of the $\\gamma$-ray emitter. Please check the [2BIGB $\\gamma$-ray catalogue](https://arxiv.org/abs/1911.08912). This is perfectly fine if one considers the position of Blazars -especially out of the galactic plane- because blazars are the dominant extragalactic $\\gamma$-ray sources. \n",
        "\n",
        "In this case, if the position of the $\\gamma$-ray candidate is known, there are only two degrees of freedom for the likelihood analysis: the ones associated with the power-law spectral modeling. \n",
        "\n",
        "A signature with TS = 19.8, deduced from a likelihood analysis with only 2 degrees of freedom, corresponds to a 4.05$\\sigma$ detection. Therefore, if one searches for gamma-ray signatures associated with blazars (as for the 2BIGB catalogue), a TS=19.8 threshold is equivalent to a TS=25 blind search detection threshold. \n",
        "\n",
        "A signature with TS = 10, deduced from a likelihood analysis with only 2 degrees of freedom, corresponds to a 2.71$\\sigma$ detection. \n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        },
        "colab_type": "code",
        "id": "OotHKgSvCkVc",
        "outputId": "8c01d217-48ce-4737-c2a4-2692b763a0ca"
      },
      "outputs": [],
      "source": [
        "# import libraries\n",
        "import numpy as np\n",
        "from scipy.stats import norm, chi2\n",
        "\n",
        "# define sigma function\n",
        "def sigma(TS, dof):\n",
        "    \"\"\"\n",
        "    Calculate the significance (in sigma) from a Test Statistic (TS) value.\n",
        "\n",
        "    Parameters:\n",
        "    TS (float): Test Statistic value, must be non-negative.\n",
        "    dof (int): Degrees of Freedom, must be positive.\n",
        "\n",
        "    Returns:\n",
        "    float: The significance in terms of sigma.\n",
        "    \n",
        "    Raises:\n",
        "    ValueError: If TS is negative or dof is not positive.\n",
        "    \"\"\"\n",
        "    # Validate inputs\n",
        "    if TS < 0:\n",
        "        raise ValueError(\"TS value must be non-negative.\")\n",
        "    if dof <= 0:\n",
        "        raise ValueError(\"Degrees of Freedom (dof) must be positive.\")\n",
        "    if not isinstance(dof, int):\n",
        "        raise ValueError(\"Degrees of Freedom (dof) must be integer.\") \n",
        "\n",
        "    # Compute sigma\n",
        "    if TS > 1400:\n",
        "        # Approximation for large TS\n",
        "        nsigma = np.sqrt(TS)\n",
        "    else:\n",
        "        # General case\n",
        "        nsigma = norm.isf(chi2.sf(TS, dof) / 2)\n",
        "    \n",
        "    return nsigma\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 34,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "For a likelihood analysis with 4 degrees of freedom,\n",
            " a signature with TS 25 corresponds to a 4.05 sigma detection.\n",
            "For a likelihood analysis with 2 degrees of freedom,\n",
            " a signature with TS 19.8 corresponds to a 4.05 sigma detection.\n",
            "For a likelihood analysis with 2 degrees of freedom,\n",
            " a signature with TS 10 corresponds to a 2.71 sigma detection.\n"
          ]
        }
      ],
      "source": [
        "# Example with TS=25 and 4 dof \n",
        "TS=25\n",
        "dof=4\n",
        "significance = sigma(TS,dof)\n",
        "print(f\"For a likelihood analysis with {dof} degrees of freedom,\\n a signature with TS {TS} corresponds to a {significance:.2f} sigma detection.\")\n",
        "\n",
        "# Example with TS=25 and 2 dof \n",
        "TS=19.8\n",
        "dof=2\n",
        "significance = sigma(TS,dof)\n",
        "print(f\"For a likelihood analysis with {dof} degrees of freedom,\\n a signature with TS {TS} corresponds to a {significance:.2f} sigma detection.\")\n",
        "\n",
        "# Example with TS=10 and 2 dof \n",
        "TS=10\n",
        "dof=2\n",
        "significance = sigma(TS,dof)\n",
        "print(f\"For a likelihood analysis with {dof} degrees of freedom,\\n a signature with TS {TS} corresponds to a {significance:.2f} sigma detection.\")\n"
      ]
    }
  ],
  "metadata": {
    "accelerator": "GPU",
    "colab": {
      "collapsed_sections": [],
      "name": "TS-DegFreedom-Sigma-Relation-FermiLAT.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}