pFGL.py

import os
import urllib.request

import numpy as np
import pandas as pd
import healpy as hp
import numpy as np
from scipy.integrate import quad
from scipy.stats import chi2
from astropy.table import Table
from tqdm import *


class plot_FGL():

    def __init__(self, catalog='3FGL'):
        
        self.catalog = catalog
        self.data_dir = os.getcwd()

        self.load_data()
        self.errors_sigma = [poisson_interval(i) for i in range(3500)]

    def load_data(self):
        if self.catalog == '3FGL':
            filename = 'gll_psc_v16.fit'
        elif self.catalog is '4FGL':
            filename = 'gll_psc_v18.fit'

        if not os.path.exists(self.data_dir + filename):
            if self.catalog == '3FGL':
                url = 'https://fermi.gsfc.nasa.gov/ssc/data/access/lat/4yr_catalog/gll_psc_v16.fit'  
                urllib.request.urlretrieve(url, self.data_dir + 'gll_psc_v16.fit')
            else:
                url = 'https://fermi.gsfc.nasa.gov/ssc/data/access/lat/8yr_catalog/gll_psc_v18.fit'  
                urllib.request.urlretrieve(url, self.data_dir + 'gll_psc_v18.fit')

        self.cat = Table.read(self.data_dir + filename, hdu=1)

    def lb2pix(self, nside, l, b, nest=False):
        """ Convert right ascension and descent to galactic
            coordinates, then get index corresponding to HEALPix array
        """
        return hp.ang2pix(nside, np.deg2rad(90 - b), np.deg2rad(l), nest=nest)

    def return_counts(self, emin=2., emax=20., flux_min=7e-12, flux_max=4e-10, flux_bins=14, nside=128, mask=None):
        """ Return histogrammed source counts from 3FGL data
        """
        self.nside = hp.npix2nside(len(mask))
        self.mask_total = mask
        self.energy_range = [emin, emax] 

        self.area_factor = 1.

        # Reduce dataframe to unmasked region
        if self.mask_total is not None:  
            to_include = []
            for i in range(len(self.cat)):
                source_pix = (self.lb2pix(self.nside, self.cat[
                              'GLON'][i], self.cat['GLAT'][i]))
                if self.mask_total[source_pix] == 0:
                    to_include.append(True)
                else:
                    to_include.append(False)
            self.cat_reduced = self.cat[to_include]
            self.area_factor = 1. - \
                np.sum(self.mask_total) / float(hp.nside2npix(self.nside))
        else:
            self.cat_reduced = self.cat

        self.fluxes_3fgl = []

        for src in tqdm_notebook(self.cat_reduced):

            spectrum = src['SpectrumType'].strip()
            
            if self.catalog == '3FGL':
                if  spectrum == 'PowerLaw':
                    flux = [quad(lambda E: y_powerlaw(E, .001 * src['Pivot_Energy'], 1000 * src['Flux_Density'], src['Spectral_Index']),
                                self.energy_range[i], self.energy_range[i + 1])[0] for i in range(len(self.energy_range) - 1)]
                elif spectrum == 'LogParabola':
                    flux = [quad(lambda E: y_logparabola(E, .001 * src['Pivot_Energy'], 1000 * src['Flux_Density'], src['Spectral_Index'],
                                                        src['beta']), self.energy_range[i], self.energy_range[i + 1])[0] for i in range(len(self.energy_range) - 1)]
                elif spectrum in ['PLExpCutoff', 'PLSuperExpCutoff','PLSuperExpCutoff2']:
                    flux = [quad(lambda E: y_expcutoff(E, .001 * src['Pivot_Energy'], 1000 * src['Flux_Density'], src['Spectral_Index'],
                                                    .001 * src['Cutoff'], src['Exp_Index']), self.energy_range[i], self.energy_range[i + 1])[0] for i in range(len(self.energy_range) - 1)]
                else: 
                    raise NotImplementedError

            elif self.catalog == '4FGL':

                if  spectrum == 'PowerLaw':
                    flux = [quad(lambda E: y_powerlaw(E, .001 * src['Pivot_Energy'], 1000 * src['PL_Flux_Density'], src['PL_Index']),
                                self.energy_range[i], self.energy_range[i + 1])[0] for i in range(len(self.energy_range) - 1)]
                elif spectrum == 'LogParabola':
                    flux = [quad(lambda E: y_logparabola(E, .001 * src['Pivot_Energy'], 1000 * src['LP_Flux_Density'], src['LP_Index'],
                                                        src['LP_beta']), self.energy_range[i], self.energy_range[i + 1])[0] for i in range(len(self.energy_range) - 1)]
                elif spectrum in ['PLExpCutoff', 'PLSuperExpCutoff','PLSuperExpCutoff2']:
                    flux = [quad(lambda E: y_expcutoff_4fgl(E, .001 * src['Pivot_Energy'], 1000 * src['PLEC_Flux_Density'], src['PLEC_Index'],
                                                    .001 * src['PLEC_Expfactor'], src['PLEC_Exp_Index']), self.energy_range[i], self.energy_range[i + 1])[0] for i in range(len(self.energy_range) - 1)]
                else: 
                    raise NotImplementedError


            self.fluxes_3fgl.append(flux)

        flux_values_reduced = self.fluxes_3fgl

        deg = 180 / np.pi
        sr = 4 * np.pi
        srdeg2 = sr * deg**2  # (180/np.pi)**2

        counts, bin_edges = np.histogram(flux_values_reduced, bins=np.logspace(
            int(np.log10(flux_min)), int(np.log10(flux_max)), flux_bins))

        bin_centres = 10**((np.log10(bin_edges[:-1]) +
                            np.log10(bin_edges[1:])) / 2.)
        bin_centers = bin_centres  # British to American
        bin_width = bin_edges[1:] - bin_edges[:-1]

        x_counts = bin_centres
        y_counts = np.array(counts / (self.area_factor * bin_width * srdeg2))

        errors = [self.errors_sigma[count] for count in counts]
        error_L = []
        error_H = []
        for i in range(len(counts)):
            error_L.append(counts[i] - errors[i][0] - 10**-8)
            error_H.append(errors[i][1] - counts[i])

        self.error_L = np.array(error_L) / \
            (self.area_factor * bin_width * srdeg2)
        self.error_H = np.array(error_H) / \
            (self.area_factor * bin_width * srdeg2)

        self.x_errors_L = np.array(
            [bin_centers[i] - bin_edges[i] for i in range(np.size(bin_centers))])
        self.x_errors_H = np.array(
            [bin_edges[i + 1] - bin_centers[i] for i in range(np.size(bin_centers))])

        return x_counts, y_counts, self.error_L, self.error_H, self.x_errors_L, self.x_errors_H

## 3FGL/4FGL fit functions

# From https://arxiv.org/pdf/1501.02003.pdf (3FGL) and https://arxiv.org/pdf/1902.10045.pdf (4FGL)

def y_powerlaw(E, E0, K, Gamma):
    return K * (E / E0)**(-Gamma)

def y_logparabola(E, E0, K, Gamma, beta):
    return K * (E / E0)**(-Gamma - beta * np.log(E / E0))

def y_expcutoff(E, E0, K, Gamma, Ec, b):
    return K * (E / E0)**(-Gamma) * np.exp((E0 / Ec)**b - (E / Ec)**b)

def y_expcutoff_4fgl(E, E0, K, Gamma, a, b):
    return K * (E / E0)**(-Gamma) * np.exp(a * ((E0)**b - (E)**b))

def poisson_interval(k, alpha=0.32): 
    """ Uses chisquared info to get the poisson interval.poisson
        Stolen from http://stackoverflow.com/questions/14813530/poisson-confidence-interval-with-numpy
    """
    a = alpha
    low, high = (chi2.ppf(a/2, 2*k) / 2, chi2.ppf(1-a/2, 2*k + 2) / 2)
    if k == 0: 
        low = 0.0
    return low, high