# by P. Kaaret 2/23/2018
import matplotlib.pyplot as plt
from numpy import array, where
import astropy.io.fits

filename = 'acisf03807N003_evt2.fits'

f = pyfits.open(filename) # read in the FITS file with event data
h = f[1].header # copy the FITS header
d = f[1].data # put all of the data into the dictionary d
f.close() # close the FITS file

 # move the data from the dictionary to numpy arrays
x = array(d['x']) # x sky coordinate in pixels
y = array(d['y']) # y sky coordinate in pixels
energy = array(d['energy']) # energy in eV
time  = array(d['time']) # energy in eV
ccd_id = array(d['ccd_id']) # which CCD
grade = array(d['grade']) # event grade

# select events according to energy
elow, ehigh = 300.0, 8000.0 # elow is the lower bound and ehigh is the upper bound
q = where((energy >= elow)&(energy <= ehigh)) # pick out the events with energies between elow and ehigh
x = x[q] # trim down the arrays, keeping only the events that pass the cut
y = y[q] # it's important to do all of the arrays so that the data for each event are kept together
energy = energy[q]
time = time[q]
ccd_id = ccd_id[q]
grade = grade[q]

# make plots
# plot the event Grades
plt.ion()
plt.figure(1)
plt.clf()
plt.xlabel('Event Grade')
plt.ylabel('Number of events')
plt.title(filename)
plt.hist(grade)
plt.show()

# plot the event energies
plt.ion()
plt.figure(2)
plt.clf()
plt.xlabel('Enegy (eV)')
plt.ylabel('Number of events')
plt.title(filename)
plt.hist(energy)
plt.show()


# plot the event positions
plt.ion()
plt.figure(3)
plt.clf()
plt.xlabel('RA')
plt.ylabel('DEC')
plt.title(filename)
plt.plot(x, y, '.')
plt.show()