All1FHL2FHLPythonCodes < Veritas

Veritas Web>All1FHL2FHLPythonCodes (2016-10-23, RichardNederlander) (raw view)
---+ 1FHL and 2FHL Python Codes: Created by Richard Nederlander, CC '18

Description: Cumulative Plot of 1FHL and 2FHL data

Created by Richard Nederlander, Columbia University

Last Use: 10/23/16

In [7]:
<pre>%matplotlib inline
import matplotlib
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import pyfits
import numpy as np

hdulist = pyfits.open('gll_psch_v07.fit')

cols = hdulist[1].columns

#print(cols.info())
#print(cols.names[0])
#print(cols.names[1])
#print(cols.names[2])
#print(cols.names[3])
#print(cols.names[4])
#print(cols.names[16])
#print(cols.names[18])
#print(cols.names[22])

tbdata = hdulist[1].data
name = tbdata.field(0)
raj = tbdata.field(1)
dej = tbdata.field(2)
glon = tbdata.field(3)
glat = tbdata.field(4)
flux = tbdata.field(16)
energy = tbdata.field(18)
spectral = tbdata.field(22)

#for 1FHL
#print(name[116])

#for 2FHL
#print(name[85])

hdulist.close()
</pre>

In [8]:
<pre>fig, ax = plt.subplots()
fig.set_size_inches(w=6, h=6)

#for 1FHL
crab = flux[116]

#for 2FHL
#crab = flux[85]

x=[]
glon1=[]

print(max(flux))
print(min(flux))

for i in range(0, len(glon)):
    if glat[i] &gt;=-5 and glat[i] &lt;= 5:
        x.append(flux[i]/crab)
    else:
        pass

n, bins, patches = plt.hist(x, bins = np.logspace(-2, 0, 23), normed=False, histtype='step', cumulative=-1, color='black', orientation='vertical')

#change first number in logspace with respect to last number in flux list 

#print(max(x))
#print(min(x))
#print(x[-1])

print(len(bins[:-1]))
print(len(n))
print(len(x))

plt.xlabel('Flux (Crab units)', fontsize=15)
plt.ylabel('Number of sources (F&gt;Flux)', fontsize=15)
plt.title('Cumulative Galactic source count of 1FHL Sources as a function of VHE gamma-ray flux', position=(0.5,1.05), fontsize=15)
plt.xscale('log')
plt.yscale('log')
plt.xlim(0.001, 1.000)
plt.ylim(1, 1000)

plt.grid(False)

plt.scatter(bins[:-1], n, alpha=1, color="black")

ax.axvspan(0.0018, 0.0038, color='blue', alpha=0.3, lw=0)
plt.axvline(x=0.085, ymax=100, color='black', linewidth=1, linestyle='--')
plt.text(0.58, 0.18, 'H.E.S.S. survey', horizontalalignment='center', verticalalignment='center',
        rotation='vertical', fontsize=10, transform=ax.transAxes, bbox=dict(facecolor='none', edgecolor='none', alpha=0.3))
plt.text(0.63, 0.18, 'completeness', horizontalalignment='right', verticalalignment='center', 
         rotation='vertical', fontsize=10, transform=ax.transAxes, bbox=dict(facecolor='none', edgecolor='none', alpha=0.3))
</pre> <pre>3.57666e-10
1.06254e-12
22
22
76
</pre>

Out[8]:
<pre>&lt;matplotlib.text.Text at 0x10726fa90&gt;</pre>

In [9]:
<pre>a=0.14
#a=0.35

ax.text(a, 0.80, "T", horizontalalignment='center', verticalalignment='center',
        fontsize=10, transform=ax.transAxes, bbox=dict(facecolor='none', edgecolor='none', alpha=0.3))
ax.text(a, 0.74, "A", horizontalalignment='center', verticalalignment='center',
        fontsize=10, transform=ax.transAxes, bbox=dict(facecolor='none', edgecolor='none', alpha=0.3))
ax.text(a, 0.68, "R", horizontalalignment='center', verticalalignment='center',
        fontsize=10, transform=ax.transAxes, bbox=dict(facecolor='none', edgecolor='none', alpha=0.3))
ax.text(a, 0.62, "G", horizontalalignment='center', verticalalignment='center',
        fontsize=10, transform=ax.transAxes, bbox=dict(facecolor='none', edgecolor='none', alpha=0.3))
ax.text(a, 0.56, "E", horizontalalignment='center', verticalalignment='center',
        fontsize=10, transform=ax.transAxes, bbox=dict(facecolor='none', edgecolor='none', alpha=0.3))
ax.text(a, 0.50, "T", horizontalalignment='center', verticalalignment='center',
        fontsize=10, transform=ax.transAxes, bbox=dict(facecolor='none', edgecolor='none', alpha=0.3))
ax.text(a, 0.44, " ", horizontalalignment='center', verticalalignment='center',
        fontsize=10, transform=ax.transAxes, bbox=dict(facecolor='none', edgecolor='none', alpha=0.3))
ax.text(a, 0.38, "R", horizontalalignment='center', verticalalignment='center',
        fontsize=10, transform=ax.transAxes, bbox=dict(facecolor='none', edgecolor='none', alpha=0.3))
ax.text(a, 0.32, "A", horizontalalignment='center', verticalalignment='center',
        fontsize=10, transform=ax.transAxes, bbox=dict(facecolor='none', edgecolor='none', alpha=0.3))
ax.text(a, 0.26, "N", horizontalalignment='center', verticalalignment='center',
        fontsize=10, transform=ax.transAxes, bbox=dict(facecolor='none', edgecolor='none', alpha=0.3))
ax.text(a, 0.20, "G", horizontalalignment='center', verticalalignment='center',
        fontsize=10, transform=ax.transAxes, bbox=dict(facecolor='none', edgecolor='none', alpha=0.3))
ax.text(a, 0.14, "E", horizontalalignment='center', verticalalignment='center',
        fontsize=10, transform=ax.transAxes, bbox=dict(facecolor='none', edgecolor='none', alpha=0.3))
</pre>

Out[9]:
<pre>&lt;matplotlib.text.Text at 0x10733a950&gt;</pre>

In [10]:
<pre>fig.savefig("cumul1FHL.png", bbox_inches='tight')
#fig.show()
print("Ready!")
</pre> <pre>Ready!
</pre>

Description: Histogram Plotter of 1FHL and 2FHL data

Created by Richard Nederlander, Columbia University

Last Use: 10/23/16

In [ ]:
<pre>import pyfits
import numpy as np
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt

hdulist = pyfits.open('gll_psch_v07.fit')
cols = hdulist[1].columns

#print(cols.info())
#print(cols.names[0])
#print(cols.names[1])
#print(cols.names[2])
#print(cols.names[3])
#print(cols.names[4])
#print(cols.names[16])
#print(cols.names[22])

tbdata = hdulist[1].data
name = tbdata.field(0)
raj = tbdata.field(1)
dej = tbdata.field(2)
glon = tbdata.field(3)
glat = tbdata.field(4)
flux = tbdata.field(16)
spectral = tbdata.field(22)

hdulist.close()
</pre>

Description: Latitude vs. Longitude Plotter of 1FHL and 2FHL data

Created by Richard Nederlander, Columbia University

Last Use: 10/23/16

In [ ]:
<pre>glon1=[]

for i in range(0, len(glon)):
    if glat[i] &gt;=-5 and glat[i] &lt;= 5 and glon[i] &gt;= 0 and glon[i] &lt;= 180:
        glon1.append(glon[i])
    elif glat[i] &gt;=-5 and glat[i] &lt;= 5 and glon[i] &gt; 180:
        glon1.append(glon[i]-360)
    else:
        pass
    
plt.hist(glon1,bins=17,alpha=1,color="green", align="mid")
plt.xlim([-200, 200])
plt.ylim([0, 12])
plt.grid(False)

plt.xlabel('Galactic longitude (deg)', fontsize=15)
plt.ylabel('Number of sources', fontsize=15)
plt.title('Histogram of Fermi-LAT 1FHL sources within 5&#x2DA; of the Galactic plane', y=1.03, fontsize=15)

print(len(glon))
print(len(glon1))

plt.savefig("hist1FHL.png", bbox_inches='tight')
#plt.show()
print("Ready!")<br /><br /></pre> <pre>%matplotlib inline
import pyfits
import matplotlib.pyplot as plt
import matplotlib.lines as mlines
from matplotlib.font_manager import FontProperties
import numpy as np

hdulist = pyfits.open('gll_psch_v07.fit')
FHL = 1

cols = hdulist[1].columns

#print(cols.info())
#print(cols.names[0])
#print(cols.names[1])
#print(cols.names[2])
#print(cols.names[3])
#print(cols.names[4])
#print(cols.names[16])
#print(cols.names[22])

tbdata = hdulist[1].data
name = tbdata.field(0)
raj = tbdata.field(1)
dej = tbdata.field(2)
glon = tbdata.field(3)
glat = tbdata.field(4)
flux = tbdata.field(16)
spectral = tbdata.field(22)
assoc = tbdata.field(32)
fhl = tbdata.field(34)

hdulist.close()
</pre>

In [2]:
<pre>marker1 = []
label1 = []
blazar = []
galactic = []
unassociated = []

aa = []
bb = []
cc = []
dd = []
ee = []
ff = []
gg = []
hh = []
ii = []
jj = []
kk = []
ll = []
mm = []
nn = []
oo = []
pp = []
qq = []
rr = []
ss = []


if FHL == 1:

    for i in range(len(glon)):
        if assoc[i].lower() == "bzb":
            label1.append("BL Lac")
            marker1.append("^")
            blazar.append(i)
            aa.append(i)
        elif assoc[i].lower() == "bzq":
            label1.append("FSRQ")
            marker1.append("^")
            blazar.append(i)
            bb.append(i)
        elif assoc[i].lower() == "agu":
            label1.append("Active Galaxy of uncertain type")
            marker1.append("^")
            blazar.append(i)
            cc.append(i)
        elif assoc[i].lower() == "hpsr":
            label1.append("Pulsar, identified by pulsations above 10 GeV")
            marker1.append("*")
            galactic.append(i)
            dd.append(i)
        elif assoc[i].lower() == "psr":
            label1.append("Pulsar, no pulsations seen in LAT yet")
            marker1.append("*")
            galactic.append(i)
            ee.append(i)
        elif assoc[i].lower() == "snr":
            label1.append("Supernova remnant")
            marker1.append("o")
            galactic.append(i)
            ff.append(i)
        elif assoc[i].lower() == "pwn":
            label1.append("Pulsar wind nebula")
            marker1.append("D")
            galactic.append(i)
            gg.append(i)
        elif assoc[i].lower() == "spp":
            label1.append("Unclear whether SNR or PWN")
            marker1.append("&lt;")
            galactic.append(i)
            hh.append(i)
        elif assoc[i].lower() == "hmb":
            label1.append("High-mass binary")
            marker1.append("d")
            galactic.append(i)
            ii.append(i)
        elif assoc[i].lower() == "gal":
            label1.append("Normal galaxy")
            marker1.append("v")
            galactic.append(i)
            jj.append(i)
        elif assoc[i].lower() == "sfr":
            label1.append("Star forming region")
            marker1.append("h")
            galactic.append(i)
            kk.append(i)
        elif assoc[i].lower() == "lvb":
            label1.append("LBV star")
            marker1.append("h")
            galactic.append(i)
            ll.append(i)
        else:
            label1.append("Unassociated source")
            marker1.append("s")
            unassociated.append(i)
            mm.append(i)

elif FHL == 2:
    
    for i in range(len(glon)):
        if fhl[i].lower() == "psr":
            label1.append("Pulsar")
            marker1.append("*")
            galactic.append(i)
            aa.append(i)
        elif fhl[i].lower() == "pwn":
            label1.append("Pulsar wind nebula")
            marker1.append("D")
            galactic.append(i)
            bb.append(i)
        elif fhl[i].lower() == "snr":
            label1.append("Supernova remnant")
            marker1.append("o")
            galactic.append(i)
            cc.append(i)
        elif fhl[i].lower() == "spp":
            label1.append("Supernova remnant / Pulsar wind nebula")
            marker1.append("o")
            galactic.append(i)
            dd.append(i)
        elif fhl[i].lower() == "hmb":
            label1.append("High-mass binary")
            marker1.append("d")
            galactic.append(i)
            ee.append(i)
        elif fhl[i].lower() == "bin":
            label1.append("Binary")
            marker1.append("d")
            galactic.append(i)
            ff.append(i)
        elif fhl[i].lower() == "sfr":
            label1.append("Star-forming region")
            marker1.append("h") 
            galactic.append(i)
            gg.append(i)
        elif fhl[i].lower() == "bll":
            label1.append("BL Lac type of blazar")
            marker1.append("^")
            blazar.append(i)
            hh.append(i)
        elif fhl[i].lower() == "bll-g":
            label1.append("BL Lac type of blazar with prominent galaxy emission")
            marker1.append("^") 
            blazar.append(i)
            ii.append(i)
        elif fhl[i].lower() == "fsrq":
            label1.append("FSRQ type of blazar")
            marker1.append("^")
            blazar.append(i)
            jj.append(i)
        elif fhl[i].lower() == "agn":
            label1.append("Non-blazar active galaxy")
            marker1.append("v")
            galactic.append(i)
            kk.append(i)
        elif fhl[i].lower() == "rdg":
            label1.append("Radio galaxy")
            marker1.append("^")
            blazar.append(i)
            ll.append(i)
        elif fhl[i].lower() == "rdg/bll":
            label1.append("Radio galaxy / BL Lac")
            marker1.append("^")
            blazar.append(i)
            mm.append(i)
        elif fhl[i].lower() == "bcu I":
            label1.append("Blazar candidate of uncertain type I")
            marker1.append("^") 
            blazar.append(i)
            nn.append(i)
        elif fhl[i].lower() == "bcu II":
            label1.append("Blazar candidate of uncertain type II")
            marker1.append("^") 
            blazar.append(i)
            oo.append(i)
        elif fhl[i].lower() == "bcu III":
            label1.append("Blazar candidate of uncertain type III")
            marker1.append("^")
            blazar.append(i)
            pp.append(i)
        elif fhl[i].lower() == "gal":
            label1.append("Normal galaxy (or part)")
            marker1.append("v")
            blazar.append(i)
            qq.append(i)
        elif fhl[i].lower() == "galclu":
            label1.append("Galaxy cluster")
            marker1.append("v")
            blazar.append(i)
            rr.append(i)
        else:
            label1.append("Unassociated source")
            marker1.append("s")
            unassociated.append(i)
            ss.append(i)

else:
    
    None
</pre>

In [3]:
<pre>#glon1=[]
#glat1=[]
#
#for i in range(0, len(glon)):
#    if glat[i] &gt;=-5 and glat[i] &lt;= 5:
#        glon1.append(glon[i])
#        glat1.append(glat[i])
#    else:
#        pass

x = glon
y = glat

fig, (ax1, ax2, ax3) = plt.subplots(3)

fig.set_size_inches(w=12, h=4)

ax1 = plt.subplot(311)
ax1.axis([180, 60, -5, 5])
ax1.set_yticks([-5,0,5])
ax1.grid(True)

for i in range(0, len(glon)):
    ax1.scatter(x[i], y[i], marker=marker1[i], color="black",
                edgecolor='black', facecolor='None', s=70, linewidth='1.5')

ax1.set_title('SkyMap of Fermi-LAT 1FHL sources within 5&#x2DA; of the Galactic plane', y=1.08, fontsize='15')

ax2 = plt.subplot(312)

glon2=[]
glat2=[]

for i in range(0, len(glon)):
    if glon[i] &gt;= 60:
        glon2.append(glon[i]-360)
        glat2.append(glat[i])
    else:
        glon2.append(glon[i])
        glat2.append(glat[i])
        
x = glon2
y = glat2
ax2.set_xlim([60, -60])
ax2.set_ylim([-5, 5])
ax2.set_yticks([-5,0,5])
ax2.grid(True)

for i in range(0, len(glon)):
    ax2.scatter(x[i], y[i], marker=marker1[i], color="black", 
                edgecolor='black', facecolor='None', s=70, linewidth='1.5')
    
ax2.set_ylabel('Galactic latitude (deg)', fontsize='15')

ax3 = plt.subplot(313)
ax3.axis([-60, -180, -5, 5])
ax3.set_yticks([-5,0,5])
ax3.grid(True)

for i in range(0, len(glon)):
    ax3.scatter(x[i], y[i], marker=marker1[i], color="black",
                edgecolor='black', facecolor='None', s=70, linewidth='1.5')

ax3.set_xlabel('Galactic longitude (deg)', fontsize='15')
</pre>

Out[3]:
<pre>&lt;matplotlib.text.Text at 0x10d1384e0&gt;</pre>

In [4]:
<pre>#print(len(blazar))
#print(len(galactic))
#print(len(unassociated))
</pre>

In [5]:
<pre>#only up to mm for 1FHL
#up to ss for 2FHL

print(len(aa))
print(len(bb))
print(len(cc))
print(len(dd))
print(len(ee))
print(len(ff))
print(len(gg))
print(len(hh))
print(len(ii))
print(len(jj))
print(len(kk))
print(len(ll))
print(len(mm))
print(len(nn))
print(len(oo))
print(len(pp))
print(len(qq))
print(len(rr))
print(len(ss))
</pre> <pre>259
71
58
20
7
11
6
6
3
1
1
0
71
0
0
0
0
0
0
</pre>

In [6]:
<pre>color='black'

no_fhl = mlines.Line2D([], [], color=color, marker='s',
                         markersize=8, label='No association', linestyle='None',
                         markerfacecolor='None', markeredgecolor=color)
pulsar = mlines.Line2D([], [], color=color, marker='*', markersize=8,
                       label='Pulsar', linestyle='None', markerfacecolor='None',
                       markeredgecolor=color)
hmb = mlines.Line2D([], [], color=color, marker='d', markersize=8,
                    label='High-mass binary', linestyle='None', markerfacecolor='None',
                    markeredgecolor=color)
pos_fhl = mlines.Line2D([], [], color=color, marker='&lt;', markersize=8,
                          label='Possible association with SNR or PWN',
                          linestyle='None', markerfacecolor='None',
                          markeredgecolor=color)
glob = mlines.Line2D([], [], color=color, marker='&gt;', markersize=8,
                    label='Globular Cluster', linestyle='None', markerfacecolor='None',
                    markeredgecolor=color)
galaxy = mlines.Line2D([], [], color=color, marker='v', markersize=8,
                       label='Galaxy', linestyle='None',
                       markerfacecolor='None', markeredgecolor=color)
starBur = mlines.Line2D([], [], color=color, marker='p', markersize=8,
                        label='Starburst Galaxy', linestyle='None', markerfacecolor='None',
                        markeredgecolor=color)
nova = mlines.Line2D([], [], color=color, marker='h', markersize=8,
                    label='Nova/LVB/SFR', linestyle='None', markerfacecolor='None',
                    markeredgecolor=color)
agn = mlines.Line2D([], [], color=color, marker='^', markersize=8,
                    label='AGN', linestyle='None', markerfacecolor='None',
                    markeredgecolor=color)
pwn = mlines.Line2D([], [], color=color, marker='D', markersize=8, 
                    label='PWN', linestyle='None', markerfacecolor='None',
                    markeredgecolor=color)
snr = mlines.Line2D([], [], color=color, marker='o', markersize=8,
                    label='SNR', linestyle='None', markerfacecolor='None',
                    markeredgecolor=color)

fontP = FontProperties()
fontP.set_size('small')

for i in range(len(glon)):
    lgd = ax1.legend(handles=[no_fhl, pulsar, hmb, 
                        pos_fhl, glob, galaxy, starBur, nova, 
                        agn, pwn, snr],loc='center left',
                        bbox_to_anchor=(-0.1,1.8), ncol=4, prop = fontP,
                        fancybox=True,shadow=False,numpoints=1)
</pre>

In [7]:
<pre>fig.savefig("grapher1FHL.png", bbox_Extra_artists=(lgd),
           bbox_inches='tight')
print("Ready!")
</pre> <pre>Ready!
</pre>

In [8]:
<pre>#fig=plt.figure()
#ax=fig.add_subplot(111)
#ax.scatter([3,2],[4,5])</pre>

-- %USERSIG{RichardNederlander - 2016-10-23}%

---++ Comments

<br />%COMMENT%
Topic revision: r1 - 2016-10-23 - RichardNederlander
Veritas
Webs
ATLAS
DOE
Main
TWiki
Veritas