The IceCube Hot-Spot Population Analysis (HPA)

This is a code to run the IceCube Hot-Spot Population Analysis. It is based on a set of scripts previously developed by René Reimann (https://github.com/renereimann/HotSpotPopulationAnalysis).

The Hot-Spot Population Analysis tests for a significant excess of small p-values in a set of p-values (in our case a skymap). These small p-values correspond to signals not significant by themselves, but which could become interesting if looked at as a population of sub-threshold signals.

Take a look to the wiki of the analysis for more details.

Requirements

To run this set of scripts you have two options:

1. Let pipenv manage everything for you:

   • Install pipenv;
   • Clone the repository on your local machine;
   • Run pipenv install: this will automatically create and manage the virtual environment for this project as it is defined in the Pipfile and Pipfile.lock files.

2. Create an environment satisfying the requirements.txt.

3. Use this virtualenv: /home/cbellenghi/.local/share/virtualenvs/HotspotPopulationAnalysis-WJcvmMvY

Note that this set of scripts needs Python>=3.7.

Analysis Workflow

1. Run_hotspot_search.py:

   Performs the search of all hotspots in one or many skymaps and then saves them. It can be run in two different modes:

   • for_parametrization == True: Run the search on many background skymaps and store the output in a list of hotspots per skymap. The output file will be used to generate the background expectation for the analysis afterwards. Default: False.
   • for_parametrization == False: Run the search on many background skymaps and store the output in a list of hotspots per skymap + an array of all hotspots with significance above the minimum -log10(p-value) threshold from all the maps together. The output list will be used to generate the background TS distribution, the output array will be used as the background p-value pool for pseudo-experiment generation.

   The two options are splitted because in principle the background parametrization and the background trials should be generated from independent sets of skymaps.

   The background skyscans for parametrization can be found here: /data/user/cbellenghi/hpa/data/background_skyscans/for_parametrization; the background skyscans for trtials generation can be found here: /data/user/cbellenghi/hpa/data/background_skyscans/for_trials_generation.

2. build_background_expectation.py:

   Constructs the counts expectation as a function of the -log10(p_local) threshold from background skymaps. By default the expectation is built scanning -log10(p_local) threhsholds from 2 to 7 in steps of 0.1.

   A list of background hotspots to build the expectation can be found here: /data/user/cbellenghi/hpa/data/background_hotspots/hs_for_parametrization.pkl.

   The built background expectation spline can be found here: /data/user/cbellenghi/hpa/data/background_hotspots/expectation.pkl.

3. generate_bkg_ts_distribution.py:

   Generates the background TS distribution and parametrize it using a gamma function. Then saves the trials and the parametrization.

   N.B. Needs the expectation generated at step 2. as input.

   A list of background hotspots to generate the background TS distribution can be founf here: /data/user/cbellenghi/hpa/data/background_hotspots/hs_for_bg_ts_distribution.pkl.

   Already computed background trials can be found here: /data/user/cbellenghi/hpa/data/analysis_output/background_hpa_result.npy

   The best fit parameters for the gamma functions can be found here: /data/user/cbellenghi/hpa/data/analysis_output/gamma_fit_result.pkl

4. prepare_signal_injection_pool.py:

   Prepares and saves a pool of p-values from single source trials from which the p-values for signal injection to generate HPA trials will be chosen. The output signal pool is a structured array with the following fields:

   • dec: injection declination in degrees;
   • gamma: spectral index for the injected power law flux model;
   • n_inj: the true number of injected neutrinos;
   • logpVal: the significance in -log10(p_local)
   • mu_per_flux: conversion factor from flux normalisation to number of injected neutrinos.

   N.B. Since the conversion factor from flux normalisation to number of injected neutrinos needs to be computed, one of the inputs for the script is the MC filepath.

   The point source trials generated injecting a spectral index 2.0 can be found here: /data/user/cbellenghi/hpa/data/single_source_signal_injection_trials/joint_trials.

   A built pool of signals for the HPA can be found here: /data/user/cbellenghi/hpa/data/signal_pools/signal_pool_E0_1GeV_mode_region.npy

5. run_trials.py:

   Generates HPA trials. Takes as input the background (/data/user/cbellenghi/hpa/data/background_hotspots/bg_hotspots_pool.npy) and signal (/data/user/cbellenghi/hpa/data/signal_pools/signal_pool_E0_1GeV_mode_region.npy) pools and the background expectation (/data/user/cbellenghi/hpa/data/background_hotspots/expectation.pkl). In order to generate trials for populations of equal strenght at Earth sources, the additional arguments nsrc and src_flux must be given; to generate population of neutrino sources using FIRESONG, the arguments density and fluxnorm must be given.

   N.B. FIRESONG is a python package to simulate populations of extragalactic neutrino sources (DOI: 10.21105/joss.03194) and is used here with the default source density evolution (Madau and Dickinson 2016 Star Formation Rate History - MD2014SFR) and a hard spectral index of 2.0.

   The output is a structured array with the following fields:

   • n_inj: the true total number of injected neutrinos (for the entire population);
   • max_logp: the largest -log10(p_local) in the population;
   • hpa_ts: the -log10 of the minimum poisson p-value (TS of this analysis);
   • log10p_thres: the -log10(p_local) corresponding to the highest TS value;
   • n_observed: number of that maximally deviates from the expectation (i.e. at log10p_thres);
   • n_expected: number of expected hotspots at log10p_thres.

   Some trials for equal strength sources can be found here: /data/user/cbellenghi/hpa/data/analysis_output/signal_trials/test_diffuse_contrib/nu_per_1deg.

6. estimate_sensitivity.py:

   Estimates the sensitivity, 3 and 5 sigma discovery potential of the analysis.

   Takes one or more trial files in input, together with the background TS distribution parametrization generated at step 3 (/data/user/cbellenghi/hpa/data/analysis_output/gamma_fit_result.pkl). The parametrization is needed to compute the critical TS values corresponding to interesting quantiles of the distribution.

   The output is a list of estimate_mu_per_ts_quantile objects (as defined in the library analysis_utils.py) having the following attributes:

   • label: one among 'sens', '3sigma disc pot', '5sigma disc pot';
   • alpha and beta: values to define the confidence level of the analysis (i.e. p-value smaller than alpha in a fraction beta of the cases);
   • ts_thres: critical TS value corresponding to the confidence level defined by alpha and beta;
   • success: boolean to assert whether the minimization performed to get the mean number of neutrinos giving the required confidence level succeeded;
   • tot_n_inj: the mean number of neutrinos got from the fit;
   • flux_per_mu: the factor from number of injected neutrinos to flux normalisation;
   • flux_units: 1/[GeV cm^2 s];
   • flux_pivot_point: 1 GeV.

   N.B. Since the conversion factor from number of injected neutrinos to flux normalisation needs to be computed, also for this script one of the inputs is the MC filepath.

Additional material

• The notebooks/ folder contains a set of ipython notebook to reprocude useful plots for this analysis.
• All the output files produced by the scripts can be found here: /data/user/cbellenghi/hpa/data. They can be used as input to run the HPA code.