p3FGL.py

###############################################################################
# p3FGL.py
###############################################################################
#
# Plot Fermi 3FGL PS catalog histogram
# Usage: 
#
# p = plot_3FGL()
# x_counts, y_counts, error_L, error_H, x_errors_L, x_errors_H = \
#                 p.return_counts(flux_min = 5e-11,
#                                 flux_max = 1e-7,
#                                 flux_bins = 12)
# # Plot F^2*dN/dF:
# 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')
#
###############################################################################

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

class plot_3FGL():

    def __init__(self):

        self.load_3FGL_data()
        self.errors_sigma = [poisson_interval(i) for i in range(1500)]

    def load_3FGL_data(self):

        # Load in the 3FGL catalog
        self.source_3fg_df = pd.read_csv('3fgl.dat', sep='|', comment='#')

        # Remove whitespace from column names
        self.source_3fg_df.rename(columns=lambda x: x.strip(), inplace=True)
        for col in self.source_3fg_df.columns.values:
            try:
                self.source_3fg_df[col] = self.source_3fg_df[
                    col].map(str.strip)
            except TypeError:
                continue

        # Convert to numeric data
        self.source_3fg_df = self.source_3fg_df.convert_objects(
            convert_numeric=True)

    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 = nside
        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.source_3fg_df)):
                source_pix = (self.lb2pix(self.nside, self.source_3fg_df[
                              '_Lii'].values[i], self.source_3fg_df['_Bii'].values[i]))
                if self.mask_total[source_pix] == 0:
                    to_include.append(i)
            self.source_3fg_df = self.source_3fg_df.iloc[to_include]
            self.area_factor = 1. - \
                np.sum(self.mask_total) / float(hp.nside2npix(self.nside))

        self.fluxes_3fgl = []

        for index, row in self.source_3fg_df.iterrows():
            if row['spectrum_type'] == 'PowerLaw':
                flux = [quad(lambda E: y_powerlaw(E, .001 * row['pivot_energy'], 1000 * row['flux_density'], row['spectral_index']),
                             self.energy_range[i], self.energy_range[i + 1])[0] for i in range(len(self.energy_range) - 1)]
            elif row['spectrum_type'] == 'LogParabola':
                flux = [quad(lambda E: y_logparabola(E, .001 * row['pivot_energy'], 1000 * row['flux_density'], row['spectral_index'],
                                                     row['beta']), self.energy_range[i], self.energy_range[i + 1])[0] for i in range(len(self.energy_range) - 1)]
            elif row['spectrum_type'] == 'PLExpCutoff':
                flux = [quad(lambda E: y_expcutoff(E, .001 * row['pivot_energy'], 1000 * row['flux_density'], row['spectral_index'],
                                                   .001 * row['cutoff'], row['exp_index']), self.energy_range[i], self.energy_range[i + 1])[0] for i in range(len(self.energy_range) - 1)]
            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))
        y_counts_err = np.array(
            np.sqrt(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 fit functions
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 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