NPTFit Documentation

NPTFit is a specialized Python/Cython package that implements Non-Poissonian Template Fitting (NPTF), originally developed for characterizing populations of unresolved point sources. The main features of the package are

• Fast evaluation of likelihoods for NPTF analyses
• Easy-to-use interface for performing non-Poissonian (as well as standard Poissonian) template fits using MultiNest or other inference tools
• Ability to include an arbitrary number of point source templates, with an arbitrary number of degrees of freedom in the modeled flux distribution
• Modules for analyzing and plotting the results of an NPTF

The most up-to-date version of the code can be found at https://github.com/bsafdi/NPTFit.

Installation

Out of the box, NPTFit relies on MultiNest for Bayesian inference, which must be installed and linked prior to use. See here for a helpful set of installation instructions.

NPTFit supports both Python 2 and 3, specifically 2.7 and 3.5. It may work with earlier 3.* versions, although this has not been tested.

Make sure Cython is installed (e.g. pip install Cython). The easiest way to install NPTFit along with it’s dependent Python packages is using pip:

$ pip install NPTFit

or using the setup script:

$ python setup.py install

which also builds the Cython modules. To just compile the Cython modules locally:

$ make build

NPTFit uses either GSL (faster) or mpmath for back-end calculations. By default, the GSL version is compiled during setup if this is available and linked, in which case the setup script concludes with "GSL compilation successful!". Otherwise, "mpmath compilation successful!" is shown.

The code is parallelizable through MPI (e.g. OpenMPI), which can considerably speed up computationally intensive scans. This requires the MPI4Py Python package for use with MultiNest, which can be installed, for example, with pip:

$ pip install mpi4py

Contents

Initializing a scan

The class nptfit.NPTF is used to instantiate, configure and run a template fit. An instance of nptfit.NPTF can be created as follows:

>>> from NPTFit import nptfit
>>> nptf = nptfit.NPTF(tag, work_dir)

Description of arguments (all optional):

Argument	Default	Purpose
tag	"Untagged"	Label associated with the scan
work_dir	None	Location where all the output from the scan is stored

Note

The current directory for work_dir is used if not provided.

Loading data and exposure

All maps are passed as 1-d numpy arrays. The data and exposure maps can be loaded into an instance of nptfit.NPTF (see Initializing a scan) as follows:

>>> nptf.load_data(data, exposure)

Adding templates

Spatial templates can be added as follows:

>>> nptf.add_template(template, key, units)

where the first argument is the template numpy array and the second argument is the template key, used to identify the template in later calls.

The argument units specifies the template units (counts or flux) or type (to be used in either Poissonian or PS models). The following values are allowed:

• 'counts': template in counts/pixel, to be used in a Poissonian model. Exposure and PSF corrected.
• 'flux': template in counts/cm^2/s/pixel, to be used in a Poissonian model. Not exposure corrected.
• 'PS': template for the underlying PS distribution, to be used in a non-Poissonian model. This shouldn’t account for exposure effects.

For example, a template specifying an underlying isotropic PS distribution for a non-Poissonian model would be added with the keyword 'PS' as a truly isotropic array, e.g. [1,1,1,...,1].

Warning

The data, exposure and template maps must be of the same length.

Note

While the data, exposure and template maps need not be in HEALPix format, creation of masks is only supported for HEALPix formatted arrays. Conversion between Cartesian and HEALPix maps can be performed using grid2healpix

Adding masks

Masks can be used to restrict an analysis within a given region of interest (ROI). Masks are boolean arrays where pixels labelled as true are masked and those labelled false are unmasked. Masks can be loaded into an instance of nptfit.NPTF (see Initializing a scan) as:

>>> nptf.load_mask(mask)

Warning

The mask must be the same length as the data.

Creating masks

For the case where the data is a HEALPix array, masks can be created by the create_mask module. For example:

>>> from NPTFit import create_mask as cm
>>> mask = cm.make_mask_total(nside, mask_ring = True, outer = 2, ring_b = -45)

where nside is the HEALPix resolution parameter. The following mask options are available. When set to True, the associated arguments determine the parameters of the mask. If multiple options are set to True then the total mask is the combination of each option.

Argument	Default	Purpose
band_mask	False	If True, masks within |b| < band_mask_ range
l_mask	False	If True, masks longitude outside l_deg_min < l < l_deg_max
b_mask	False	If True, masks latitude outside b_deg_min < b < b_deg_max
mask_ring	False	If True, masks outside inner < r < outer, of a ring centred at (ring_b, ring_l )
custom_mask	None	Optional user-provided mask

Note

By convention, True pixels are masked and False unmasked.

Note

Creation of masks by create_mask is only supported for HEALPix formatted maps. For other inputs, masks otherwise made can be used as usual with nptf.load_mask().

Tip

See the Example 2: Creating Masks for an exposition of the options described here.

Adding models

From the list of templates that nptfit.NPTF knows about (see Loading data and exposure), we can define an arbitrary number of Poissonian (smooth/diffuse) and non-Poissonian (point source) models.

Poissonian models

Poissonian models only have one parameter associated with them: the template normalisation. Poissonian models corresponding to an added template can be loaded as:

nptf.add_poiss_model(template_name='iso',
                     model_tag='$A_{iso}$',
                     prior_range=[-2,2],
                     log_prior=False,
                     fixed=False,
                     fixed_norm=1.0)

Arguments for Poissonian model:

Argument	Default	Purpose
template_name	-	Key corresponding to previously loaded template
model_tag	-	LaTeX-ready string, for plots
prior_range	[min,max]	Prior range to scan over
log_prior	False	Whether to scan template in log-space
fixed	False	Whether the model is fixed
fixed_norm	1.0	Template normalization if the model is fixed

Note

When log_prior=True, the associated values of prior_range and fixed_norm must also be the logs of the values being used.

Warning

When log_prior=True, the log used for the parameters is base 10.

Non-Poissonian models

The flux distribution of non-Poissonian (point source) models is modeled as a multiply broken power law with a specified number of breaks \(l\), the best-fit parameters of which can then be inferred. This source count distribution, which gives the differential number of sources per unit of flux, takes the form

\[\begin{split}\frac{dN}{dF} = A \left\{ \begin{array}{lc} \left( \frac{F}{F_{b,1}} \right)^{-n_1}, & F \geq F_{b,1} \\ \left(\frac{F}{F_{b,1}}\right)^{-n_2}, & F_{b,1} > F \geq F_{b,2} \\ \left( \frac{F_{b,2}}{F_{b,1}} \right)^{-n_2} \left(\frac{F}{F_{b,2}}\right)^{-n_3}, & F_{b,2} > F \geq F_{b,3} \\ \left( \frac{F_{b,2}}{F_{b,1}} \right)^{-n_2} \left( \frac{F_{b,3}}{F_{b,2}} \right)^{-n_3} \left(\frac{F}{F_{b,3}}\right)^{-n_4}, & F_{b,3} > F \geq F_{b,4} \\ \\ \ldots & \ldots \\ \\ \left[ \prod_{i=1}^{\ell-1} \left( \frac{F_{b,i+1}}{F_{b,i}} \right)^{-n_{i+1}} \right] \left( \frac{F}{F_{b,\ell}} \right)^{-n_{\ell+1}} & F_{b,\ell} > F \end{array} \right.\end{split}\]

It is sometimes convenient to specify the breaks in terms of counts instead of flux. However, if the exposure map is non-uniform over the ROI, then the notion of counts in pixel dependent. While the NPTFit code properly accounts for the pixel-dependent exposure correction, we also allow the user the specify the breaks \(F_{b,i}\) in terms of an effective number of counts \(S_{b,i} \equiv F_{b,i} \cdot \text{mean}_\text{ROI}(E_p)\), where \(\text{mean}_\text{ROI}(E_p)\) is the mean of the exposure map \(E_p\) over the ROI.

Non-Poissonian models corresponding to an added template can be loaded as:

nptf.add_non_poiss_model(template_name='iso',
                         model_tag=['$A_{ps}$','$n_1$','$n_2$','$S_b^{(1)}$'],
                         prior_range=[[-6,6],[2.05,30],[-2,1.95],[0.05,30.0]],
                         log_prior=[True,False,False,False],
                         dnds_model='specify_breaks',
                         fixed_params=None,units='counts')

Arguments for non-Poissonian model:

Argument	Default	Purpose
template_name	-	Key corresponding to loaded template
model_tag	-	LaTeX-ready string of nb-broken power-law parameters [A, n_1, … , n_{nb+1}, Sb_1, … , Sb_nb] with Sb_1/n_1 the highest break/slope
prior_range	[]	Prior range to scan over, given as a list of [min,max] each model parameter
log_prior	False	Whether to scan each parameter in log-space
dnds_model	specify_breaks	Whether to use absolute or relative breaks (see \(dN/dS\) model specifications below)
fixed_params	None	Which parameters to keep fixed (see Fixed parameter specifications below)
units	counts	Whether the breaks in dN/dF are specified in terms of \(F_{b}\) or \(S_{b}\), which is defnied above. (see units specifications below)

Note

The number of breaks in the non-Poissonian model is inferred from the length of the model_tag array.

Warning

Non-Poissonian (or PS) models must use a template loaded with units='PS', while non-Poissonian models should use units='counts' or units='flux'.

\(dN/dF\) model specifications

where \(n_1\) is the highest index and \(S_b^{(1)}\) the highest break.

The following options are allowed for dnds_model:

• specify_breaks: all breaks are specified in absolute counts, \(\left[ A, n_1, \ldots, n_{\ell+1}, S_b^{(1)}, \ldots, S_b^{(\ell)} \right]\)
• specify_relative_breaks: the highest break is specified in counts, with each subsequent lower break specified relative to the subsequent higher break. \(\left[ A, n_1, \ldots, n_{\ell+1}, S_b^{(1)}, \lambda^{(2)}, \ldots, \lambda^{(\ell - 1)}, \lambda^{(\ell)} \right]\) where \(\lambda^{(i)} = S_b^{(i)}/S_b^{(i-1)}\).

Fixed parameter specifications

Fixed parameters should be passed as an array with syntax

fixed_params = [[param_index_1,fixed_value_1],[param_index_2,fixed_value_2]]

where parameter indexing starts from 0.

Units specifications

The following options are allowed for units:

• counts: The code assumes the user has specified the break in counts (\(S_b\)) and will infer the breaks in flux (\(F_b\)) by dividing the mean exposure: \(F_{b,i} \equiv S_{b,i} / \text{mean}_\text{ROI}(E_p)\).
• flux: The code assumes the breaks are already specified in terms of flux

Tip

See Example5_Running_nonPoissonian_Scans and Example7_Galactic_Center_nonPoissonian for examples of using these options in an analysis.

Running the scan

The scan can be configured as follows:

>>> nptf.configure_for_scan(f_ary, df_rho_div_f_ary, nexp)

Note

For no PSF correction, simply leave f_ary and df_rho_div_f_ary to their default values.

Note

If no non-Poissonian models are added when configure_for_scan is called, the likelihood will default to a pure Poissonian scan.

From here the scan can be run using:

>>> nptf.perform_scan(run_tag='example_run', nlive, pymultinest_options=None)

Scan options:

Argument	Default	Purpose
run_tag	None	An optional custom label for the folder within work_dir/chains/tag where the MultiNest output is stored
nlive	100	Number of live sampling points in MultiNest
pymultinest_options	None	A custom set of MCMC options for pymultinest (see below)

Setting custom MultiNest options

A custom set of MCMC options for pymultinest can be provided to nptfit.NPTF.perform_scan() with the argument pymultinest_options. This should be a dictionary of the form {'option': value ...}. The default options are:

pymultinest_options = {importance_nested_sampling: False,
                       resume: False, verbose: True,
                       make_or_load_psf_corr()sampling_efficiency: model,
                       init_MPI: False, evidence_tolerance: 0.5,
                       const_efficiency_mode: False}

Parallel implementation

MultiNest uses MPI for parallel sampling, which can significantly speed up the NPTF. This can be seamlessly implemented as with NPTFit. For example, to run on 12 CPU cores:

>>> mpirun -np 12 nptf_run.py

Analyzing results of a scan

Tip

See the Example 9: Analyzing the Results of an NPTFit Run for an interactive version of the analysis options described here.

While the chain samples of a non-Poissonian fit performed using MultiNest can be readily accessed, a basic analysis module dnds_analysis.py is provided which includes the basic functions described below.

Initializing the analysis module

Having performed a scan using an instance of nptfit.NPTF, the first thing to do is load the scan parameters in. This is done with

>>> nptf.load_scan()

where nptf is an instance of nptfit.NPTF.

Note

The analysis can be performed on an already existing scan. To do this an instance of nptfit.NPTF should be created with the same template and model configuration as used to perform the scan. Then instead of running the scan, simply load it and proceed with the analysis as below.

An instance of the analysis module can then be created as follows:

>>> an = dnds_analysis.Analysis(nptf, mask=None, pixarea=0.)

where mask specifies the ROI used for the analysis if this is different from that used for the scan, and pixarea is the area of a pixel in sr if using non-HEALPix maps.

Making triangle plots

Triangle/corner plots can be used to visualize multidimensional samples using a scatterplot matrix. A triangle plot with the default options using the corner package can be made using:

>>> an.make_triangle()

To use your own custom plotting options, use corner as follows

>>> corner.corner(an.nptf.samples, labels=an.nptf.params, range=[1 for i in range(an.nptf.n_params)])

with additional arguments as specified in http://corner.readthedocs.io/en/latest/.

Getting template intensities

Template intensities [counts/cm²/s/sr] can be calculated with

>>> an.return_intensity_arrays_poiss(comp)
>>> an.return_intensity_arrays_non_poiss(comp)

for Poissonian and non-Poissonian templates with the key comp respectively. This returns a list of intensities corresponding to the posterior parameters for the given template.

The non-Poissonian templates (NPT) intensity is calculated by integrating up \(\int_{S_{min}}^{S_{max}} dS~S~dN/dS\). This is approximated as a sum between \(S_{min}\) and \(S_{max}\). The options associated with the non-Poissonian template intensity are:

Argument	Default	Purpose
comp	-	The NPT key
smin	0.01	Minimum counts to sum up from
smax	10000	Maximum counts to sum up to
nsteps	10000	Number of bins in s while summing up

Getting source count distributions

The posterior arrays for the source count distributions \(dN/dF\) [counts^-1 cm² s deg^-2] associated with a given template comp at a given flux (in counts/cm²/s) can be obtained using

>>> an.return_dndf_arrays(comp,flux)

The source count distribution can be plotted with

>>> an.plot_source_count_median(comp, smin, smax, nsteps, spow, **kwargs)
>>> an.plot_source_count_band(comp, smin, smax, nsteps, spow, qs, **kwargs)

The options being the same as for obtaining the NPT intensity above. Additionally, spow is the power \(n\) in \(F^ndN/dF\) to return while plotting, and qs is an array of quantiles for which to return the dN/dF band.

Plotting intensity fractions

Intensity fractions (fraction of template intensity to total intensity) for Poissonian and non-Poissonian templates respectively can be plotting using

>>> an.plot_intensity_fraction_poiss(comp, bins, **kwargs)
>>> an.plot_intensity_fraction_non_poiss(comp, bins, **kwargs)

where comp is the template key, bins is the number of bins between 0 and 100 and **kwargs specify plotting options.

Accessing posteriors

While the posteriors can be accessed with nptf.samples (or an.nptf.samples) as above, the following functions provide a useful interfact to access individual parameters:

>>> an.return_poiss_parameter_posteriors(comp)
>>> an.return_non_poiss_parameter_posteriors(comp)

where comp is the (non-)Poissonian template key.

For Poissonian models, this returns a list of posterior normalizaion parameters for that model. For non-Poissonian models, this returns three arrays:

>>> A_non_poiss_post, n_non_poiss_post, Sb_non_poiss_post = an.return_non_poiss_parameter_posteriors(comp)

where

• A_non_poiss_post is an array of non-Poissonian normalization parameter posteriors
• n_non_poiss_post is a 2-D array, each sub-array containing posteriors for a given slope parameter, starting from the highest to the lowest
• Sb_non_poiss_post is a 2-D array, each sub-array containing posteriors for a given break parameter, starting from the highest to the lowest

Getting Bayesian log-evidences

The Bayesian log-evidence and associated error can be accessed as follows:

>>> l_be, l_be_err = an.get_log_evidence()

Example notebooks

Example 1: Overview of the Fermi Data

This example provides an overview of the Fermi-LAT gamma-ray data we have included with the code, which provide the testing ground for the examples in later notebooks.

We emphasize that NPTFit is applicable much more generally. The code package can be implemented in any scenario in which the data is measured in discrete counts. We use the example of Fermi data here, as it represented some of the earliest applications of this code and also it is fully publicly available, see http://fermi.gsfc.nasa.gov/ssc/data/access/

The specifics of the dataset we use are given below.

Parameter	Value
Energy Range	2-20 GeV
Time Period	Aug 4, 2008 to July 7, 2016 (413 weeks)
Event Class	UltracleanVeto (1024) - highest cosmic ray rejection
Event Type	PSF3 (32) - top quartile graded by angular reconstruction
Quality Cuts	DATA_QUAL==1 && LAT_CONFIG==1
Max Zenith Angle	90 degrees

# Import relevant modules

%matplotlib inline
%load_ext autoreload
%autoreload 2

import os

import numpy as np
import healpy as hp

Downloading the data

The Fermi data used here can be downloaded from https://dspace.mit.edu/handle/1721.1/105492 along with a README detailing the data. The following commands automatically download and set up the data on a UNIX system:

# Assumes wget is available! Otherwise, use curl or download manually
# from https://dspace.mit.edu/handle/1721.1/105492
os.system("wget https://dspace.mit.edu/bitstream/handle/1721.1/105492/fermi_data.tar.gz?sequence=5");
os.system("tar -xvf fermi_data.tar.gz?sequence=5");
os.system("rm -r fermi_data.tar.gz*");

The Data - Map of Gamma-ray Counts

The data we have provided is in a single energy bin, so only spatial information of the photons remains. We bin this into an nside=128 healpix map. This map has units of [photon counts/pixel]. The data is shown below on two scales, the second making it clear there are pixels with large numbers of counts - there are very bright gamma-ray sources in the data.

counts = np.load('fermi_data/fermidata_counts.npy')

hp.mollview(counts,title='Data counts',max=100)

hp.mollview(counts,title='Data counts (no max)')

To see the detailed structure in the map, we also mock up a logarithmic version of the data.

nonzero = np.where(counts != 0)[0]
zero = np.where(counts == 0)[0]
counts_log = np.zeros(len(counts))
counts_log[nonzero] = np.log10(counts[nonzero])
counts_log[zero] = -1
hp.mollview(counts_log,title='Data counts - Log')

Exposure Map - Map of where Fermi-LAT looks

Fermi-LAT is a space based telescope that takes data from the entire sky. Nevertheless it does not look at every part of the sky equally, and the exposure map keeps track of this.

The natural units for diffuse astrophysical sources is intensity [counts/cm:math:^2/s/sr], whereas for point sources we use flux [counts/cm:math:^2/s]. Neither of these knows Fermi was pointing. But Fermi measures counts or counts/pixel, and this is also the space in which we perform our statistical analysis. The mapping from [counts/cm:math:^2/s(/sr)] to [counts(/pixel)] is performed by the exposure map, which has units of [cm:math:^2 s].

exposure = np.load('fermi_data/fermidata_exposure.npy')

hp.mollview(exposure,title='Fermi Exposure [cm$^2$ s]')

When performing an NPTF, technically the non-Poissonian templates should be separately exposure corrected in every pixel. Doing this exactly is extremely computationally demanding, and so instead we approximate this by breaking the exposure map up into regions of approximately similar exposure values.

For the Fermi instrument a small number of exposure values (set by nexp at the point of configuring a scan) is often sufficient, as the exposure is quite uniform over the sky. For datasets with less uniform exposure, however, larger values of nexp are recommended. We show a run performed with nexp != 1 in Example 4 and 10.

Below we show how the sky is divided into different exposure regions - try changing nexp.

NB: In the actual analysis the exposure region division is done within the specified ROI, not the entire sky

# Number of exposure regions - change this to see what the regions look like when dividing the full sky
nexp = 2

# Divide the exposure map into a series of masks
exp_array = np.array([[i, exposure[i]] for i in range(len(exposure))])
array_sorted = exp_array[np.argsort(exp_array[:, 1])]
array_split = np.array_split(array_sorted, nexp)
expreg_array = [[int(array_split[i][j][0]) for j in range(len(array_split[i]))] for i in range(len(array_split))]
temp_expreg_mask = []
for i in range(nexp):
    temp_mask = np.logical_not(np.zeros(len(exposure)))
    for j in range(len(expreg_array[i])):
        temp_mask[expreg_array[i][j]] = False
    temp_expreg_mask.append(temp_mask)
expreg_mask = temp_expreg_mask

for ne in range(nexp):
    hp.mollview(expreg_mask[ne],title='Fermi Exposure Region '+str(ne+1),min=0,max=1)

Point Source Catalog Mask

We also include a mask of all point sources in the 3FGL, as well as large extended objects such as the Large Magellanic Cloud. All point sources are masked at 95% containment as appropriate for this nside. The map below is a mask, so just a boolean array

Note that with the NPTF it is not always desirable to mask point sources - we can often simply use a non-Poissonian template to model them.

pscmask = np.load('fermi_data/fermidata_pscmask.npy')

hp.mollview(pscmask, title='Point Source Catalog Mask')

Templates - Spatial Models for the Fermi Data

Next we show the different spatial templates that we will use to model the above Fermi data. These also represent examples for the types of models we can use as a basis for either Poissonian or non-Poissonian templates.

Note that templates given to the NPTF must be exposure corrected. That is they should be models of counts, not flux. Furthermore they should also be smoothed to account for the PSF if necessary.

In addition to exposure correcting the maps, for each template below we have also adjusted it so that it has mean 1 within an ROI defined by \(|b|>2^{\circ}\) and \(r<30^{\circ}\).

Diffuse Emission

Firstly we show a model for the diffuse emission, which arises mainly from three sources: 1. protons hitting the gas, giving rise to pions which then decay to photons (\(pp \to \pi^0 \to \gamma \gamma\)); 2. inverse compton scattering from electrons upscattering starlight or the CMB; and 3. bremsstrahlung off of ambient gas.

This model accounts for the majority of the Fermi data. We use the Fermi p6v11 model for the purpose (as it does not also include a template for the Fermi bubbles which we model separately).

Below we show a log and linear version of the map, as we did for the data.

dif = np.load('fermi_data/template_dif.npy')

hp.mollview(dif,title='Diffuse Model Template (p6v11)',max=10)

hp.mollview(np.log10(dif),title='Diffuse Model Template (p6v11) - Log')

Isotropic Emission

There is also an approximately isotropic contribution to the data from extragalactic emission and also cosmic ray contamination. Note that this map makes the fact the template has been exposure corrected manifest.

iso = np.load('fermi_data/template_iso.npy')

hp.mollview(iso,title='Isotropic Emission Template')

Fermi Bubbles

We also provide a separate model for emission from the Fermi bubbles. Emission from the bubbles is taken to be uniform in intensity, which becomes non-uniform in counts after exposure correction.

bub = np.load('fermi_data/template_bub.npy')

hp.mollview(bub,title='Fermi Bubbles Template')

Point Source Catalog Model

As seen in the initial data, the gamma-ray sky includes some extremely bright point sources. As such even if a mask is used to largely cover these, it is still often a good idea to model the point sources as well. Below we show the template for these point sources.

A linear plot of this map shows that these point sources are quite localised on the sky, but in the log plot below we can clearly see their spread due to the Fermi PSF.

psc = np.load('fermi_data/template_psc.npy')

hp.mollview(psc,title='Point Source Catalog Template',max=50)

hp.mollview(np.log10(psc),title='Point Source Catalog Template - Log')

Model for the Galactic Center Excess

Finally we include a model to describe the Galactic Center Excess (GCE). Regardless of the origin of this excess, it has been shown to be spatially distributed as approximately a squared generalized Navarro–Frenk–White (NFW) profile integrated along the line of sight. The generalized NFW for the Milky Way has the form:

\[\rho(r) \propto \frac{(r/r_s)^{-\gamma}}{(1+r/r_s)^{3-\gamma}}\,,\]

where \(r\) is the distance from the Galactic center. We take \(r_s = 20\) kpc and \(\gamma = 1.0\). The flux GCE template is then formed as:

\[J(\psi) = \int_{\rm los} \rho^2(r) ds\,,\]

where \(s\) parameterizes the line of sight distance, which is integrated over, and \(\psi\) is the angle away from the Galactic center.

gce = np.load('fermi_data/template_gce.npy')

hp.mollview(gce,title='Galactic Center Excess Template',max=50)

hp.mollview(np.log10(gce),title='Galactic Center Excess Template - Log')

Model for Disk Correlated Point Sources

When studying the point source origin of the GCE - done in Example 8 - we will also include a model for point sources correlated with the disk of the Milky Way. For this purpose we use the following thick disk double exponential model for the point source population:

\[n(z,R) \propto \exp \left( - R / 5~\textrm{kpc} \right) \exp \left( - |z| /1~\textrm{kpc} \right)\,,\]

where \(R\) and \(z\) are cylindrical polar coordinates measured from the Galactic Center.

disk = np.load('fermi_data/template_dsk.npy')

hp.mollview(disk,title='Thick Disk')

hp.mollview(np.log10(disk),title='Thick Disk - Log')

Example 2: Creating Masks

In this example we show how to create masks using create_mask.py.

Often it is convenient to consider only a reduced Region of Interest (ROI) when analyzing the data. In order to do this we need to create a mask. The masks are boolean arrays where pixels labelled as True are masked and those labelled False are unmasked. In this notebook we give examples of how to create various masks.

The masks are created by create_mask.py and can be passed to an instance of nptfit via the function load_mask for a run, or an instance of dnds_analysis via load_mask_analysis for an analysis. If no mask is specified the code defaults to the full sky as the ROI.

NB: Before you can call functions from NPTFit, you must have it installed. Instructions to do so can be found here:

http://nptfit.readthedocs.io/

# Import relevant modules

%matplotlib inline
%load_ext autoreload
%autoreload 2

import numpy as np
import healpy as hp

from NPTFit import create_mask as cm # Module for creating masks

Example 1: Mask Nothing

If no options are specified, create mask returns an empty mask. In the plot here and for those below, purple represents unmasked, yellow masked.

example1 = cm.make_mask_total()
hp.mollview(example1, title='', cbar=False, min=0,max=1)

Example 2: Band Mask

Here we show an example of how to mask a region either side of the plane - specifically we mask 30 degrees either side

example2 = cm.make_mask_total(band_mask = True, band_mask_range = 30)
hp.mollview(example2, title='', cbar = False, min=0, max=1)

Example 3: Mask outside a band in b and l

This example shows several methods of masking outside specified regions in galactic longitude (l) and latitude (b). The third example shows how when two different masks are specified, the mask returned is the combination of both.

example3a = cm.make_mask_total(l_mask = False, l_deg_min = -30, l_deg_max = 30,
                               b_mask = True, b_deg_min = -30, b_deg_max = 30)
hp.mollview(example3a,title='',cbar=False,min=0,max=1)

example3b = cm.make_mask_total(l_mask = True, l_deg_min = -30, l_deg_max = 30,
                               b_mask = False, b_deg_min = -30, b_deg_max = 30)
hp.mollview(example3b,title='',cbar=False,min=0,max=1)

example3c = cm.make_mask_total(l_mask = True, l_deg_min = -30, l_deg_max = 30,
                              b_mask = True, b_deg_min = -30, b_deg_max = 30)
hp.mollview(example3c,title='',cbar=False,min=0,max=1)

Example 4: Ring and Annulus Mask

Next we show examples of masking outside a ring or annulus. The final example demonstrates that the ring need not be at the galactic center.

example4a = cm.make_mask_total(mask_ring = True, inner = 0, outer = 30, ring_b = 0, ring_l = 0)
hp.mollview(example4a,title='',cbar=False,min=0,max=1)

example4b = cm.make_mask_total(mask_ring = True, inner = 30, outer = 180, ring_b = 0, ring_l = 0)
hp.mollview(example4b,title='',cbar=False,min=0,max=1)

example4c = cm.make_mask_total(mask_ring = True, inner = 30, outer = 90, ring_b = 0, ring_l = 0)
hp.mollview(example4c,title='',cbar=False,min=0,max=1)

example4d = cm.make_mask_total(mask_ring = True, inner = 0, outer = 30, ring_b = 45, ring_l = 45)
hp.mollview(example4d,title='',cbar=False,min=0,max=1)

Example 5: Custom Mask

In addition to the options above, we can also add in custom masks. In this example we highlight this by adding a random mask.

random_custom_mask = np.random.choice(np.array([True, False]), hp.nside2npix(128))
example5 = cm.make_mask_total(custom_mask = random_custom_mask)
hp.mollview(example5,title='',cbar=False,min=0,max=1)

Example 6: Full Analysis Mask including Custom Point Source Catalog Mask

Finally we show an example of a full analysis mask that we will use for an analysis of the Galactic Center Excess in Example 3 and 8. Here we mask the plane with a band mask, mask outside a ring and also include a custom point source mask. The details of the point source mask are given in Example 1.

NB: before the point source mask can be loaded, the Fermi Data needs to be downloaded. See details in Example 1.

pscmask=np.array(np.load('fermi_data/fermidata_pscmask.npy'), dtype=bool)
example6 = cm.make_mask_total(band_mask = True, band_mask_range = 2,
                              mask_ring = True, inner = 0, outer = 30,
                              custom_mask = pscmask)
hp.mollview(example6,title='',cbar=False,min=0,max=1)

Example 3: Running Poissonian Scans with MultiNest

In this example we demonstrate how to run a scan using only templates that follow Poisson statistics. Nevertheless many aspects of how the code works in general, such as initialization, loading data, masks and templates, and running the code with MultiNest carry over to the non-Poissonian case.

In detail we will perform an analysis of the inner galaxy involving all five background templates discussed in Example 1. We will show that the fit prefers a non-zero value for the GCE template.

NB: This example makes use of the Fermi Data, which needs to already be installed. See Example 1 for details.

# Import relevant modules

%matplotlib inline
%load_ext autoreload
%autoreload 2

import numpy as np
import corner
import matplotlib.pyplot as plt

from NPTFit import nptfit # module for performing scan
from NPTFit import create_mask as cm # module for creating the mask
from NPTFit import dnds_analysis # module for analysing the output

Step 1: Setting up an instance of NPTFit

To begin with we need to create an instance of NPTF from nptfit.py. We will load it with the tag set to “Poissonian_Example”, which is the name attached to the folder within the chains directory where the output will be stored. Note for long runs the chains output can become large, so periodically deleting runs you are no longer using is recommended.

n = nptfit.NPTF(tag='Poissonian_Example')

The full list of parameters that can be set with the initialization are as follows (all are optional).

Argument	Defaults	Purpose
tag	“Untagged”	The label of the file where the output of MultiNest will be stored, specifically they are stored at work_dir/chains/tag /.
work_dir	$pwd	The directory where all outputs from the NPTF will be stored. This defaults to the notebook directory, but an alternative can be specified.
psf_dir	work_dir/psf_dir/	Where the psf corrections will be stored (this correction is discussed in the next notebook).

Step 2: Add in Data, a Mask and Background Templates

Next we need to pass the code some data to analyze. For this purpose we use the Fermi Data described in Example 1. The format for load_data is data and then exposure.

NB: we emphasize that although we use the example of HEALPix maps here, the code more generally works on any 1-d arrays, as long as the data, exposure, mask, and templates all have the same length.

Code to embed data on a regular Cartesian grid into a HEALPix map can be found here: https://github.com/nickrodd/grid2healpix.

fermi_data = np.load('fermi_data/fermidata_counts.npy').astype(np.int32)
fermi_exposure = np.load('fermi_data/fermidata_exposure.npy')
n.load_data(fermi_data, fermi_exposure)

In order to study the inner galaxy, we restrict ourselves to a smaller ROI defined by the analysis mask discussed in Example 2. The mask must be the same length as the data and exposure.

pscmask=np.array(np.load('fermi_data/fermidata_pscmask.npy'), dtype=bool)
analysis_mask = cm.make_mask_total(band_mask = True, band_mask_range = 2,
                                   mask_ring = True, inner = 0, outer = 30,
                                   custom_mask = pscmask)
n.load_mask(analysis_mask)

Add in the templates we will want to use as background models. When adding templates, the first entry is the template itself and the second the string by which it is identified. The length for each template must again match the data.

dif = np.load('fermi_data/template_dif.npy')
iso = np.load('fermi_data/template_iso.npy')
bub = np.load('fermi_data/template_bub.npy')
psc = np.load('fermi_data/template_psc.npy')
gce = np.load('fermi_data/template_gce.npy')

n.add_template(dif, 'dif')
n.add_template(iso, 'iso')
n.add_template(bub, 'bub')
n.add_template(psc, 'psc')
n.add_template(gce, 'gce')

Step 3: Add Background Models to the Fit

Now from this list of templates the NPTF now knows about, we add in a series of background models which will be passed to MultiNest. In Example 7 we will show how to evaluate the likelihood without MultiNest, so that it can be interfaced with alternative inference packages.

Poissonian templates only have one parameter associated with them: \(A\) the template normalisation. Poissonian models are added to the fit via add_poiss_model. The first argument sets the spatial template for this background model, and should match the string used in add_template. The second argument is a LaTeX ready string used to identify the floated parameter later on.

By default added models will be floated. For floated templates the next two parameters are the prior range, added in the form [param_min, param_max] and then whether the prior is log flat (True) or linear flat (False). For log flat priors the priors are specified as indices, so that [-2,1] floats over a linear range [0.01,10].

Templates can also be added with a fixed normalisation. In this case no prior need be specified and instead fixed=True should be specified as well as fixed_norm=value, where value is \(A\) the template normalisation.

We use each of these possibilities in the example below.

n.add_poiss_model('dif', '$A_\mathrm{dif}$', False, fixed=True, fixed_norm=13.20)
n.add_poiss_model('iso', '$A_\mathrm{iso}$', [-2,1], True)
n.add_poiss_model('bub', '$A_\mathrm{bub}$', [0,2], False)
n.add_poiss_model('psc', '$A_\mathrm{psc}$', [0,2], False)
n.add_poiss_model('gce', '$A_\mathrm{gce}$', [0,2], False)

Note the diffuse model is normalised to a much larger value than the maximum prior of the other templates. This is because the diffuse model explains the majority of the flux in our ROI. The value of 13 was determined from a fit where the diffuse model was not fixed.

Step 4: Configure the Scan

Now the scan knows what models we want to fit to the data, we can configure the scan. In essence this step prepares all the information given above into an efficient format for calculating the likelihood. The main actions performed are: 1. Take the data and templates, and reduce them to only the ROI we will use as defined by the mask; 2. Further for a non-Poissonian scan an accounting for the number of exposure regions requested is made; and 3. Take the priors and parameters and prepare them into an efficient form for calculating the likelihood function that can then be used directly or passed to MultiNest.

n.configure_for_scan()

The number of parameters to be fit is 4

Step 5: Perform the Scan

Having setup all the parameters, we can now perform the scan using MultiNest. We will show an example of how to manually calculate the likelihood in Example 7.

Argument	Default Value	Purpose
run_tag	None	An additional tag can be specified to create a subdirectory of work_dir/chains/tag/ in which the output is stored.
nlive	100	Number of live points to be used during the MultiNest scan. A higher value thatn 100 is recommended for most runs, although larger values correspond to increased run time.
pymultinest_options	None	When set to None our default choices for MultiNest will be used (explained below). To alter these options, a dictionary of parameters and their values should be placed here.

Our default MultiNest options are defined as follows:

pymultinest_options = {'importance_nested_sampling': False,
                       'resume': False, 'verbose': True,
                       'sampling_efficiency': 'model',
                       'init_MPI': False, 'evidence_tolerance': 0.5,
                       'const_efficiency_mode': False}

For variations on these, a dictionary in the same format should be passed to perform_scan. A detailed explanation of the MultiNest options can be found here: https://johannesbuchner.github.io/PyMultiNest/pymultinest_run.html

n.perform_scan(nlive=500)

Step 6: Analyze the Output

Here we show a simple example of the output - the triangle plot. The full list of possible analysis options is explained in more detail in Example 9.

In order to do this we need to first load the scan using load_scan, which takes as an optional argument the same run_tag as used for the run. Note that load_scan can be used to load a run performed in a previous instance of NPTF, as long as the various parameters match.

After the scan is loaded we then create an instance of dnds_analysis, which takes an instance of nptfit.NPTF as an argument - which must already have a scan loaded. From here we simply make a triangle plot.

n.load_scan()
an = dnds_analysis.Analysis(n)
an.make_triangle()

analysing data from /zfs/nrodd/NPTFRemakeExamples/chains/Poissonian_Example/.txt

The triangle plot makes it clear that a non-zero value of the GCE template is preferred by the fit. Note also that as we gave the isotropic template a log flat prior, the parameter in the triangle plot is \(\log_{10} A_\mathrm{iso}\).

We also show the relative fraction of the Flux obtained by the GCE as compared to other templates. Note the majority of the flux is absorbed by the diffuse model.

an.plot_intensity_fraction_poiss('gce', bins=800, color='tomato', label='GCE')
an.plot_intensity_fraction_poiss('iso', bins=800, color='cornflowerblue', label='Iso')
an.plot_intensity_fraction_poiss('bub', bins=800, color='plum', label='Bub')
plt.xlabel('Flux fraction (%)')
plt.legend(fancybox = True)
plt.xlim(0,8);

Example 4: Simple NPTF example

In this example we perform a non-Poissonian template fit in a simplified setting. Specifically we will restrict ourselves to randomly generated nside=2 maps, which means our data consists only of 48 pixels. Nevertheless in this simple setting we will be able to clearly see the difference between Poissonian and non-Poissonian statistics as well as basic features of how non-Poissonian template fitting is performed with the code.

Throughout this example we will assume there is no smearing of the counts coming from all point sources. The effect of a finite point spread function on the statistics and how to account for it is discussed in Example 5.

# Import relevant modules

%matplotlib inline
%load_ext autoreload
%autoreload 2

import numpy as np
import healpy as hp
import matplotlib.pyplot as plt
from matplotlib import rcParams

from NPTFit import nptfit # module for performing scan
from NPTFit import dnds_analysis # module for analysing the output

from __future__ import print_function

Example 1: A map without point sources

We start out by analyzing a map without any point sources, using a uniform exposure map. First let’s create and plot our random data.

nside = 2
npix = hp.nside2npix(nside)
data = np.random.poisson(1,npix).astype(np.int32)
exposure = np.ones(npix)

hp.mollview(data,title='Fake Data')
hp.mollview(exposure,title='Exposure Map')

Next we set up an instance of NPTFit and add in the data. We’ll analyze the entire sky at once, so we won’t add in a mask.

n = nptfit.NPTF(tag='SimpleNPTF_Example')
n.load_data(data, exposure)

Now we add in templates, one to describe isotropic Poissonian emission and one for isotropic point sources. Note the different syntax requires for each.

iso = np.ones(npix)
n.add_template(iso, 'iso_p', units='flux')
n.add_template(iso, 'iso_np', units='PS')

We add in both models, being careful to select the right template. Here we model the PS point spread function as a singly broken power law, which requires four parameters to describe it: the normalization A, the indices above and below the breaks n1 and n2, and the location of the break Sb.

More details on the forms for the non-Poissonian model can be found in Example 6.

n.add_poiss_model('iso_p', '$A_\mathrm{iso}$', [0,2], False)
n.add_non_poiss_model('iso_np',
                      ['$A^\mathrm{ps}_\mathrm{iso}$','$n_1$','$n_2$','$S_b$'],
                      [[-10,1],[2.05,60],[-60,1.95],[0.01,200]],
                      [True,False,False,False])

Once everything is setup, we configure and perform the scan, and then show the triangle plot and flux fraction plot.

n.configure_for_scan()
n.perform_scan(nlive=500)

n.load_scan()
an = dnds_analysis.Analysis(n)
an.make_triangle()
plt.show()
plt.close()

an.plot_intensity_fraction_poiss('iso_p', bins=20, color='cornflowerblue', label='Poissonian')
an.plot_intensity_fraction_non_poiss('iso_np', bins=20, color='firebrick', label='non-Poissonian')
plt.xlabel('Flux fraction (\%)')
plt.legend(fancybox = True)
plt.xlim(0,100);
plt.ylim(0,0.4);

No mask set; defaulting to a blank mask
The number of parameters to be fit is 5
  analysing data from /zfs/nrodd/NPTFRemakeExamples/chains/SimpleNPTF_Example/.txt

We see that the Poissonian template has absorbed essentially everything, whilst the non-Poissonian parameters are poorly converged - expected as there were no point sources injected.

Example 2: A map with point sources

We now repeat the analysis above, but now add in 10 mean 50 count point sources. First lets take the data from above and add the point sources.

for ips in range(10):
    data[np.random.randint(npix)] += np.random.poisson(50)

hp.mollview(data,title='Fake Data with point sources')

Now we repeat all the steps used in the example without point sources.

n = nptfit.NPTF(tag='SimpleNPTF_Example')
n.load_data(data,exposure)

iso = np.ones(npix)
n.add_template(iso, 'iso_p',units='flux')
n.add_template(iso, 'iso_np',units='PS')

n.add_poiss_model('iso_p', '$A_\mathrm{iso}$', [0,2], False)
n.add_non_poiss_model('iso_np',
                      ['$A^\mathrm{ps}_\mathrm{iso}$','$n_1$','$n_2$','$S_b$'],
                      [[-10,1],[2.05,60],[-60,1.95],[0.01,200]],
                      [True,False,False,False])

n.configure_for_scan()
n.perform_scan(nlive=500)

n.load_scan()
an = dnds_analysis.Analysis(n)
an.make_triangle()
plt.show()
plt.close()

No mask set; defaulting to a blank mask
The number of parameters to be fit is 5
  analysing data from /zfs/nrodd/NPTFRemakeExamples/chains/SimpleNPTF_Example/.txt

We now see that both the Poissonian and non-Poissonian parameters are quite well converged. Note that the indices both want to have a large magnitude, which makes sense as we have effectively injected a delta function in flux, and the singly broken power law is trying to mimic that. Note that Sb is well converged near 50 counts per source, which is what we injected.

Example 3: A map with point sources and non-uniform exposure map

We will now repeat the above exercise but on a map without uniform exposure. This will highlight the importance of exposure regions.

To begin with let’s again create the data, but now we pretend that one side of the sky is expected to obtain twice as many counts as the other (which could occur if the instrument looked at that half of the sky twice as long for example).

nside = 2
npix = hp.nside2npix(nside)
exposure = np.zeros(npix)
exposure[0:int(npix/2)] = 1.0
exposure[int(npix/2):npix] = 2.0
data = np.random.poisson(exposure).astype(np.int32)

for ips in range(10):
    loc = np.random.randint(npix)
    data[loc] += np.random.poisson(50*exposure[loc])

hp.mollview(data,title='Fake Data with point sources')
hp.mollview(exposure,title='non-uniform Exposure Map')

Now we again analyze this data. Critically, note that when we configure the scan we set nexp=2 to indicate that the code should run with 2 exposure regions. In this simple example we know that 2 is all we need, but in real situations it is worth trying various values of nexp to see where results converge.

n = nptfit.NPTF(tag='SimpleNPTF_Example')
n.load_data(data,exposure)

iso = np.ones(npix)
n.add_template(iso, 'iso_p',units='flux')
n.add_template(iso, 'iso_np',units='PS')

n.add_poiss_model('iso_p', '$A_\mathrm{iso}$', [0,2], False)
n.add_non_poiss_model('iso_np',
                      ['$A^\mathrm{ps}_\mathrm{iso}$','$n_1$','$n_2$','$S_b$'],
                      [[-10,1],[2.05,60],[-60,1.95],[0.01,200]],
                      [True,False,False,False])

n.configure_for_scan(nexp=2)
n.perform_scan(nlive=500)

n.load_scan()
an = dnds_analysis.Analysis(n)
an.make_triangle()
plt.show()
plt.close()

No mask set; defaulting to a blank mask
The number of parameters to be fit is 5
  analysing data from /zfs/nrodd/NPTFRemakeExamples/chains/SimpleNPTF_Example/.txt

Everything is again well converged. Note that this time Sb has converged near 75, not 50. This is exactly what should be expected though, as the mean number of injected counts per PS over the sky is 50 x mean(exposure) = 75.

To highlight the importance of the exposure regions, let’s repeat this using only one exposure region which we emphasize is the wrong thing to do.

n = nptfit.NPTF(tag='SimpleNPTF_Example')
n.load_data(data,exposure)

iso = np.ones(npix)
n.add_template(iso, 'iso_p',units='flux')
n.add_template(iso, 'iso_np',units='PS')

n.add_poiss_model('iso_p', '$A_\mathrm{iso}$', [0,2], False)
n.add_non_poiss_model('iso_np',
                      ['$A^\mathrm{ps}_\mathrm{iso}$','$n_1$','$n_2$','$S_b$'],
                      [[-10,1],[2.05,60],[-60,1.95],[0.01,200]],
                      [True,False,False,False])

n.configure_for_scan(nexp=1)
n.perform_scan(nlive=500)

n.load_scan()
an = dnds_analysis.Analysis(n)
an.make_triangle()
plt.show()
plt.close()

No mask set; defaulting to a blank mask
The number of parameters to be fit is 5
  analysing data from /zfs/nrodd/NPTFRemakeExamples/chains/SimpleNPTF_Example/.txt

We see that the non-Poissonian parameters are not as well converged, and in particular Sb has not converged to its mean value over the sky. Note that Sb appears in some instance to be bimodal, at 50 and 100 representing the two point source brightnesses if the exposure correction is not accounted for.

Example 5: NPTF Correction for the Point Spread Function (PSF)

In this example we show how to account for the PSF correction using psf_correction.py

Fundamentally the presence of a non-zero PSF implies that the photons from any point source will be smeared out into some region around its true location. This effect must be accounted for in the NPTF. This is achieved via a function \(\rho(f)\). In the code we discretize \(\rho(f)\) as an approximation to the full function.

The two outputs of an instance of psf_correction are: 1. f_ary, an array of f values; and 2. df_rho_div_f_ary, an associated array of \(\Delta f \rho(f)/f\) values, where \(\Delta f\) is the width of the f_ary bins.

If the angular reconstruction of the data is perfect, then \(\rho(f) = \delta(f-1)\). In many situations, such as for the Fermi data at higher energies, a Gaussian approximation of the PSF will suffice. Even then there are a number of variables that go into evaluating the correction, as shown below. Finally we will show how the code can be used for the case of non-Gaussian PSFs.

As the calculation of \(\rho(f)\) can be time consuming, we always save the output to avoid recomputing the same correction twice. Consequently it can be convenient to have a common psf_dir where all PSF corrections for the runs are stored.

# Import relevant modules

%matplotlib inline
%load_ext autoreload
%autoreload 2

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams

from NPTFit import psf_correction as pc # Module for determining the PSF correction

from __future__ import print_function

Example 1: Default Gaussian PSF

We start by showing the PSF correction for a Gaussian PSF - that is the PSF as a function of \(r\) is \(\exp \left[ -r^2 / (2\sigma^2) \right]\) - with \(\sigma\) set to the value of the 68% containment radius for the PSF of the Fermi dataset we will use in later examples.

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.1812)
f_ary_1 = pc_inst.f_ary
df_rho_div_f_ary_1 = pc_inst.df_rho_div_f_ary

print('f_ary:', f_ary_1)
print('df_rho_div_f_ary:', df_rho_div_f_ary_1)

plt.plot(f_ary_1,f_ary_1**2*df_rho_div_f_ary_1/(f_ary_1[1]-f_ary_1[0]),color='black', lw = 1.5)
plt.xlabel('$f$')
plt.ylabel('$f \\times \\rho(f)$')
plt.title('Gaussian PSF, $\sigma_\mathrm{PSF} = 0.1812$', y=1.04)

File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.181_10_50000_1000_0.01.npy
f_ary: [0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95]
df_rho_div_f_ary: [65.19815984  6.88897747  2.52908225  1.28920055  0.80522461  0.54317562
  0.09145394  0.          0.          0.        ]

Text(0.5,1.04,'Gaussian PSF, $\sigma_\mathrm{PSF} = 0.1812$')

Example 2: Impact of changing \(\sigma\)

Here we show the impact on the PSF of changing \(\sigma\). From the plot we can see that for a small PSF, \(\rho(f)\) approaches the no PSF case of \(\delta(f-1)\) implying that the flux fractions are concentrated at a single large value. As \(\sigma\) increases we move away from this idealized scenario and the flux becomes more spread out, leading to a \(\rho(f)\) peaked at lower flux values.

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.05)
f_ary_2 = pc_inst.f_ary
df_rho_div_f_ary_2 = pc_inst.df_rho_div_f_ary

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.4)
f_ary_3 = pc_inst.f_ary
df_rho_div_f_ary_3 = pc_inst.df_rho_div_f_ary

plt.plot(f_ary_1,f_ary_1**2*df_rho_div_f_ary_1/(f_ary_1[1]-f_ary_1[0]),color='cornflowerblue',label='0.18', lw = 1.5)
plt.plot(f_ary_2,f_ary_2**2*df_rho_div_f_ary_2/(f_ary_2[1]-f_ary_2[0]),color='forestgreen',label='0.05', lw = 1.5)
plt.plot(f_ary_3,f_ary_3**2*df_rho_div_f_ary_3/(f_ary_3[1]-f_ary_3[0]),color='maroon',label='0.4', lw = 1.5)
plt.xlabel('$f$')
plt.ylabel('$f \\times \\rho(f)$')
plt.legend(loc='upper right', fancybox=True)
plt.title('Varying $\sigma_\mathrm{PSF}$', y=1.04)

File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.05_10_50000_1000_0.01.npy
File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.4_10_50000_1000_0.01.npy

Text(0.5,1.04,'Varying $\sigma_\mathrm{PSF}$')

Example 3: Changing the default options for determining \(\rho(f)\)

In this example we show how for a given PSF, the other parameters associated with how accurately we calculate \(\rho(f)\) can impact what we get back. The parameters that can be changed are are:

Argument	Defaults	Purpose
num_f_bins	10	number of f_bins used
n_psf	50000	number of PSFs placed down when calculating
n_pts_per_psf	1000	number of points to place per psf in the calculation
f_trunc	0.01	minimum flux fraction to keep track of
nside	128	nside of the map the PSF is used on

The default parameters have been chosen to be accurate enough for the Fermi analyses we will be performed later. But if the user changes the PSF (even just \(\sigma\)), it is important to be sure that the above parameters are chosen so that \(\rho(f)\) is evaluated accurately enough.

In general increasing num_f_bins, n_psf, and n_pts_per_psf, whilst decreasing f_trunc leads to a more accurate \(\rho(f)\). But each will also slow down the evaluation of \(\rho(f)\), and in the case of num_f_bin, slow down the subsequent non-Poissonian likelihood evaluation.

nside should be set to the value of the map being analysed, but we also highlight the impact of changing it below. For an analysis on a non-HEALPix grid, the PSF can often be approximated by an appropriate HEALPix binning. If this is not the case, however, a different approach must be pursued in calculating \(\rho(f)\).

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.1812,num_f_bins=20)
f_ary_4 = pc_inst.f_ary
df_rho_div_f_ary_4 = pc_inst.df_rho_div_f_ary

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.1812,n_psf=5000,n_pts_per_psf=100)
f_ary_5 = pc_inst.f_ary
df_rho_div_f_ary_5 = pc_inst.df_rho_div_f_ary

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.1812,f_trunc=0.1)
f_ary_6 = pc_inst.f_ary
df_rho_div_f_ary_6 = pc_inst.df_rho_div_f_ary

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.1812,nside=64)
f_ary_7 = pc_inst.f_ary
df_rho_div_f_ary_7 = pc_inst.df_rho_div_f_ary

plt.plot(f_ary_1,f_ary_1**2*df_rho_div_f_ary_1/(f_ary_1[1]-f_ary_1[0]),color='black',label=r'Default', lw=2.2)
plt.plot(f_ary_4,f_ary_4**2*df_rho_div_f_ary_4/(f_ary_4[1]-f_ary_4[0]),color='forestgreen',label=r'more f\_bins', lw = 1.5)
plt.plot(f_ary_5,f_ary_5**2*df_rho_div_f_ary_5/(f_ary_5[1]-f_ary_5[0]),color='cornflowerblue',label=r'fewer points', lw = 1.5)
plt.plot(f_ary_6,f_ary_6**2*df_rho_div_f_ary_6/(f_ary_6[1]-f_ary_6[0]),color='salmon',label=r'larger f\_trunc', lw = 1.5)
plt.plot(f_ary_7,f_ary_7**2*df_rho_div_f_ary_7/(f_ary_7[1]-f_ary_7[0]),color='orchid',label=r'lower nside', lw = 1.5)
plt.xlabel('$f$')
plt.ylabel('$f \\times \\rho(f)$')
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5), fancybox=True)

File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.181_20_50000_1000_0.01.npy
File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.181_10_5000_100_0.01.npy
File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.181_10_50000_1000_0.1.npy
File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_64_0.181_10_50000_1000_0.01.npy

<matplotlib.legend.Legend at 0x7fe4c00337b8>

Example 4: PSF on a Cartesian Grid

For some applications, particularly when analyzing smaller regions of the sky, it may be desirable to work with data on a Cartesian grid rather than a healpix map. Note generally for larger regions, in order to account for curvature on the sky a healpix pixelization is recommended. Code to convert from Cartesian grids to healpix can be found here: https://github.com/nickrodd/grid2healpix

In order to calculate the appropriate PSF correction for Cartesian maps the general syntax is the same, except now the healpix_map keyword should be set to False and the pixarea keyword set to the area in sr of each pixel of the Cartesian map. In addition the gridsize keyword determines how large the map is, and flux that falls outside the map is lost in the Cartesian case.

As an example of this syntax we calculate the PSF correction on a Cartesian map that has pixels the same size as an nside=128 healpix map, and compare the two PSF corrections. Note they are essentially identical.

pixarea = 4*np.pi/(12*128*128)
pc_inst = pc.PSFCorrection(psf_sigma_deg=0.1812, healpix_map=False, pixarea=pixarea, gridsize=100)
f_ary_8 = pc_inst.f_ary
df_rho_div_f_ary_8 = pc_inst.df_rho_div_f_ary

plt.plot(f_ary_1,f_ary_1**2*df_rho_div_f_ary_1/(f_ary_1[1]-f_ary_1[0]),color='black', label=r'healpix', lw = 1.5)
plt.plot(f_ary_8,f_ary_8**2*df_rho_div_f_ary_8/(f_ary_8[1]-f_ary_8[0]),color='forestgreen', label=r'cartesian', lw = 1.5)
plt.xlabel('$f$')
plt.ylabel('$f \\times \\rho(f)$')

File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_0.21_0.181_10_50000_1000_0.01.npy

Text(0,0.5,'$f \times \rho(f)$')

Example 5: Custom PSF

In addition to the default Gausian PSF, psf_correction.py also has the option of taking in a custom PSF. In order to use this ability, the user needs to initialise psf_correction with delay_compute=True, manually define the parameters that define the PSF and then call make_or_load_psf_corr.

The variables that need to be redefined in the instance of psf_correction are:

Argument	Purpose
psf_r_func	the psf as a function of r, distance in radians from the center of the point source
sample_psf_max	maximum distance to sample the psf from the center, should be larger for psfs with long tails
psf_samples	number of samples to make from the psf (linearly spaced) from 0 to sample_psf_max, should be large enough to adequately represent the full psf
psf_tag	label the psf is saved with

As an example of a more complicated PSF we consider the full Fermi-LAT PSF. The PSF of Fermi is approximately Gaussian near the core, but has larger tails. To model this a pair of King functions are used to describe the radial distribution. Below we show a comparison between the Gaussian approximation and the full PSF, for two different energies. As shown, for low energies where the Fermi PSF is larger, the difference between the two can be significant. For higher energies where the PSF becomes smaller, however, the difference is marginal.

For the full details of the Fermi-LAT PSF, see: http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_LAT_IRFs/IRF_PSF.html

# Fermi-LAT PSF at 2 GeV

# Calculate the appropriate Gaussian approximation to the PSF for 2 GeV
pc_inst = pc.PSFCorrection(psf_sigma_deg=0.2354)
f_ary_9 = pc_inst.f_ary
df_rho_div_f_ary_9 = pc_inst.df_rho_div_f_ary

# Define parameters that specify the Fermi-LAT PSF at 2 GeV
fcore = 0.748988248179
score = 0.428653790656
gcore = 7.82363229341
stail = 0.715962650769
gtail = 3.61883748683
spe = 0.00456544262478

# Define the full PSF in terms of two King functions
def king_fn(x, sigma, gamma):
    return 1./(2.*np.pi*sigma**2.)*(1.-1./gamma)*(1.+(x**2./(2.*gamma*sigma**2.)))**(-gamma)

def Fermi_PSF(r):
    return fcore*king_fn(r/spe,score,gcore) + (1-fcore)*king_fn(r/spe,stail,gtail)

# Modify the relevant parameters in pc_inst and then make or load the PSF
pc_inst = pc.PSFCorrection(delay_compute=True)
pc_inst.psf_r_func = lambda r: Fermi_PSF(r)
pc_inst.sample_psf_max = 10.*spe*(score+stail)/2.
pc_inst.psf_samples = 10000
pc_inst.psf_tag = 'Fermi_PSF_2GeV'
pc_inst.make_or_load_psf_corr()

# Extract f_ary and df_rho_div_f_ary as usual
f_ary_10 = pc_inst.f_ary
df_rho_div_f_ary_10 = pc_inst.df_rho_div_f_ary

plt.plot(f_ary_9,f_ary_9**2*df_rho_div_f_ary_9/(f_ary_9[1]-f_ary_9[0]),color='maroon',label='Gauss PSF', lw = 1.5)
plt.plot(f_ary_10,f_ary_10**2*df_rho_div_f_ary_10/(f_ary_10[1]-f_ary_10[0]),color='forestgreen',label='Fermi PSF', lw = 1.5)
plt.xlabel('$f$')
plt.ylabel('$f \\times \\rho(f)$')
plt.legend(loc='upper right', fancybox=True)

File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.235_10_50000_1000_0.01.npy
File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/Fermi_PSF_2GeV.npy

<matplotlib.legend.Legend at 0x7fe4b5217e10>

# Fermi-LAT PSF at 20 GeV

# Calculate the appropriate Gaussian approximation to the PSF for 20 GeV
pc_inst = pc.PSFCorrection(psf_sigma_deg=0.05529)
f_ary_11 = pc_inst.f_ary
df_rho_div_f_ary_11 = pc_inst.df_rho_div_f_ary

# Define parameters that specify the Fermi-LAT PSF at 20 GeV
fcore = 0.834725201378
score = 0.498192326976
gcore = 6.32075520959
stail = 1.06648424558
gtail = 4.49677834267
spe = 0.000943339426754

# Define the full PSF in terms of two King functions
def king_fn(x, sigma, gamma):
    return 1./(2.*np.pi*sigma**2.)*(1.-1./gamma)*(1.+(x**2./(2.*gamma*sigma**2.)))**(-gamma)

def Fermi_PSF(r):
    return fcore*king_fn(r/spe,score,gcore) + (1-fcore)*king_fn(r/spe,stail,gtail)

# Modify the relevant parameters in pc_inst and then make or load the PSF
pc_inst = pc.PSFCorrection(delay_compute=True)
pc_inst.psf_r_func = lambda r: Fermi_PSF(r)
pc_inst.sample_psf_max = 10.*spe*(score+stail)/2.
pc_inst.psf_samples = 10000
pc_inst.psf_tag = 'Fermi_PSF_20GeV'
pc_inst.make_or_load_psf_corr()

# Extract f_ary and df_rho_div_f_ary as usual
f_ary_12 = pc_inst.f_ary
df_rho_div_f_ary_12 = pc_inst.df_rho_div_f_ary

plt.plot(f_ary_11,f_ary_11**2*df_rho_div_f_ary_11/(f_ary_11[1]-f_ary_11[0]),color='maroon',label='Gauss PSF', lw = 1.5)
plt.plot(f_ary_12,f_ary_12**2*df_rho_div_f_ary_12/(f_ary_12[1]-f_ary_12[0]),color='forestgreen',label='Fermi PSF', lw = 1.5)
plt.xlabel('$f$')
plt.ylabel('$f \\times \\rho(f)$')
plt.legend(loc='upper left', fancybox=True)

File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.055_10_50000_1000_0.01.npy
File saved as: /zfs/nrodd/NPTFRemakeExamples/psf_dir/Fermi_PSF_20GeV.npy

<matplotlib.legend.Legend at 0x7fe4b51cccc0>

The above example also serves as a tutorial on how to combine various PSFs into a single PSF. In the case of the Fermi PSF the full radial dependence is the sum of two King functions. More generally if the full PSF is a combination of multiple individual ones (for example from multiple energy bins), then this can be formed by just adding these functions with an appropriate weighting to get a single psf_r_func.

Example 6: Running non-Poissonian Scans with MultiNest

In this example we show a simple example of how to run a scan using templates that follow non-Poisson statistics.

To this end we perform a simple analysis of a small region around the north galactic pole, looking for isotropically distributed point sources.

NB: This example makes use of the Fermi Data, which needs to already be installed. See Example 1 for details.

# Import relevant modules

%matplotlib inline
%load_ext autoreload
%autoreload 2

import numpy as np
import corner
import matplotlib.pyplot as plt
from matplotlib import rcParams

from NPTFit import nptfit # module for performing scan
from NPTFit import create_mask as cm # module for creating the mask
from NPTFit import psf_correction as pc # module for determining the PSF correction
from NPTFit import dnds_analysis # module for analysing the output

Step 1: Setup an instance of NPTFit, add data, mask and templates

As in the Poissonian case (see Example 3), the first step is to create an instance of nptfit.NPTF and add to it data, a mask and templates. Here we choose an ROI that is a small ring around the galactic north pole. In this small region we will only use isotropically distributed templates.

n = nptfit.NPTF(tag='non-Poissonian_Example')

fermi_data = np.load('fermi_data/fermidata_counts.npy').astype(np.int32)
fermi_exposure = np.load('fermi_data/fermidata_exposure.npy')
n.load_data(fermi_data, fermi_exposure)

analysis_mask = cm.make_mask_total(mask_ring = True, inner = 0, outer = 5, ring_b = 90, ring_l = 0)
n.load_mask(analysis_mask)

Now we add the required templates. The iso_p template here is in units of photon counts and is exposure corrected, which we anticipate using as such for the Poissonian model, while for the non-Poissonian model the underlying PS distribution is truly isotropic.

iso_p = np.load('fermi_data/template_iso.npy')
n.add_template(iso_p, 'iso_p')
iso_np = np.ones(len(iso_p))
n.add_template(iso_np, 'iso_np',units='PS')

Step 2: Add Background Models to the Fit

We now add background models. For the details of adding a Poissonian model, see Example 3.

For the non-poissonian model the format is similar, except for the fact we must account for the additional parameters that describe an non-Poissonian (NP) template. In nptfit, NP templates are determined by source count functions which we allow to be specified by a broken power law with an arbitrary number of breaks.

The simplest possible example consistent with a finite number of sources is a singly broken power law, which is specified by four parameters: the template normalisation \(A\), indices above and below the break \(n_1\) and \(n_2\), and the location of the break \(S_b\).

In general, for a source count function with \(\ell\) breaks, the \(2\ell+2\) parameters are specified as follows:

\[\left[ A, n_1, \ldots, n_{\ell+1}, S_b^{(1)}, \ldots, S_b^{(\ell)} \right]\,,\]

where \(n_1\) is the highest index and \(S_b^{(1)}\) the highest break.

Priors must be entered as an array where each element is an array of the priors for each unfixed parameter. For multiply broken power laws it is possible to specify the breaks in terms of the highest break, in which case the option dnds_model=specifiy_relative_breaks should be used.

Fixed parameters are similarly entered as an array, where the first element is the location of the parameter to be fixed (an integer), and the second element is the value to which it should be fixed.

In the example below we add an isotropic distributed non-Poissonian template, with a log flat normalisation, linear flat indices, and a fixed break.

n.add_poiss_model('iso_p','$A_\mathrm{iso}$', False, fixed=True, fixed_norm=1.51)
n.add_non_poiss_model('iso_np',
                      ['$A^\mathrm{ps}_\mathrm{iso}$','$n_1$','$n_2$','$S_b$'],
                      [[-6,1],[2.05,30],[-2,1.95]],
                      [True,False,False],
                      fixed_params = [[3,172.52]])

Step 3: Configure the Scan with the PSF correction

For a non-Poissonian fit, we need to specify the PSF correction at the stage of configuring the scan. The details of this are described in Example 5. These are calculated using psf_correction.py and then passed to the NPTF via configure_for_scan.

At this stage we also specify the number of exposure regions to be used. Here we take nexp=1 for a simple example. Generally increasing nexp leads to more accurate results, but also increases the runtime of the code.

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.1812)
f_ary = pc_inst.f_ary
df_rho_div_f_ary = pc_inst.df_rho_div_f_ary

n.configure_for_scan(f_ary=f_ary, df_rho_div_f_ary=df_rho_div_f_ary, nexp=1)

Loading the psf correction from: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.181_10_50000_1000_0.01.npy
The number of parameters to be fit is 3

Step 4: Perform the Scan

Next we perform the scan. The syntax is identical to the Poissonian case and described in Example 3. Note even though we float fewer parameters than in the Poissonian example, the runtime is longer here. This is due to the fact that the NPTF likelihood is inherently more complicated and so takes longer to evaluate.

n.perform_scan(nlive=800)

Step 5: Analyze the Output

Here we analyze the output using the same commands as in the Poissonian example.

n.load_scan()
cs=dnds_analysis.Analysis(n)
cs.make_triangle()

analysing data from /zfs/nrodd/NPTFRemakeExamples/chains/non-Poissonian_Example/.txt

/zfs/nrodd/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_base.py:2961: UserWarning: Attempting to set identical left==right results
in singular transformations; automatically expanding.
left=172.52, right=172.52
  'left=%s, right=%s') % (left, right))
/zfs/nrodd/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_base.py:3285: UserWarning: Attempting to set identical bottom==top results
in singular transformations; automatically expanding.
bottom=172.52, top=172.52
  'bottom=%s, top=%s') % (bottom, top))

We also show a plot of the source count function, although a careful explanation of the details here are deferred until Example 9.

cs.plot_source_count_median('iso_np',smin=0.01,smax=10000,nsteps=1000,spow=2,color='forestgreen')
cs.plot_source_count_band('iso_np',smin=0.01,smax=10000,nsteps=1000,qs=[0.16,0.5,0.84],spow=2,color='forestgreen',alpha=0.3)

plt.yscale('log')
plt.xscale('log')
plt.xlim([1e-10,1e-7])
plt.ylim([1e-15,1e-9])
plt.tick_params(axis='x', length=5,width=2,labelsize=18)
plt.tick_params(axis='y',length=5,width=2,labelsize=18)
plt.ylabel('$F^2 dN/dF$ [counts cm$^{-2}$s$^{-1}$deg$^{-2}$]', fontsize=18)
plt.xlabel('$F$  [counts cm$^{-2}$ s$^{-1}$]', fontsize=18)

Text(0.5,0,'$F$  [counts cm$^{-2}$ s$^{-1}$]')

Example 7: Manual evaluation of non-Poissonian Likelihood

In this example we show to manually evaluate the non-Poissonian likelihood. This can be used, for example, to interface nptfit with parameter estimation packages other than MultiNest. We also show how to extract the prior cube.

We will take the exact same analysis as considered in the previous example, and show the likelihood peaks at exactly the same location for the normalisation of the non-Poissonian template.

NB: This example makes use of the Fermi Data, which needs to already be installed. See Example 1 for details.

# Import relevant modules

%matplotlib inline
%load_ext autoreload
%autoreload 2

import numpy as np
import healpy as hp
import matplotlib.pyplot as plt

from NPTFit import nptfit # module for performing scan
from NPTFit import create_mask as cm # module for creating the mask
from NPTFit import psf_correction as pc # module for determining the PSF correction
from NPTFit import dnds_analysis # module for analysing the output

from __future__ import print_function

Setup an identical instance of NPTFit to Example 6

Firstly we initialize an instance of nptfit identical to that used in the previous example.

n = nptfit.NPTF(tag='non-Poissonian_Example')

fermi_data = np.load('fermi_data/fermidata_counts.npy').astype(np.int32)
fermi_exposure = np.load('fermi_data/fermidata_exposure.npy')
n.load_data(fermi_data, fermi_exposure)

analysis_mask = cm.make_mask_total(mask_ring = True, inner = 0, outer = 5, ring_b = 90, ring_l = 0)
n.load_mask(analysis_mask)

iso_p = np.load('fermi_data/template_iso.npy')
n.add_template(iso_p, 'iso_p')
iso_np = np.ones(len(iso_p))
n.add_template(iso_np, 'iso_np',units='PS')

n.add_poiss_model('iso_p','$A_\mathrm{iso}$', False, fixed=True, fixed_norm=1.51)
n.add_non_poiss_model('iso_np',
                      ['$A^\mathrm{ps}_\mathrm{iso}$','$n_1$','$n_2$','$S_b$'],
                      [[-6,1],[2.05,30],[-2,1.95]],
                      [True,False,False],
                      fixed_params = [[3,172.52]])

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.1812)
f_ary = pc_inst.f_ary
df_rho_div_f_ary = pc_inst.df_rho_div_f_ary

Loading the psf correction from: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.181_10_50000_1000_0.01.npy

n.configure_for_scan(f_ary=f_ary, df_rho_div_f_ary=df_rho_div_f_ary, nexp=1)

The number of parameters to be fit is 3

Evaluate the Likelihood Manually

After configuring for the scan, the instance of nptfit.NPTF now has an associated function ll. This function was passed to MultiNest in the previous example, but we can also manually evaluate it.

The log likelihood function is called as: ll(theta), where theta is a flattened array of parameters. In the case above:

\[\theta = \left[ \log_{10} \left( A^\mathrm{ps}_\mathrm{iso} \right), n_1, n_2 \right]\]

As an example we can evaluate it at a few points around the best fit parameters:

print('Vary A: ', n.ll([-4.76+0.32,18.26,0.06]), n.ll([-4.76,18.26,0.06]), n.ll([-4.76-0.37,18.26,0.06]))
print('Vary n1:', n.ll([-4.76,18.26+7.98,0.06]), n.ll([-4.76,18.26,0.06]), n.ll([-4.76,18.26-9.46,0.06]))
print('Vary n2:', n.ll([-4.76,18.26,0.06+0.93]), n.ll([-4.76,18.26,0.06]), n.ll([-4.76,18.26,0.06-1.31]))

Vary A:  -587.1163770023244 -586.2113098768381 -588.2406707060421
Vary n1: -586.1824552299396 -586.2113098768381 -586.3726977861785
Vary n2: -587.0768970698804 -586.2113098768381 -587.5645486197602

To make the point clearer we can fix \(n_1\) and \(n_2\) to their best fit values, and calculate a Test Statistics (TS) array as we vary \(\log_{10} \left( A^\mathrm{ps}_\mathrm{iso} \right)\). As shown the likelihood is maximised at approximated where MultiNest told us was the best fit point for this parameter.

Avals = np.arange(-6.0,-2.0,0.01)
TSvals_A = np.array([2*(n.ll([-4.76,18.26,0.06])-n.ll([Avals[i],18.26,0.06])) for i in range(len(Avals))])

plt.plot(Avals,TSvals_A,color='black', lw=1.5)
plt.axvline(-4.76+0.32,ls='dashed',color='black')
plt.axvline(-4.76,ls='dashed',color='black')
plt.axvline(-4.76-0.37,ls='dashed',color='black')
plt.axhline(0,ls='dashed',color='black')
plt.xlim([-5.5,-4.0])
plt.ylim([-5.0,15.0])
plt.xlabel('$A^\mathrm{ps}_\mathrm{iso}$')
plt.ylabel('$\mathrm{TS}$')
plt.show()

Next we do the same thing for \(n_2\). This time we see that this parameter is much more poorly constrained than the value of the normalisation, as the TS is very flat.

NB: it is important not to evaluate breaks exactly at a value of \(n=1\). The reason for this is the analytic form of the likelihood involves \((n-1)^{-1}\).

n2vals = np.arange(-1.995,1.945,0.01)
TSvals_n2 = np.array([2*(n.ll([-4.76,18.26,0.06])-n.ll([-4.76,18.26,n2vals[i]])) for i in range(len(n2vals))])

plt.plot(n2vals,TSvals_n2,color='black', lw=1.5)
plt.axvline(0.06+0.93,ls='dashed',color='black')
plt.axvline(0.06,ls='dashed',color='black')
plt.axvline(0.06-1.31,ls='dashed',color='black')
plt.axhline(0,ls='dashed',color='black')
plt.xlim([-2.0,1.5])
plt.ylim([-5.0,15.0])
plt.xlabel('$n_2$')
plt.ylabel('$\mathrm{TS}$')
plt.show()

In general \(\theta\) will always be a flattened array of the floated parameters. Poisson parameters always occur first, in the order in which they were added (via add_poiss_model), following by non-Poissonian parameters in the order they were added (via add_non_poiss_model). To be explicit if we have \(m\) Poissonian templates and \(n\) non-Poissonian templates with breaks \(\ell_n\), then:

\[\theta = \left[ A_\mathrm{P}^1, \ldots, A_\mathrm{P}^m, A_\mathrm{NP}^1, n_1^1, \ldots, n_{\ell_1+1}^1, S_b^{(1)~1}, \ldots, S_b^{(\ell_1)~1}, \ldots, A_\mathrm{NP}^n, n_1^n, \ldots, n_{\ell_n+1}^n, S_b^{(1)~n}, \ldots, S_b^{(\ell_n)~n} \right]\]

Fixed parameters are deleted from the list, and any parameter entered with a log flat prior is replaced by \(\log_{10}\) of itself.

Extract the Prior Cube Manually

To extract the prior cube, we use the internal function log_prior_cube. This requires two arguments: 1. cube, the unit cube of dimension equal to the number of floated parameters; and 2. ndim, the number of floated parameters.

print(n.prior_cube(cube=[1,1,1],ndim=3))

[1.0, 30.0, 1.9500000000000002]

Example 8: Application of NPTFit to the Galactic Center Excess

It was found in Example 3 that a non-zero value of the GCE template is preferred by a fit in the galactic center. In this example we will test the point source interpretation of this excess by including, in addition to the Poissonian templates considered, non-Poissonian point source templates of various morphologies.

Here we simply perform the run, a detailed analysis of the results can be found in the next example. A Python script version of this notebook can be found as Example8_Galactic_Center_nonPoissonian.py, which can be run faster on multiple processors with MPI (see example in Example8_Galactic_Center_Batch.batch.

NB: Even with nlive=100, this notebook can take up to an hour to complete. This highlights that for realistic non-Poissonian runs, running on multiple cores becomes necessary. We show an explicit application of this in Example 10.

NB: This example makes use of the Fermi Data, which needs to already be installed. See Example 1 for details.

# Import relevant modules

%matplotlib inline
%load_ext autoreload
%autoreload 2

import numpy as np

from NPTFit import nptfit # module for performing scan
from NPTFit import create_mask as cm # module for creating the mask
from NPTFit import dnds_analysis # module for analysing the output
from NPTFit import psf_correction as pc # module for determining the PSF correction

Step 1: Set up the Scan

We first need to 1. Set up an instance of NPTF from npfit.py 2. Load in the data and exposure maps 3. Set up and load the mask used for the scan 4. Load in the spatial templates

These are done identically to Example 3, and we refer to that notebook for details.

n = nptfit.NPTF(tag='GCE_Example')

fermi_data = np.load('fermi_data/fermidata_counts.npy').astype(np.int32)
fermi_exposure = np.load('fermi_data/fermidata_exposure.npy')
n.load_data(fermi_data, fermi_exposure)

pscmask=np.array(np.load('fermi_data/fermidata_pscmask.npy'), dtype=bool)
analysis_mask = cm.make_mask_total(band_mask = True, band_mask_range = 2,
                                   mask_ring = True, inner = 0, outer = 30,
                                   custom_mask = pscmask)
n.load_mask(analysis_mask)

dif = np.load('fermi_data/template_dif.npy')
iso = np.load('fermi_data/template_iso.npy')
bub = np.load('fermi_data/template_bub.npy')
gce = np.load('fermi_data/template_gce.npy')
dsk = np.load('fermi_data/template_dsk.npy')

n.add_template(dif, 'dif')
n.add_template(iso, 'iso')
n.add_template(bub, 'bub')
n.add_template(gce, 'gce')
n.add_template(dsk, 'dsk')

# Remove the exposure correction for PS templates
rescale = fermi_exposure/np.mean(fermi_exposure)
n.add_template(gce/rescale, 'gce_np', units='PS')
n.add_template(dsk/rescale, 'dsk_np', units='PS')

Step 2: Add Models

n.add_poiss_model('dif', '$A_\mathrm{dif}$', fixed=True, fixed_norm=12.85)
n.add_poiss_model('iso', '$A_\mathrm{iso}$', [0,2], False)
n.add_poiss_model('gce', '$A_\mathrm{gce}$', [0,2], False)
n.add_poiss_model('bub', '$A_\mathrm{bub}$', [0,2], False)

This time we add a non-Poissonian template correlated with the Galactic Center Excess and also one spatially distributed as a thin disk. The latter is designed to account for the unresolved point sources attributed to the disk of the Milky Way (known sources in the 3FGL are masked).

n.add_non_poiss_model('gce_np',
                      ['$A_\mathrm{gce}^\mathrm{ps}$','$n_1^\mathrm{gce}$','$n_2^\mathrm{gce}$','$S_b^{(1), \mathrm{gce}}$'],
                      [[-6,1],[2.05,30],[-2,1.95],[0.05,40]],
                      [True,False,False,False])
n.add_non_poiss_model('dsk_np',
                      ['$A_\mathrm{dsk}^\mathrm{ps}$','$n_1^\mathrm{dsk}$','$n_2^\mathrm{dsk}$','$S_b^{(1), \mathrm{dsk}}$'],
                      [[-6,1],[2.05,30],[-2,1.95],[0.05,40]],
                      [True,False,False,False])

Step 3: Configure Scan with PSF correction

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.1812)
f_ary, df_rho_div_f_ary = pc_inst.f_ary, pc_inst.df_rho_div_f_ary

Loading the psf correction from: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.181_10_50000_1000_0.01.npy

n.configure_for_scan(f_ary, df_rho_div_f_ary, nexp=1)

The number of parameters to be fit is 11

Step 4: Perform the Scan

As noted above, we take a small value of nlive simply to ensure the run finishes in a reasonable time on a single core.

n.perform_scan(nlive=100)

This can take up to an hour to run. The output of this run will be analyzed in detail in the next example.

from IPython.display import Image
Image(url = "https://imgs.xkcd.com/comics/compiling.png")

Example 9: Analyzing the Results of an NPTFit Run

While the chain samples of a non-Poissonian fit performed using MultiNest can be readily accessed, we provide a basic analysis module, dnds_analysis.py that contains helper functions to: 1. Make triangle plots; 2. Get template intensities; 3. Plot source count distributions; 4. Plot flux fractions; 5. Access individual posteriors; and 6. Get Bayesian log-evidences.

In this example we provide the details of how to use each function.

NB: Example 8 must be run before this notebook. Note that the run performed there was with a low nside and a fixed diffuse model, so the results below should be interpreted only as approximate.

# Import relevant modules

%matplotlib inline
%load_ext autoreload
%autoreload 2

import numpy as np
import corner
import matplotlib.pyplot as plt

from NPTFit import nptfit # module for performing scan
from NPTFit import create_mask as cm # module for creating the mask
from NPTFit import dnds_analysis # module for analysing the output
from NPTFit import psf_correction as pc # module for determining the PSF correction

from __future__ import print_function

Analysis

At the outset, an instance of nptfit.NPTF must be created as was done when initiating and performing the scan. The process here is the same as in Example 8, up to configuring the scan. Finally, the scan is loaded with n.load_scan().

n = nptfit.NPTF(tag='GCE_Example')

fermi_data = np.load('fermi_data/fermidata_counts.npy').astype(np.int32)
fermi_exposure = np.load('fermi_data/fermidata_exposure.npy')
n.load_data(fermi_data, fermi_exposure)

pscmask=np.array(np.load('fermi_data/fermidata_pscmask.npy'), dtype=bool)
analysis_mask = cm.make_mask_total(band_mask = True, band_mask_range = 2,
                                   mask_ring = True, inner = 0, outer = 30,
                                   custom_mask = pscmask)
n.load_mask(analysis_mask)

dif = np.load('fermi_data/template_dif.npy')
iso = np.load('fermi_data/template_iso.npy')
bub = np.load('fermi_data/template_bub.npy')
gce = np.load('fermi_data/template_gce.npy')
dsk = np.load('fermi_data/template_dsk.npy')

n.add_template(dif, 'dif')
n.add_template(iso, 'iso')
n.add_template(bub, 'bub')
n.add_template(gce, 'gce')
n.add_template(dsk, 'dsk')

# Remove the exposure correction for PS templates
rescale = fermi_exposure/np.mean(fermi_exposure)
n.add_template(gce/rescale, 'gce_np', units='PS')
n.add_template(dsk/rescale, 'dsk_np', units='PS')

n.add_poiss_model('dif', '$A_\mathrm{dif}$', fixed=True, fixed_norm=12.85)
n.add_poiss_model('iso', '$A_\mathrm{iso}$', [0,2], False)
n.add_poiss_model('gce', '$A_\mathrm{gce}$', [0,2], False)
n.add_poiss_model('bub', '$A_\mathrm{bub}$', [0,2], False)

n.add_non_poiss_model('gce_np',
                      ['$A_\mathrm{gce}^\mathrm{ps}$','$n_1^\mathrm{gce}$','$n_2^\mathrm{gce}$','$S_b^{(1), \mathrm{gce}}$'],
                      [[-6,1],[2.05,30],[-2,1.95],[0.05,40]],
                      [True,False,False,False])
n.add_non_poiss_model('dsk_np',
                      ['$A_\mathrm{dsk}^\mathrm{ps}$','$n_1^\mathrm{dsk}$','$n_2^\mathrm{dsk}$','$S_b^{(1), \mathrm{dsk}}$'],
                      [[-6,1],[2.05,30],[-2,1.95],[0.05,40]],
                      [True,False,False,False])

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.1812)
f_ary, df_rho_div_f_ary = pc_inst.f_ary, pc_inst.df_rho_div_f_ary

Loading the psf correction from: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.181_10_50000_1000_0.01.npy

n.configure_for_scan(f_ary, df_rho_div_f_ary, nexp=1)

The number of parameters to be fit is 11

Finally, instead of running the scan we simply load the completed scan performed in Example 8.

n.load_scan()

analysing data from /zfs/nrodd/NPTFRemakeExamples/chains/GCE_Example/.txt

Analysis

An instance of nptf.NPTF with a loaded scan as above can already be used to access the posterior chains with n.samples:

print(np.shape(n.samples))
print(n.samples)

(727, 11)
[[ 3.02530027e-01  1.28844616e-02  7.61010110e-01 ...  3.33197176e+00
  -1.75810216e+00  1.21244569e+01]
 [ 1.42579079e-01  4.47359376e-02  8.78739424e-01 ...  3.23668719e+00
  -1.01507893e-02  2.16130996e+01]
 [ 3.55121759e-01  6.53227405e-02  7.85106865e-01 ...  2.99195590e+00
   7.37752368e-01  2.04954322e+01]
 ...
 [ 2.29948101e-01  1.61386174e-02  8.31147674e-01 ...  2.69403412e+00
  -3.40231551e-01  9.37923721e+00]
 [ 2.12020248e-01  1.13076016e-02  8.34925625e-01 ...  3.54014211e+00
  -1.02525168e+00  2.08287337e+01]
 [ 1.35120585e-01  1.12892027e-03  8.63017349e-01 ...  3.15637545e+00
   5.02206794e-01  2.85446996e+01]]

In the analysis module described next we provide basic helper functions to load in and manipulate these chain samples.

0. Initialize Analysis Module

The first thing to do is initialize an instance of the analysis module, dnds_analysis from dnds_analysis.py with a provided instance of nptfit.NPTF. The NPTF instance should have a scan already loaded in, as done with n.load_scan() above.

an = dnds_analysis.Analysis(n)

dnds_analysis has an optional argument mask, which if unset defaults to the mask in the passed instance of NPTF. If a mask is given, however, then the analysis will be performed in a different ROI to the main run.

1. Make triangle plots

Triangle/corner plots let us visualize multidimensional samples using a scatterplot matrix. A triangle plot with the default options can be made as follows.

an.make_triangle()

To use your own custom plotting options, use corner as follows

corner.corner(an.nptf.samples, labels=an.nptf.params, range=[1 for i in range(an.nptf.n_params)])

with additional arguments as specified in http://corner.readthedocs.io/en/latest/.

2. Get Intensities

Template intensities can be calculated with

dnds_analysis.return_intensity_arrays_poiss(comp)
dnds_analysis.return_intensity_arrays_non_poiss(comp)

for the Poissonian and non-Poissonian templates respectively. This returns an intensity array corresponding to each chain sample associated with the template comp.

The NPT intensity is calculated by integrating up \(\int_{S_{min}}^{S_{max}} dS~S~dN/dS\). This is approximated as a sum between \(S_{min}\) and \(S_{max}\). The options associated with the non-Poissonian template intensity are:

Argument	Default Value	Purpose
comp		The NPT key
smin	0.01	Minimum counts to sum up from
smax	10000	Maximum counts to sum up to
nsteps	10000	Number of bins in s while summing up

We can then look at the quantiles of this distribution, for example to see the middle 68% along with the medians of the GCE and disk NPT as well as that of the GCE PT:

print("GCE NPT Intensity", corner.quantile(an.return_intensity_arrays_non_poiss('gce_np'),[0.16,0.5,0.84]), "ph/cm^2/s")
print("Disk NPT Intensity", corner.quantile(an.return_intensity_arrays_non_poiss('dsk_np'),[0.16,0.5,0.84]), "ph/cm^2/s")
print("GCE PT Intensity", corner.quantile(an.return_intensity_arrays_poiss('gce'),[0.16,0.5,0.84]), "ph/cm^2/s")

GCE NPT Intensity [1.15207873e-07 1.27800928e-07 1.41013411e-07] ph/cm^2/s
Disk NPT Intensity [1.94686971e-08 3.06909244e-08 4.59075661e-08] ph/cm^2/s
GCE PT Intensity [7.26099026e-10 3.05887937e-09 7.86384174e-09] ph/cm^2/s

3. Plot Source Count Distributions

The posterior arrays for the source count distributions \(dN/dF\) [counts:math:^{-1} cm\(^2\) s deg\(^{-2}\)] associated with a given template comp at a given flux [counts/cm:math:^2/s] can be obtained using

dnds.return_dndf_arrays(comp,flux)

The quantiles of this can then be obtained as before. For example, the middle 68% and medians for the GCE and disk non-Poissonian templates:

print(corner.quantile(an.return_dndf_arrays('gce_np',1e-12),[0.16,0.5,0.84]))
print(corner.quantile(an.return_dndf_arrays('dsk_np',1e-12),[0.16,0.5,0.84]))

[3.11754163e+05 6.13090927e+06 2.96631422e+08]
[1.26336916e+04 2.10184377e+06 9.07314393e+08]

The following arrays are used to show the resolved 3FGL points sources and associated Poisson errors as appropriate for the plots below. For how these were obtained, see this snippet.

x_counts, y_counts, error_L, error_H, x_errors_L, x_errors_H = \
[np.array([  1.36887451e-10,   2.56502091e-10,   4.80638086e-10,
          9.00628020e-10,   1.68761248e-09,   3.16227766e-09,
          5.92553098e-09,   1.11033632e-08,   2.08056754e-08,
          3.89860370e-08,   7.30527154e-08]),
 np.array([  1.04000127e+08,   1.83397053e+08,   9.65856820e+07,
          1.51198295e+07,   4.76804443e+06,   9.78677656e+05,
          2.08916332e+05,   0.00000000e+00,   0.00000000e+00,
          0.00000000e+00,   0.00000000e+00]),
 np.array([  2.14237668e+07,   2.08831658e+07,   1.10708578e+07,
          3.18362798e+06,   1.29929969e+06,   4.21069315e+05,
          1.34538182e+05,  -5.57461814e-04,  -2.97500603e-04,
         -1.58767124e-04,  -8.47292389e-05]),
 np.array([  2.63822671e+07,   2.34164673e+07,   1.24232945e+07,
          3.93887993e+06,   1.71404939e+06,   6.58746511e+05,
          2.74201404e+05,   1.02159419e+05,   5.45194091e+04,
          2.90953689e+04,   1.55273233e+04]),
 np.array([  3.68874510e-11,   6.91203483e-11,   1.29518913e-10,
          2.42694796e-10,   4.54765736e-10,   8.52147960e-10,
          1.59676969e-09,   2.99205487e-09,   5.60656455e-09,
          1.05056783e-08,   1.96857231e-08]),
 np.array([  5.04942913e-11,   9.46170829e-11,   1.77295138e-10,
          3.32218719e-10,   6.22517224e-10,   1.16648362e-09,
          2.18577733e-09,   4.09574765e-09,   7.67468330e-09,
          1.43809553e-08,   2.69472846e-08])]

The source count distribution can be plotted with

dnds.plot_source_count_median(comp, smin, smax, nsteps, spow, **kwargs)
dnds.plot_source_count_band(comp, smin, smax, nsteps, spow, qs, **kwargs)

The options being the same as for obtaining the NPT intensity above. Additionally, spow is the power \(n\) in \(F^ndN/dF\) to return while plotting, and qs is an array of quantiles for which to return the dN/dF band. We plot here the median in addition to 68% and 95% confidence intervals.

plt.figure(figsize=[6,5])

an.plot_source_count_median('dsk_np',smin=0.01,smax=1000,nsteps=1000,color='royalblue',spow=2,label='Disk PS')
an.plot_source_count_band('dsk_np',smin=0.01,smax=1000,nsteps=1000,qs=[0.16,0.5,0.84],color='royalblue',alpha=0.15,spow=2)
an.plot_source_count_band('dsk_np',smin=0.01,smax=1000,nsteps=1000,qs=[0.025,0.5,0.975],color='royalblue',alpha=0.1,spow=2)


an.plot_source_count_median('gce_np',smin=0.01,smax=1000,nsteps=1000,color='firebrick',spow=2,label='GCE PS')
an.plot_source_count_band('gce_np',smin=0.01,smax=1000,nsteps=1000,qs=[0.16,0.5,0.84],color='firebrick',alpha=0.15,spow=2)
an.plot_source_count_band('gce_np',smin=0.01,smax=1000,nsteps=1000,qs=[0.025,0.5,0.975],color='firebrick',alpha=0.1,spow=2)

plt.errorbar(x_counts,x_counts**2*y_counts,xerr=[x_errors_L,x_errors_H],yerr=x_counts**2*np.array([error_L,error_H]), fmt='o', color='black', label='3FGL PS')


plt.yscale('log')
plt.xscale('log')
plt.xlim([1e-10,1e-8])
plt.ylim([2e-13,1e-10])

plt.tick_params(axis='x', length=5, width=2, labelsize=18)
plt.tick_params(axis='y', length=5, width=2, labelsize=18)
plt.ylabel('$F^2 dN/dF$ [counts cm$^{-2}$s$^{-1}$deg$^{-2}$]', fontsize=18)
plt.xlabel('$F$  [counts cm$^{-2}$ s$^{-1}$]', fontsize=18)
plt.title(r'Galactic Center NPTF', y=1.02)
plt.legend(fancybox=True, loc='lower right');
plt.tight_layout()


# plt.savefig("dnds_masked.pdf")

As some references also show \(dN/dF\), and we give an example of it below, also demonstrating the use of spow.

plt.figure(figsize=[6,5])

an.plot_source_count_median('dsk_np',smin=0.01,smax=1000,nsteps=1000,color='royalblue',spow=0,label='Disk PS')
an.plot_source_count_band('dsk_np',smin=0.01,smax=1000,nsteps=1000,qs=[0.16,0.5,0.84],color='royalblue',alpha=0.15,spow=0)
an.plot_source_count_band('dsk_np',smin=0.01,smax=1000,nsteps=1000,qs=[0.025,0.5,0.975],color='royalblue',alpha=0.1,spow=0)


an.plot_source_count_median('gce_np',smin=0.01,smax=1000,nsteps=1000,color='firebrick',spow=0,label='GCE PS')
an.plot_source_count_band('gce_np',smin=0.01,smax=1000,nsteps=1000,qs=[0.16,0.5,0.84],color='firebrick',alpha=0.15,spow=0)
an.plot_source_count_band('gce_np',smin=0.01,smax=1000,nsteps=1000,qs=[0.025,0.5,0.975],color='firebrick',alpha=0.1,spow=0)

plt.errorbar(x_counts, y_counts,xerr=[x_errors_L,x_errors_H],yerr=np.array([error_L,error_H]), fmt='o', color='black', label='3FGL PS')


plt.yscale('log')
plt.xscale('log')
plt.xlim([5e-11,5e-9])
plt.ylim([2e5,2e9])
plt.tick_params(axis='x', length=5, width=2, labelsize=18)
plt.tick_params(axis='y', length=5, width=2, labelsize=18)
plt.ylabel('$dN/dF$ [counts$^{-1}$cm$^2$ s deg$^{-2}$]', fontsize=18)
plt.xlabel('$F$  [counts cm$^{-2}$ s$^{-1}$]', fontsize=18)
plt.title('Galactic Center NPTF', y=1.02)
plt.legend(fancybox=True);

4. Plot Intensity Fractions

Intensity fractions (fraction of template intensity to total intensity) for Poissonian and non-Poissonian templates respectively can be plotting using

dnds.plot_intensity_fraction_poiss(comp, bins, **kwargs)
dnds.plot_intensity_fraction_non_poiss(comp, bins, **kwargs)

where comp is the template key, bins is the number of bins between 0 and 100 and **kwargs specify plotting options.

an.plot_intensity_fraction_non_poiss('gce_np', bins=800, color='firebrick', label='GCE PS')
an.plot_intensity_fraction_poiss('gce', bins=800, color='darkgrey', label='GCE DM')
plt.xlabel('Flux fraction (%)')
plt.legend(fancybox = True)
plt.xlim(0,6);
plt.ylim(0,.3);

This plot makes it clear, that when given the choice, the fit prefers to put the GCE flux into point sources rather than diffuse emission.

5. Access Parameter Posteriors

While the posteriors can be accessed with n.samples (or an.nptf.samples) as above, the following functions provide a useful interfact to access individual parameters:

dnds_analysis.return_poiss_parameter_posteriors(comp)
dnds_analysis.return_poiss_parameter_posteriors(comp)

where comp is the (non-)Poissonian template key.

Poissonian parameters

Posterior normalizations of Poissonian parameters can be loaded simply as:

Aiso_poiss_post = an.return_poiss_parameter_posteriors('iso')
Agce_poiss_post = an.return_poiss_parameter_posteriors('gce')
Abub_poiss_post = an.return_poiss_parameter_posteriors('bub')

These can then be use in any way required, for example simply plotted:

f, axarr = plt.subplots(nrows = 1, ncols=3)
f.set_figwidth(12)
f.set_figheight(4)

axarr[0].hist(Aiso_poiss_post, histtype='stepfilled', color='cornflowerblue', bins=np.linspace(0,1.,30), alpha=.4);
axarr[0].set_title('$A_\mathrm{iso}$')
axarr[1].hist(Agce_poiss_post, histtype='stepfilled', color='lightsalmon', bins=np.linspace(0,.2,30), alpha=.4);
axarr[1].set_title('$A_\mathrm{gce}$')
axarr[2].hist(Abub_poiss_post, histtype='stepfilled', color='plum', bins=np.linspace(.5,1.5,30), alpha=.4);
axarr[2].set_title('$A_\mathrm{bub}$')

plt.setp([a.get_yticklabels() for a in axarr], visible=False);

plt.tight_layout()

Non-poissonian parameters

A similar syntax can be used to extract the non-Poissonian parameters.

Agce_non_poiss_post, n_non_poiss_post, Sb_non_poiss_post = an.return_non_poiss_parameter_posteriors('gce_np')

f, axarr = plt.subplots(2, 2);
f.set_figwidth(8)
f.set_figheight(8)


axarr[0, 0].hist(Agce_non_poiss_post, histtype='stepfilled', color='cornflowerblue', bins=np.linspace(0,0.02,30), alpha=.4);
axarr[0, 0].set_title('$A_\mathrm{gce}^\mathrm{ps}$')
axarr[0, 1].hist(n_non_poiss_post[0], histtype='stepfilled', color='lightsalmon', bins=np.linspace(2,30,30), alpha=.4);
axarr[0, 1].set_title('$n_1^\mathrm{gce}$')
axarr[1, 0].hist(n_non_poiss_post[1], histtype='stepfilled', color='lightsalmon', bins=np.linspace(-2,2,30), alpha=.4);
axarr[1, 0].set_title('$n_2^\mathrm{gce}$')
axarr[1, 1].hist(Sb_non_poiss_post, histtype='stepfilled', color='plum', bins=np.linspace(0,40,30), alpha=.4);
axarr[1, 1].set_title('$S_b^{(1), \mathrm{gce}}$')

plt.setp(axarr[0, 0], xticks=[x*0.01 for x in range(5)])
plt.setp(axarr[1, 0], xticks=[x*1.0-2 for x in range(5)])
plt.setp(axarr[1, 1], xticks=[x*10 for x in range(6)])
plt.setp([a.get_yticklabels() for a in axarr[:, 1]], visible=False);
plt.setp([a.get_yticklabels() for a in axarr[:, 0]], visible=False);

plt.tight_layout()

6. Bayesian log-evidence

Finally the Bayesian log-evidence and associated error can be accessed as follows.

l_be, l_be_err = an.get_log_evidence()
print(l_be, l_be_err)

-30956.359589319734 0.43278107806496635

Example 10: Analyzing the Results of the High-Lat Run

In this Example we analyze the results of an MPI run of NPTFit performed over high latitudes.

The example batch script we provide, Example10_HighLat_Batch.batch is for SLURM. This calls the run file Example10_HighLat_Run.py using MPI, and is an example of how to perform a more realistic analysis using NPTFit.

NB: The batch file must be run before this notebook.

NB: This example makes use of the Fermi Data, which needs to already be installed. See Example 1 for details.

In this example, we model the source-count function as a triply-broken power law. In detail, the source count function is then:

\[\begin{split}\frac{dN}{dS} = A \left\{ \begin{array}{lc} \left( \frac{S}{S_{b,1}} \right)^{-n_1}, & S \geq S_{b,1} \\ \left(\frac{S}{S_{b,1}}\right)^{-n_2}, & S_{b,1} > S \geq S_{b,2} \\ \left( \frac{S_{b,2}}{S_{b,1}} \right)^{-n_2} \left(\frac{S}{S_{b,2}}\right)^{-n_3}, & S_{b,2} > S \geq S_{b,3} \\ \left( \frac{S_{b,2}}{S_{b,1}} \right)^{-n_2} \left( \frac{S_{b,3}}{S_{b,2}} \right)^{-n_3} \left(\frac{S}{S_{b,3}}\right)^{-n_4}, & S_{b,3} > S \end{array} \right.\end{split}\]

and is thereby described by the following eight parameters:

\[\theta = \left[ A, n_1, n_2, n_3, n_4, S_b^{(1)}, S_b^{(2)}, S_b^{(3)} \right]\,.\]

This provides an example of a more complicated source count function, and also explains why the run needs MPI.

# Import relevant modules

%matplotlib inline
%load_ext autoreload
%autoreload 2

import numpy as np
import corner
import matplotlib.pyplot as plt

from NPTFit import nptfit # module for performing scan
from NPTFit import create_mask as cm # module for creating the mask
from NPTFit import dnds_analysis # module for analysing the output
from NPTFit import psf_correction as pc # module for determining the PSF correction

from __future__ import print_function

Load in scan

We need to create an instance of nptfit.NPTF and load in the scan performed using MPI.

n = nptfit.NPTF(tag='HighLat_Example')

fermi_data = np.load('fermi_data/fermidata_counts.npy').astype(np.int32)
fermi_exposure = np.load('fermi_data/fermidata_exposure.npy')
n.load_data(fermi_data, fermi_exposure)

analysis_mask = cm.make_mask_total(band_mask = True, band_mask_range = 50)
n.load_mask(analysis_mask)

dif = np.load('fermi_data/template_dif.npy')
iso = np.load('fermi_data/template_iso.npy')

n.add_template(dif, 'dif')
n.add_template(iso, 'iso')
n.add_template(np.ones(len(iso)), 'iso_np', units='PS')

n.add_poiss_model('dif','$A_\mathrm{dif}$', [0,20], False)
n.add_poiss_model('iso','$A_\mathrm{iso}$', [0,5], False)

n.add_non_poiss_model('iso_np',
                      ['$A^\mathrm{ps}_\mathrm{iso}$',
                      '$n_1$','$n_2$','$n_3$','$n_4$',
                      '$S_b^{(1)}$','$S_b^{(2)}$','$S_b^{(3)}$'],
                      [[-6,2],
                      [2.05,5],[1.0,3.5],[1.0,3.5],[-1.99,1.99],
                      [30,80],[1,30],[0.1,1]],
                      [True,False,False,False,False,False,False,False])

pc_inst = pc.PSFCorrection(psf_sigma_deg=0.1812)
f_ary, df_rho_div_f_ary = pc_inst.f_ary, pc_inst.df_rho_div_f_ary

Loading the psf correction from: /zfs/nrodd/NPTFRemakeExamples/psf_dir/gauss_128_0.181_10_50000_1000_0.01.npy

n.configure_for_scan(f_ary=f_ary, df_rho_div_f_ary=df_rho_div_f_ary, nexp=5)

The number of parameters to be fit is 10

Finally, load the completed scan performed using MPI.

n.load_scan()

analysing data from /zfs/nrodd/NPTFRemakeExamples/chains/HighLat_Example/.txt

Analysis

As in Example 9 we first initialize the analysis module. We will provide the same basic plots as in that notebook, where more details on each option is provided.

an = dnds_analysis.Analysis(n)

1. Make triangle plots

an.make_triangle()

2. Get Intensities

print("Iso NPT Intensity",corner.quantile(an.return_intensity_arrays_non_poiss('iso_np'),[0.16,0.5,0.84]), "ph/cm^2/s")
print("Iso PT Intensity",corner.quantile(an.return_intensity_arrays_poiss('iso'),[0.16,0.5,0.84]), "ph/cm^2/s")
print("Dif PT Intensity",corner.quantile(an.return_intensity_arrays_poiss('dif'),[0.16,0.5,0.84]), "ph/cm^2/s")

Iso NPT Intensity [1.13023625e-07 1.21451391e-07 1.29983760e-07] ph/cm^2/s
Iso PT Intensity [1.62091903e-07 1.68566455e-07 1.74563161e-07] ph/cm^2/s
Dif PT Intensity [1.88574352e-07 1.93195545e-07 1.97764592e-07] ph/cm^2/s

3. Plot Source Count Distributions

an.plot_source_count_median('iso_np',smin=0.01,smax=1000000,nsteps=10000,color='tomato',spow=2)
an.plot_source_count_band('iso_np',smin=0.01,smax=1000000,nsteps=10000,qs=[0.16,0.5,0.84],color='tomato',alpha=0.3,spow=2)

plt.yscale('log')
plt.xscale('log')
plt.xlim([1e-12,5e-6])
plt.ylim([5e-14,1e-11])
plt.tick_params(axis='x', length=5, width=2, labelsize=18)
plt.tick_params(axis='y', length=5, width=2, labelsize=18)
plt.ylabel('$F^2 dN/dF$ [counts cm$^{-2}$s$^{-1}$deg$^{-2}$]', fontsize=18)
plt.xlabel('$F$  [counts cm$^{-2}$ s$^{-1}$]', fontsize=18)
plt.title('High Latitudes Isotropic NPTF', y=1.02)

Text(0.5,1.02,'High Latitudes Isotropic NPTF')

4. Plot Intensity Fractions

an.plot_intensity_fraction_non_poiss('iso_np', bins=100, color='tomato', label='Iso PS')
an.plot_intensity_fraction_poiss('iso', bins=100, color='cornflowerblue', label='Iso')
an.plot_intensity_fraction_poiss('dif', bins=100, color='plum', label='Dif')
plt.xlabel('Flux fraction (%)')
plt.legend(fancybox = True)
plt.xlim(0,80);

5. Access Parameter Posteriors

Poissonian parameters

Aiso_poiss_post = an.return_poiss_parameter_posteriors('iso')
Adif_poiss_post = an.return_poiss_parameter_posteriors('dif')

f, axarr = plt.subplots(1, 2);
f.set_figwidth(8)
f.set_figheight(4)


axarr[0].hist(Aiso_poiss_post, histtype='stepfilled', color='cornflowerblue', bins=np.linspace(.5,1,30),alpha=0.4);
axarr[0].set_title('$A_\mathrm{iso}$')
axarr[1].hist(Adif_poiss_post, histtype='stepfilled', color='lightsalmon', bins=np.linspace(10,15,30),alpha=0.4);
axarr[1].set_title('$A_\mathrm{dif}$')

plt.setp([a.get_yticklabels() for a in axarr[:]], visible=False);

plt.tight_layout()

Non-poissonian parameters

Aiso_non_poiss_post, n_non_poiss_post, Sb_non_poiss_post = an.return_non_poiss_parameter_posteriors('iso_np')

f, axarr = plt.subplots(2, 4);
f.set_figwidth(16)
f.set_figheight(8)

axarr[0, 0].hist(Aiso_non_poiss_post, histtype='stepfilled', color='cornflowerblue', bins=np.linspace(0,.0001,30),alpha=0.4);
axarr[0, 0].set_title('$A_\mathrm{iso}^\mathrm{ps}$')
axarr[0, 1].hist(n_non_poiss_post[0], histtype='stepfilled', color='lightsalmon', bins=np.linspace(2,4,30),alpha=0.4);
axarr[0, 1].set_title('$n_1^\mathrm{iso}$')
axarr[0, 2].hist(n_non_poiss_post[1], histtype='stepfilled', color='lightsalmon', bins=np.linspace(1,3.5,30),alpha=0.4);
axarr[0, 2].set_title('$n_2^\mathrm{iso}$')
axarr[0, 3].hist(n_non_poiss_post[2], histtype='stepfilled', color='lightsalmon', bins=np.linspace(1,3.5,30),alpha=0.4);
axarr[0, 3].set_title('$n_3^\mathrm{iso}$')
axarr[1, 0].hist(n_non_poiss_post[3], histtype='stepfilled', color='lightsalmon', bins=np.linspace(-2,2,30),alpha=0.4);
axarr[1, 0].set_title('$n_4^\mathrm{iso}$')
axarr[1, 1].hist(Sb_non_poiss_post[0], histtype='stepfilled', color='plum', bins=np.linspace(30,80,30),alpha=0.4);
axarr[1, 1].set_title('$S_b^{(1), \mathrm{iso}}$')
axarr[1, 2].hist(Sb_non_poiss_post[1], histtype='stepfilled', color='plum', bins=np.linspace(1,30,30),alpha=0.4);
axarr[1, 2].set_title('$S_b^{(2), \mathrm{iso}}$')
axarr[1, 3].hist(Sb_non_poiss_post[2], histtype='stepfilled', color='plum', bins=np.linspace(0.1,1,30),alpha=0.4);
axarr[1, 3].set_title('$S_b^{(3), \mathrm{iso}}$')

plt.setp(axarr[0, 0], xticks=[x*.00005 for x in range(3)])
plt.setp(axarr[1, 0], xticks=[x*1.-2.0 for x in range(4)])
plt.setp(axarr[1, 3], xticks=[x*0.2+0.2 for x in range(5)])
plt.setp([a.get_yticklabels() for a in axarr[:, 0]], visible=False);
plt.setp([a.get_yticklabels() for a in axarr[:, 1]], visible=False);
plt.setp([a.get_yticklabels() for a in axarr[:, 2]], visible=False);
plt.setp([a.get_yticklabels() for a in axarr[:, 3]], visible=False);
plt.tight_layout()

Tip

When using NPTFit we recommend the use of National Pipe Thread (NPT) rather than the more common British Standard Pipe (BSP). For details, see NPT.