"The gammapy package is used to compute Feldman-Cousins statistics\n",
"\n",
"## Prerequisits\n",
"`pip install gammapy, matplotlib, openkm3`"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"import pandas as pd\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import math\n",
"from scipy import stats\n",
"import gammapy.stats as gstats\n",
"import openkm3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Helper functions\n",
"\n",
"### Background evaluation\n",
"The function *get_number_background(decl, roi)* permits to evaluate the expected background for a give sky declination and a given region of interest ('half width' of the declination band around the source declination). The right ascension is not used here, as the background distribution depends only on the visibility of the given sky position from the location of ANTARES, and therefore on the declination only.\n",
"\n",
"This function uses a polynomial interpolation to smoothen the background distibution otherwise obtained from the data. This function returns the number of expected background events that in the following will be indicated as nbkg."
"The functions *get_upper_limit(nbkg, nobs)* and *get_lower_limit(nbkg, nobs)* compute the upper and lower limits in presence of nbkg backgroud events, if the measurements is nobs observed events inside the region of interest."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"# compute upper limit according to Feldamn-Cousins\n",
"\n",
"def get_upper_limit(nbkg, nobs):\n",
" \n",
" if nobs <= int(nbkg):\n",
" ul = Sens(nbkg) #if underfluctuation, UL set as the sensitivity \n",
"The function *get_sensitivity(nbkg)* is used to calculate the sensitivity, defined as average upper limit, that is reached in presence of nbkg events in the declination band of the source."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"confidence_level = 0.9 # set confidence level to 90% \n",
" # count events whose angle from the source falls inside the RoI\n",
" if alpha < RoI:\n",
" obs = obs+1\n",
" return obs"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Math helpers\n",
"\n",
"The two functions *make_vector(Decl, RA)* and *calc_angle(v1, v2)* are auxilary functions used later. \n",
"\n",
"*make_vector* fills a 3-dimensional vector using the sky coordinates of an event (Decl, RA), and assigning length 1. \n",
"\n",
"*calc_angle* computes the angular distance (in units degrees) between two points in the sky defined as the angle between their corresponding vectors v1 and v2."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"def make_vector(Decl, RA): \n",
" x = math.cos(math.radians(Decl))*math.cos(math.radians(RA))\n",
" y = math.cos(math.radians(Decl))*math.sin(math.radians(RA))\n",
"Here, the flux upper limits for the observation of ANTARES events from the 2007-2017 sample for a given direction in the sky is calculated.\n",
"\n",
"### Setting the point source coordinates\n",
"The limits are calculated providing \n",
"* the source coordinates **right ascension** and **declination**, ra_source and dec_source, and \n",
"* the two additional source parameters **gamma** (spectral index) and \n",
"* **RoI** (region of interest, representing the source extension)."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"ra_source = 0.1 # between 0 and 360\n",
"dec_source = -29 # Declination must be between -80 and 50 degree\n",
"gamma = 2.0 # spectral index between 1.5 and 3.0\n",
"RoI = 2.6 # Region of interest > 0.1 and < 1000 (in degree)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Getting ANTARES data and acceptance\n",
"This is a proposal of how retrieving the relevant ANTARES data should work.\n",
"The KM3Store should offer\n",
"* the event list\n",
"* the acceptance function\n",
"* the background estimate from interpolation\n",
"\n",
"The following functions are now dummies, but they can be added to work accordingly in the first release of openkm3.\n",
"\n",
"See http://pi1154.physik.uni-erlangen.de:82/data/collections/1/view for what is already provided"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# read the ANTARES data tracks into a dataframe - PROPOSAL!\n",
"ks = openkm3.KM3Store()\n",
"\n",
"ant07_17_collection = ks.get(\"ant2017\")\n",
"\n",
"ant07_17_collection.list()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"should return a table e.g. \n",
"\n",
"| identifier | name | description | type |\n",
"| --------| --------- | ------- | ------- |\n",
"| \"events\" | ANTARES 2007-2017 neutrino events | Event list of neutrino candidates | km3.data.l4.vo.scs |\n",
"| \"acceptance\" | ANTARES 2007-2017 detector acceptance | acceptance table (for different spectral indices and sin(decl) | km3.data.l4.table.d2 |\n",
"| \"background_distribution\" | ANTARES 2007-2017 background distribution | interpolation of expected background events for given declination and region of interest | km3.data.l4.function.d6 |"
The gammapy package is used to compute Feldman-Cousins statistics
## Prerequisits
`pip install gammapy, matplotlib, openkm3`
%% Cell type:code id: tags:
``` python
importpandasaspd
importnumpyasnp
importmatplotlib.pyplotasplt
importmath
fromscipyimportstats
importgammapy.statsasgstats
importopenkm3
```
%% Cell type:markdown id: tags:
## Helper functions
### Background evaluation
The function *get_number_background(decl, roi)* permits to evaluate the expected background for a give sky declination and a given region of interest ('half width' of the declination band around the source declination). The right ascension is not used here, as the background distribution depends only on the visibility of the given sky position from the location of ANTARES, and therefore on the declination only.
This function uses a polynomial interpolation to smoothen the background distibution otherwise obtained from the data. This function returns the number of expected background events that in the following will be indicated as nbkg.
%% Cell type:code id: tags:
``` python
defget_number_background(interpolation,decl,roi):
dec_low=decl-roi;
dec_high=decl+roi;
# estimation of background rate over the full declination band
The functions *get_upper_limit(nbkg, nobs)* and *get_lower_limit(nbkg, nobs)* compute the upper and lower limits in presence of nbkg backgroud events, if the measurements is nobs observed events inside the region of interest.
%% Cell type:code id: tags:
``` python
# compute upper limit according to Feldamn-Cousins
defget_upper_limit(nbkg,nobs):
ifnobs<=int(nbkg):
ul=Sens(nbkg)#if underfluctuation, UL set as the sensitivity
The function *get_sensitivity(nbkg)* is used to calculate the sensitivity, defined as average upper limit, that is reached in presence of nbkg events in the declination band of the source.
# count events whose angle from the source falls inside the RoI
ifalpha<RoI:
obs=obs+1
returnobs
```
%% Cell type:markdown id: tags:
### Math helpers
The two functions *make_vector(Decl, RA)* and *calc_angle(v1, v2)* are auxilary functions used later.
*make_vector* fills a 3-dimensional vector using the sky coordinates of an event (Decl, RA), and assigning length 1.
*calc_angle* computes the angular distance (in units degrees) between two points in the sky defined as the angle between their corresponding vectors v1 and v2.
Here, the flux upper limits for the observation of ANTARES events from the 2007-2017 sample for a given direction in the sky is calculated.
### Setting the point source coordinates
The limits are calculated providing
* the source coordinates **right ascension** and **declination**, ra_source and dec_source, and
* the two additional source parameters **gamma** (spectral index) and
***RoI** (region of interest, representing the source extension).
%% Cell type:code id: tags:
``` python
ra_source=0.1# between 0 and 360
dec_source=-29# Declination must be between -80 and 50 degree
gamma=2.0# spectral index between 1.5 and 3.0
RoI=2.6# Region of interest > 0.1 and < 1000 (in degree)
```
%% Cell type:markdown id: tags:
### Getting ANTARES data and acceptance
This is a proposal of how retrieving the relevant ANTARES data should work.
The KM3Store should offer
* the event list
* the acceptance function
* the background estimate from interpolation
The following functions are now dummies, but they can be added to work accordingly in the first release of openkm3.
See http://pi1154.physik.uni-erlangen.de:82/data/collections/1/view for what is already provided
%% Cell type:code id: tags:
``` python
# read the ANTARES data tracks into a dataframe - PROPOSAL!
ks=openkm3.KM3Store()
ant07_17_collection=ks.get("ant2017")
ant07_17_collection.list()
```
%% Cell type:markdown id: tags:
should return a table e.g.
| identifier | name | description | type |
| --------| --------- | ------- | ------- |
| "events" | ANTARES 2007-2017 neutrino events | Event list of neutrino candidates | km3.data.l4.vo.scs |
| "acceptance" | ANTARES 2007-2017 detector acceptance | acceptance table (for different spectral indices and sin(decl) | km3.data.l4.table.d2 |
| "background_distribution" | ANTARES 2007-2017 background distribution | interpolation of expected background events for given declination and region of interest | km3.data.l4.function.d6 |
@@ -40,5 +40,5 @@ The FAIR principles provide a solid set of requirements for the development of a
#### Reusable data
***[Licensing standards](https://open-data.pages.km3net.de/licensing/)** for data, software and supplementary material have been introduced.
* Basic **provenance** information is provided with the data, which serves as a development start point to propagate provenance management through the complex the full data processing workflow in the future.
* Basic **provenance** information is provided with the data, which serves as a development start point to propagate provenance management through the complex data processing workflow in the future.
One of the primary goals for ANTARES is the identification of neutrino sources, whose signature would appear in ANTARES data (in form of neutrino arrival directions) as clusters of events at given sky coordinates.
The significance of a neutrino excess from a given sky position must be assessed over the expected background fluctuations using, for instance the Feldman-Cousins statistics.
This ANTARES use case allows to inspect a sample of neutrino arrival directions in equatorial coordinates (RA, dec), evaluate from it the expected background rate for a user-selected sky position, and finally assess the significance of a cluster defined as 'all arrival directions that fall inside a given radius, selected by the user and indicated here 'region of interest' (RoI).
Here follow the step-by-step description of the code provided at https://git.km3net.de/srgozzini/workflow_ps/-/blob/master/PS.ipynb
- The function *bkg_evaluation(decl, roi)* permits to evaluate the expected background for a give sky declination and a given region of interest ('half width' of the declination band around the source declination). The right ascension is not used here, as the background distribution depends only on the visibility of the given sky position from the location of ANTARES, and therefore on the declination only.
This function uses a polynomial interpolation to smoothen the background distibution otherwise obtained from the data. This function returns the number of expected background events that in the following will be indicated as nbkg.
- the functions *UL(nbkg, nobs)* and *LL(nbkg, nobs)* compute the upper and lower limits in presence of nbkg backgroud events, if the measurements is nobs observed events inside the region of interest.
- the function *Sens(nbkg)* is used to calculate the sensitivity, defined as average upper limit, that is reached in presence of nbkg events in the declination band of the source.
- the two functions *make_vector(Decl, RA)* and *calc_angle(v1, v2)* are auxilary functions used later. *make_vector* fills a 3-dimensional vector using the sky coordinates of an event (Decl, RA), and assigning length 1. *calc_angle* computes the angular distance (in units degrees) between two points in the sky defined as the angle between their corresponding vectors v1 and v2.
Here follow the step-by-step description of the code provided as [Jupyter notebook](notebooks/A01_recorded_rate)
- The function *compute_limits(ra_source, dec_source, gamma, RoI)* is the main function of this script and computes the flux upper limits for a given direction in the sky defined from the coordinates ra_source and dec_source, and for the two additional source parameters gamma (spectral index) and RoI (region of interest, representing the source extension).
### Data set
## Using ANTARES data
ANTARES data has already been published to the VO for two data sets by using the services of the German Astrophysical Virtual Observatory (GAVO) which run the DaCHS software. The most recent public data sample of the 2007-2017 point source search is available through the ANTARES website, however, it is thus not findable through VO and does not match the FAIR criteria. Including ANTARES data in the development of future VO data types allowes to increase the chance for a long-term availability of high-quality ANTARES data. On the other hand, the KM3NeT VO server could be registered to the VO and protocols be tested using the ANTARES 2007-2017 sample.
ANTARES data has already been published to the VO for two data sets by using the services of the German Astrophysical Virtual Observatory (GAVO) which run the DaCHS software. The most recent public data sample is available through the ANTARES website, however, it is thus not findable through VO and does not match the FAIR criteria. Including ANTARES data in the development of future VO data types allowes to increase the chance for a long-term availability of high-quality ANTARES data. On the other hand, the KM3NeT VO server could be registered to the VO and protocols be tested using the ANTARES 2007-2017 sample.
The provided data set includes
* The **full event list** of 2007-2017 selected astrophysics neutrino candidates, provided through the VO server,
* Supplementary distributions from simulations provided via the ODC including
* the **detector acceptance** for a given source spectral index and declination
* the interpolated **distribution of background events** for a given declination and region of interest
* the **effective area** for an E^-2 source spectrum in three different zenith bands
