In [1]:
from astropy.io import fits
import astropy.io.ascii # separate to not overwrite namespace
from astropy.table import Column
from astropy import units as u
from scipy import optimize
from os.path import expanduser
# from ROOT import TRolke
%pylab inline
#%matplotlib inline
In [2]:
home = expanduser("~")
gc_dir = home + "/Dropbox/GalacticCenter/"
In [3]:
erg2TeV = (u.erg).to(u.TeV)
print(erg2TeV)
pylab.rcParams['figure.figsize'] = (12.0, 6.0)
#matplotlib.rcParams['figure.figsize'] = (12.0, 6.0)
In [79]:
# define our line fitting function
fitfunc = lambda p, x: p[0] + p[1] * (x)
errfunc = lambda p, x, y, err: (y - fitfunc(p, x)) / err
# define our () fitting function
fitfuncECPL = lambda p, x: p[0] + p[1] * np.log(x) - (x) / p[2] # np.log is natural log
errfuncECPL = lambda p, x, y, err: (np.log(y) - fitfuncECPL(p, x)) / (err)
fitfuncECPL_CF = lambda N0, gamma, beta, E: N0 + gamma*E - 1.*np.exp(E) / beta
#these are just copied from http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/python_tutorial.html
f = lambda E, N0, E0, gamma: N0*(E/E0)**(-1.*gamma)
ferr = lambda E, F, N0, N0err, E0, cov_gg: \
F*np.sqrt(N0err**2/N0**2 + ((np.log(E/E0))**2)*cov_gg)
f_ecpl = lambda E,N0,E0,gamma,beta: N0*(E/E0)**(-1.*gamma)*np.exp(-1.*E/beta)
ferr_ecpl = lambda E, F, N0, N0err, E0, cov_gg, b, cov_bb: \
F*np.sqrt(N0err**2/N0**2 + ((np.log(E/E0))**2) * cov_gg + (E/E0)**2 / b**4 * cov_bb)
def plotPythonFit(energy, flux, flux_err_arg, color, name, ECPL=False, E0=1., ls="--"):#, power=0.):
"""fit flux points to a curve then plot
by default it's labeled in units of TeV, and flux / m^2
fluxes are multiplied by E^power """
logx = np.log(energy/E0)
logy = np.log(flux)
#logyerr = np.log(flux_err)
if isinstance(flux_err_arg, tuple):
flux_err = (flux_err_arg[1] + flux_err_arg[0]) / 2
else:
flux_err = flux_err_arg
logyerr = flux_err / flux
if ECPL:
pinit = [-26, -2.25, 10]
out = optimize.leastsq(errfuncECPL, pinit,
args=(energy/E0, flux, flux_err / flux),
full_output=1)
else:
pinit = [-26, -2.25] # nb ln
out = optimize.leastsq(errfunc, pinit,
args=(logx, logy, logyerr),
full_output=1)
# end if else ECPL
pfinal = out[0]
covar = out[1]
print("parameters:")
print(pfinal)
print("covariance matrix")
print(covar)
N0 = np.exp(pfinal[0])
gamma = pfinal[1]
E = np.linspace(energy[0], energy[-1], num=100)
if ECPL:
beta = pfinal[2]
F = f_ecpl(E, N0, E0, -1.*gamma, beta)
chi2 = np.sum((flux - f_ecpl(energy, N0, E0, -1.*gamma, beta))**2/flux_err**2) / (len(energy) - 3)
print("chi^2: "+str(chi2)+'\n')
beta_err = np.sqrt( covar[2][2] ) * chi2 #* N0
gamma_err = np.sqrt( covar[0][0] ) * chi2
N0_err = np.sqrt( covar[1][1] ) * N0 * chi2
cov_gg = gamma_err**2
cov_bb = beta_err**2
Ferr = ferr_ecpl(E, F, N0, N0_err, E0, cov_gg, beta, cov_bb)
fitTitle = (name + ' - N0: {0:.2e} +- {2:.2e}, '\
'G: {1:.2f} +- {3:.2f}, '\
'Ec: {4:.2f} +- {5:.2f}, '
'E0: {6:.0f}').format(float(N0), float(gamma),
float(N0_err), float(gamma_err),
float(beta), float(beta_err), float(E0))
else:
F = f(E, N0, E0, -1.*gamma)
chi2 = np.sum((flux - f(energy, N0, E0, -1.*gamma))**2 / flux_err**2) / (len(energy) - 2)
print(chi2)
gamma_err = np.sqrt( covar[0][0] ) * chi2
N0_err = np.sqrt( covar[1][1] ) * N0 * chi2
cov_gg = gamma_err**2
Ferr = ferr(E, F, N0, N0_err, E0, cov_gg)
fitTitle = (name + ' - N0= {0:.2e} +- {2:.2e}, '\
'gamma= {1:.2f} +- {3:.2f}, '\
'E0: {4:.2f}').format(float(N0), float(gamma),
float(N0_err), float(gamma_err), float(E0))
plt.plot(E, F * (E)**power, color=color, ls=ls, marker="", label=fitTitle)
plt.fill_between(E, (E)**power*(F+Ferr), (E)**power*(F-Ferr), color=color, alpha='0.25')
plt.loglog(nonposy="clip")
plt.errorbar(energy, flux*energy**power, flux_err_arg*energy**power,
color=color, ls='', marker='_') # ,label=name
plt.loglog(nonposy="clip")
plt.xlabel("Energy [TeV]")
# end plotPythonFit
In [5]:
def extract_points(filename):
"""extracts points from prepared spectral file"""
with open(filename) as infile:
# make dict
for line in iter(infile):
line = line.split()
try:
float(line[0]) and float(line[3]) and float(line[4]) and float(line[5])
H_energy = np.append(H_energy, float(line[0]))
H_flux = np.append(H_flux, float(line[3]))
H_err_up = np.append(H_err_up, float(line[4]))
H_err_dwn = np.append(H_err_dwn, float(line[5]))
except ValueError:
continue
# end extract SgrA spectral points from file
In [70]:
def extract_spectral_points_from_log(filename):
"""takes in filename of standard root log file
and returns array of tuples representing spectral points"""
verbose = False
points_filename = filename.replace('_stage6', '_spectral-points')
points_file = open(points_filename, 'w')
fitparams = [0., 0., 0., 0.] # norm, index, norm_err, index_err
#handle, ulfilename = mkstemp()
#UL_file = os.fdopen(handle, 'w')
#UL_file = open(filename.replace('stage6', 'ULs'))
with open(filename) as infile:
foundspec = False
specover = False
for line in infile:
if foundspec and not specover:
if line[0] == '+': # this signifies a point
ls = line.split()
newline = ' '.join(ls[1:]) + '\n'
if verbose:
print(line)
print(newline)
points_file.write(newline)
else:
try:
float(line[0])
#UL_file.write(line + '\n')
except ValueError:
specover = True
elif specover:
ls = line.split()
if ls[0] == '1' and ls[1] == "Norm":
fitparams[0] = float(ls[2])
fitparams[2] = float(ls[3]) # err
elif ls[0] == '2' and ls[1] == "Index":
fitparams[1] = float(ls[2])
fitparams[3] = float(ls[3]) # err
print(fitparams)
break
#elif line == " Bin Energy error Flux error Non Noff Nexcess RawOff Alpha Sig Low Edge High Edge":
else:
ls = line.split()
if len(ls) >= 3 and ls[0] == "Bin" and ls[1] == "Energy" and ls[2] == "error":
foundspec = True
points_file.close()
flux_points = np.genfromtxt(points_filename)
#UL_array = np.genfromtxt(ulfilename)
#UL_file.close()
return flux_points, fitparams #, UL_points
# end extract_spectral_points_from_log
In [6]:
# SgrA spectral points and errors
with open(gc_dir+"spectralPoints/HESS_SgrAstar_SpectralPoints_TeV-cm2.txt") as infile:
H_energy = np.array([])
H_flux = np.array([])
H_err_up = np.array([])
H_err_dwn = np.array([])
for line in iter(infile):
line = line.split()
try:
float(line[0]) and float(line[3]) and float(line[4]) and float(line[5])
H_energy = np.append(H_energy, float(line[0]))
H_flux = np.append(H_flux, float(line[3]))
H_err_up = np.append(H_err_up, float(line[4]))
H_err_dwn = np.append(H_err_dwn, float(line[5]))
except ValueError:
continue
# end extract SgrA spectral points from file
#"/spectralPoints/HESS_Diffuse_SpectralPoints_Bins.txt"
In [7]:
print(H_err_dwn)
In [8]:
VEGAS_Points = """
2.499 6.62e-09 5.64e-10
3.96 1.61e-09 1.01e-10
6.273 4.9e-10 3.7e-11
9.935 1.49e-10 1.43e-11
15.73 4.3e-11 5.43e-12
24.87 6.44e-12 1.55e-12
35.37 5.29e-13 5.86e-13"""
V_Points_fine = astropy.io.ascii.read(VEGAS_Points)
# convert from m^-2 to cm^-2
V_Points_fine['col2'] *= 1e-4
V_Points_fine['col3'] *= 1e-4
V_Points_fine
Out[8]:
In [9]:
power = 0.
plotPythonFit(V_Points_fine['col1'], V_Points_fine['col2'], V_Points_fine['col3'], "blue", "HESS", ECPL=True, E0=1.25)
plt.errorbar(V_Points_fine['col1'], V_Points_fine['col2'] * V_Points_fine['col1']**power,
yerr = V_Points_fine['col3'] * V_Points_fine['col1']**power,
label = "VERITAS 2016 Paper", ls="", marker="+", color="red")
plt.ylim(ymin=1e-17)
plt.ylabel(r" dN/dE [TeV m$^{-2}$ s$^{-1}$]")
Out[9]:
Combine
In [10]:
# updated VEGAS points
# highest energy bins
#
VEGAS_Points = """
2.498 5.71e-09 4.03e-10
3.96 1.7e-09 1.08e-10
6.276 5.48e-10 4.15e-11
9.946 1.69e-10 1.69e-11
15.76 6.07e-11 7.63e-12
24.98 6.81e-12 2.67e-12
39.6 2.39e-13 7.91e-13 """
VEGAS_Points = """
2.498 6.82e-09 5.97e-10
3.959 1.69e-09 1.08e-10
6.274 5.42e-10 4.14e-11
9.943 1.67e-10 1.69e-11
15.76 6.00e-11 7.63e-12
24.98 1.02e-11 2.59e-12
39.59 1.03e-12 7.56e-13 """
V_Points = astropy.io.ascii.read(VEGAS_Points)
# convert from m^-2 to cm^-2
V_Points['col2'] *= 1e-4 #* erg2TeV
V_Points['col3'] *= 1e-4 #* erg2TeV
V_Points
Out[10]:
In [11]:
Andy_Points = '''
2.813 3.44e-13 4.52e-14
3.541 2.17e-13 2.43e-14
4.458 1.23e-13 1.37e-14
5.613 5.13e-14 7.4e-15
7.066 2.72e-14 4.3e-15
8.896 1.27e-14 2.48e-15
12.49 5.8e-15 8.86e-16
19.8 1.44e-15 3.36e-16
31.39 1.22e-16 8.14e-17'''
A_Points = astropy.io.ascii.read(Andy_Points)
A_Points['col2'] *= 1 / erg2TeV
A_Points['col3'] *= 1 / erg2TeV
print(log10(A_Points['col1']*erg2TeV))
#A_Points
In [26]:
print(And_p)
print(A_Points)
print(1.114e-15*(17.3929)**2)
In [12]:
power = 2.
# my VEGAS points
#plotPythonFit(V_Points['col1'], V_Points['col2'], V_Points['col3'], "M2016", "red", ECPL=True, E0=1.0)
plt.errorbar(V_Points['col1'], V_Points['col2'] * V_Points['col1']**power,
yerr = V_Points['col3'] * V_Points['col1']**power,
label = "Matt Buchovecky 2016 Update", ls="", marker="+", color="red")
# HESS points
plotPythonFit(H_energy[:-3], H_flux[:-3], (H_err_dwn[:-3],H_err_up[:-3]),
"blue", "HESS", ECPL=True)
plt.errorbar(H_energy[:-3], H_energy[:-3]**2*H_flux[:-3],
yerr=H_energy[:-3]**2*(H_err_dwn[:-3], H_err_up[:-3]),
marker="+", ls="", color="blue", label="HESS")
plt.errorbar(H_energy[-3:], H_energy[-3:]**2*H_flux[-3:],
yerr=(H_energy[-3:]**2*H_err_up[-3:], H_energy[-3:]**2*H_err_dwn[-3:]),
marker="_", ls="", uplims=True, color="blue")
# Andy's points
plt.errorbar(A_Points['col1'], A_Points['col2'] * A_Points['col1']**power,
yerr = A_Points['col3'] * A_Points['col1']**power,
label = "VERITAS 2016 Paper", ls="", marker="+", color="gray")
# plot format and save
plt.title("Sgr A* Spectrum ")
plt.loglog(nonposy="clip")
plt.ylim(ymin=1e-15)
plt.legend(loc="best")
plt.xlabel("Energy [TeV]")
plt.ylabel(r"E$^2$ dN/dE [TeV cm$^{-2}$ s$^{-1}$]")
plt.savefig(gc_dir+"/plots/spectra/SgrA_spectra_HESSoverlay_wAndy_ECPL.png")
In [49]:
power = 2.
# pulled
#plt.plot(And_p[:,0], And_p[:,1]*And_p[:,0]**(power-2.)*erg2TeV, label="Andy 2016 - pulled from paper", ls="", marker="+")
# A_Points[:,2]*And_p[:,0]**(power-2.)*erg2TeV,
# sent
#plt.errorbar(A_Points['col1'], A_Points['col2']*1e4*A_Points['col1']**(power)*erg2TeV, yerr=A_Points['col3']*1e4*A_Points['col1']**(power)*erg2TeV, label="VERITAS 2016 Paper", ls="", marker="_")
#print(A_Points)
#plt.errorbar(H_energy[:-3], H_energy[:-3]**2*H_flux[:-3]*1e4, yerr=H_energy[:-3]**2*(H_err_dwn[:-3], H_err_up[:-3])*1e4, marker="_", ls="", label="HESS - points sent to me")
#/erg2TeV
#plt.errorbar(V5[0], V5[1], V5[2], marker='+', label='V5', ls='')
#plt.errorbar(V6[0], V6[1], V6[2], marker='+', label='V6', ls='')
plt.errorbar(allOff[0], allOff[1], allOff[2], marker='+', label='My 2016 analysis', ls='')
plt.errorbar(mine_rl[0], mine_rl[1], mine_rl[2], marker='+', label='My analsyis w/ Andys runlist', ls='')
msb_nd = np.asarray(msb)
plt.errorbar(msb_nd[0], msb[1], msb[2], marker='+', label='My 2016 analysis w/ diff spectral binning', ls='')
plt.plot(A_c_a[:,0], A_c_a[:,1]*1e4*erg2TeV, label='Andys results', ls='', marker='+')
#plt.plot(A_c_m[:,0], A_c_m[:,1]*1e4*erg2TeV, label='matt runlist', ls='', marker='+')
plt.title("Sgr A* Spectrum ")
plt.loglog()
#plt.ylim(ymin=1e-15)
plt.legend(loc="best")
plt.xlabel("Energy [TeV]")
plt.ylabel(r"E$^2$ dN/dE [TeV m$^{-2}$ s$^{-1}$]")
Out[49]:
In [13]:
diffuse_points = np.genfromtxt(gc_dir+"/spectralPoints/HESS_diffuse_spectrum_points_transpose.csv")
diffuse_points *= 1e3 # to go from cm^2 to m^2, and account for factor of 10 on plot
diffuse_points[0] /= 1e3
# not needed anymore
diffuse_err_up = diffuse_points[2] - diffuse_points[1]
diffuse_err_down = diffuse_points[1] - diffuse_points[3]
diffuse_points[2] = diffuse_err_down
diffuse_points[3] = diffuse_err_up
print(diffuse_points)
#np.savetxt(home+"/Downloads/HESS_diffuse_spectrum_E2flux_TeV-m2.csv", diffuse_points, delimiter='\t')
In [14]:
power = 2.
mult_factor = 10 # to put diffuse and Sgr B2 closer
# transpose so each variable is a list
diffuse_points = np.genfromtxt(gc_dir+"/spectralPoints/HESS_diffuse_spectrum_E2flux_TeV-m2.csv")
diffuse_points *= mult_factor # to put it closer for comparison
diffuse_points[0] /= mult_factor
# values are E^2 * flux
SgrB2_points = np.genfromtxt(gc_dir+"/spectralPoints/SgrB2_spectral_flux_TeV-m2.txt")
#SgrB2_points *= 1e-4 #
#SgrB2_points[:,0] *= 1e4 # don't want to adjust energy
# values are just flux
plt.errorbar(diffuse_points[0], diffuse_points[1]*diffuse_points[0]**(power-2.),
yerr=(diffuse_points[2]*diffuse_points[0]**(power-2.),
diffuse_points[3]*diffuse_points[0]**(power-2.)),
marker='+', ls='', color='red', label='HESS Diffuse')
plt.errorbar(SgrB2_points[:,0], SgrB2_points[:,1]*SgrB2_points[:,0]**power,
yerr=SgrB2_points[:,2]*SgrB2_points[:,0]**power,
marker='_', ls='', color='blue', label='SgrB2')
plotPythonFit(diffuse_points[0], diffuse_points[1]/diffuse_points[0]**2,
(diffuse_points[2]/diffuse_points[0]**2,diffuse_points[3]/diffuse_points[0]**2),
name="HESS Diffuse", color='red', ls='')
plotPythonFit(SgrB2_points[:,0], SgrB2_points[:,1], SgrB2_points[:,2],
name='SgrB2', color='blue', ls='')
E_SgrB2 = np.linspace(SgrB2_points[0,0], SgrB2_points[-1,0], 100)
flux_SgrB2 = 3.522e-9 * np.power(E_SgrB2, -1.932+power)
plt.plot(E_SgrB2, flux_SgrB2, color='blue', ls='-', marker='',
label="Sgr B2: N0=3.522e-9+-1.178e-9 G=-1.932+-0.1672")
E_diffuse = np.linspace(diffuse_points[0,0], diffuse_points[0,-1], 100)
flux_diffuse = 1.92e-8 * np.power(E_diffuse, -2.32+power)
plt.plot(E_diffuse, flux_diffuse*mult_factor, color='red', ls='-', marker='',
label="HESS: N0=(1.92+-0.08stat+-0.28sys)e-8 G=-2.32+-0.05stat+-0.11sys")
plt.title("Sgr B2 / Diffuse Spectrum")
plt.loglog()
#plt.ylim(ymin=3e-14)
plt.legend(loc="best")
plt.xlabel("Energy [TeV]")
plt.ylabel(r"E$^2$ dN/dE [TeV m$^{-2}$ s$^{-1}$]")
plt.savefig(gc_dir+"/plots/spectra/SgrB2_diffuse_spectra_HESSoverlay.png")
In [15]:
power = 0.
# all energies in TeV
G09_points_M2016 = astropy.io.ascii.read(gc_dir+"/spectralPoints/G09+01_allOff_flux_TeV-m2.txt")
# need to fix HESS points, to give actual size of error bars
G09_points_HESS = astropy.io.ascii.read(gc_dir+"/spectralPoints/G09+01_HESS_2005_flux_TeV-cm2.csv")
G09_points_Andy = astropy.io.ascii.read(gc_dir+"spectralPoints/G09+01_Andy_email_flux_m2.txt")
#G09_points_Andy = astropy.io.ascii.read(gc_dir+"spectralPoints/G09+01_Andy_spectral_points_E2-ergs.txt")
# convert cm^-2 to m^-2
G09_points_HESS['col2'] *= 1e4
G09_points_HESS['col3'] *= 1e4
G09_points_HESS['col4'] *= 1e4
G09_points_Andy['col2'] *= 1e4
G09_points_Andy['col3'] *= 1e4
print(G09_points_HESS)
#plt.errorbar(G09_points_M2016['col1'], G09_points_M2016['col2'], G09_points_M2016['col3'],
# label='M2016', ls='', marker='_')
#plt.errorbar(G09_points_HESS['col1'], G09_points_HESS['col2'],
# (G09_points_HESS['col2']-G09_points_HESS['col4'], G09_points_HESS['col3']-G09_points_HESS['col2']),
# label="HESS", ls='', marker='_')
#plt.errorbar()
#plt.errorbar(G09_points_Andy['col1'], G09_points_Andy['col2']/erg2TeV, G09_points_Andy['col3']/erg2TeV, label="Andy", ls='')
plotPythonFit(G09_points_HESS['col1'], G09_points_HESS['col2'],
(G09_points_HESS['col2']-G09_points_HESS['col4'], G09_points_HESS['col3']-G09_points_HESS['col2']),
color='red', name='HESS')
plotPythonFit(G09_points_M2016['col1'], G09_points_M2016['col2'],
G09_points_M2016['col3'],
color='blue', name='M2016')
plotPythonFit(G09_points_Andy['col1'], G09_points_Andy['col2']*1e-4/erg2TeV, G09_points_Andy['col3']*1e-4/erg2TeV,
name="Andy", color='green')
plt.title("G0.9+0.1 Spectrum ")
plt.loglog(nonposy="clip")
plt.legend(loc="best")
#plt.xlabel("Energy [TeV]")
plt.xlim(xmin=0.15, xmax=20)
# think this is just dN/dE
plt.ylabel(r"dN/dE [TeV m$^{-2}$ s$^{-1}$]")
plt.savefig(gc_dir+"/plots/spectra/G09_spectra_HESSoverlay_wAndy.png")
In [16]:
plotPythonFit(G09_points_M2016['col1'], G09_points_M2016['col2'], G09_points_M2016['col3'],
'blue', name='M2016')
In [98]:
#def
power = 2.
plt.rcParams["figure.figsize"] = (16, 9)
plt.ylabel(r"E^2 dN/dE [TeV m$^{-2}$ s$^{-1}$]")
crab_dir = home + "/Dropbox/VEGAS/Crab"
logfile = home + "/Dropbox/VEGAS/NERSC/validation/stage6/Crab_validation_V5_medium_rc6_stage6.txt"
# disp 5t
sza_points = np.genfromtxt(crab_dir+"/spectralPoints/spectral_points_Crab_SZA.txt")
lza_points = np.genfromtxt(crab_dir+"/spectralPoints/spectral_points_Crab_LZA.txt")
# fit parameters from VEGAS output
sza_params = [3.133e-7, -2.427, 1.470e-8, 0.04705] # norm, index, norm_err, index_err
lza_params = [3.157e-7, -2.525, 1.584e-8, 0.04649]
flux_sza = sza_params[0] * np.power(E_sza, sza_params[1]+power)
flux_lza = lza_params[0] * np.power(E_lza, lza_params[1]+power)
# standard analysis
standard_points, params_std = extract_spectral_points_from_log(logfile)
E_std = standard_points[:,1]
Epow = np.power(E_std, power) # e.g. E^2 dN/dE
y_std = standard_points[:,3] * Epow
yerr_std = standard_points[:,4] * Epow
std_label = ("SZA-std - N0={0:.2e} +- {1:.2e} gamma={2:.2f} +- {3:.2f}")
std_label = std_label.format(params_std[0], params_std[2], params_std[1], params_std[3])
plt.errorbar(E_std, y_std, yerr_std, ls='', color='red', label=std_label)
E = np.linspace(E_std[0], E_std[-1], num=100)
plt.plot(E, params_std[0]*np.power(E, params_std[1]+power), color='red', ls='-')
flux_upper = (params_std[0]+params_std[2])*np.power(E, params_std[1]+params_std[3]+power)
flux_lower = (params_std[0]-params_std[2])*np.power(E, params_std[1]-params_std[3]+power)
plt.fill_between(E, flux_upper, flux_lower, color='red', alpha='0.25')
plt.title("Crab Spectrum, Disp 5t vs standard")
plt.xlabel("Energy [TeV]")
plt.xlim(sza_points[0,0]/1.5, lza_points[-1,0]*1.5)
plt.ylim(2e-8, 1e-6)
plt.ylabel(r"E^"+str(power)+"dN/dE [TeV m$^{-2}$ s$^{-1}$]")
plt.loglog()
E_sza = np.linspace(sza_points[0,0], sza_points[-1,0], 100)
E_lza = np.linspace(lza_points[0,0], lza_points[-1,0], 100)
plotPythonFit(sza_points[:,0], sza_points[:,1], sza_points[:,2], name='SZA-disp5t', color='blue')
plotPythonFit(lza_points[:,0], lza_points[:,1], lza_points[:,2], name='LZA-disp5t', color='green')
#plotPythonFit(standard_points[:,1], standard_points[:,3], standard_points[:,4], name='standard', color='red')
#plt.plot(E_sza, flux_sza, color='blue', ls='-', marker='', label="SZA: N0=3.133e-7+-1.47e-8 G=-2.427+-0.04705")
#plt.plot(E_lza, flux_lza, color='green', ls='-', marker='', label="LZA: N0=3.157e-7+-1.584e-8 G=-2.525+-0.04649")
#plt.fill_between(E, (E)**power*(sza_params[+Ferr), (E)**power*(F-Ferr), color=color, alpha='0.25')
plt.legend(loc="best")
plt.savefig(home+"/Dropbox/VEGAS/Crab/plots/Crab_disp5t_SZAvLZA_spectrum_E"+str(power)+"dNdE.png")
# add upper limit
# :,0 gives energy - then flux, error
#plt.errorbar(sza_points[:,0],
# sza_points[:,1]*sza_points[:,0]**power,
# sza_points[:,2]*sza_points[:,0]**power,
# label='SZA', ls='', color='blue', marker='_')
#plt.errorbar(lza_points[:,0],
# lza_points[:,1]*lza_points[:,0]**power,
# lza_points[:,2]*lza_points[:,0]**power,
# label='LZA', ls='', color='green', marker='_')
In [38]:
from matplotlib import pyplot as plt
%matplotlib inline
def plot_all_epochs(cuts):
""""""
epochs = ('V4', 'V5', 'V6')
logdir = home + "/Dropbox/VEGAS/NERSC/validation/stage6"
plotdir = home + "/Dropbox/VEGAS/NERSC/validation/plots"
plt.clf()
plt.loglog()
plt.title("Crab spectrum: " + cuts + " cuts")
plt.xlabel("Energy (TeV)")
plt.ylabel("Flux [g's/m^2/TeV/s]")
for epoch in epochs:
base = "Crab_validation_" + epoch + '_' + cuts + "_rc6"
fn = logdir + "/" + base + "_stage6.txt"
print(fn)
flux_points, fitparams = extract_spectral_points_from_log(fn)
label = "Norm: " + str(fitparams[0]) + " Index: " + str(fitparams[1])
bins = flux_points[:,0].astype(np.int)
energy = flux_points[:,1]
#energyerr = flux_points[:,2]
flux = flux_points[:, 3]
fluxerr = flux_points[:, 4]
plot = plt.errorbar(energy, flux, fluxerr, ls='', label=label)
# loop over epochs
plt.legend(loc='best', ncol=1)
plotname = plotdir + "/Crab_validation_rc6_" + cuts + ".png"
plt.savefig(plotname)
# plot_all_epochs
In [39]:
all_cuts = ('medium', 'hard', 'soft', 'loose')
for cut in all_cuts:
plot_all_epochs(cut)
In [ ]:
import collections
from collections import namedtuple
from tempfile import mkstemp
import os
flux_point = collections.namedtuple('flux_point', " bin energy energyerr flux fluxerr Non Noff Nexcess RawOff alpha sig eLow eHigh")
#file = open()
#points_file = os.open(pfilename)
#handle, pfilename = mkstemp()
#points_file = os.fdopen(handle, 'w')
#points_file.seek(0)
#points_file.read()
#flux_points = np.genfromtxt(open(points_file))
#flux_points = np.genfromtxt(pfilename)
#points_file.delete()
#
#print(flux_points)
#print(UL_points)
In [43]:
# define our (line) fitting function
fitfunc = lambda p, x: p[0] + p[1] * (x)
errfunc = lambda p, x, y, err: (y - fitfunc(p, x)) / err
# define our (line) fitting function
fitfuncECPL = lambda p, x: p[0] + p[1] * np.log(x) - (x) / p[2]
errfuncECPL = lambda p, x, y, err: (np.log(y) - fitfuncECPL(p, x)) / (err)
fitfuncECPL_CF = lambda N0, gamma, beta, E: N0 + gamma * E -1.*np.exp(E) / beta
#these are just copied from http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/python_tutorial.html
f = lambda E,N0,E0,gamma: N0*(E/E0)**(-1.*gamma)
ferr = lambda E,F,N0,N0err,E0,cov_gg: F*np.sqrt(N0err**2/N0**2 + ((np.log(E/E0))**2)*cov_gg)
f_ecpl = lambda E,N0,E0,gamma,beta: N0*(E/E0)**(-1.*gamma)*np.exp(-1.*E/beta)
ferr_ecpl = lambda E,F,N0,N0err,E0,cov_gg,b,cov_bb: F*np.sqrt(N0err**2/N0**2 \
+ ((np.log(E/E0))**2) * cov_gg \
+ (E/E0)**2 / b**4 * cov_bb)
def plotPythonFit_ralph(energy, flux, fluxEr, color, cuts='nocuts', ECPL = False, E0 = 1):
logx = np.log(energy/E0)
logy = np.log(flux)
logyerr = fluxEr / flux
if ECPL:
pinit = [-26, -2.25, 10]
out = optimize.leastsq(errfuncECPL, pinit,
args=(energy/E0, flux, fluxEr / flux),
full_output=1)
# # print out
# # pinit = [out[0][0], -1.*out[0][1], out[0][2]]
# # print pinit
# out = optimize.curve_fit(fitfuncECPL_CF, logx, logy,
# p0=pinit,
# sigma=fluxEr,
# bounds = ([-30.,-3.,1.],[-25.,-2.,1e2]))
# absolute_sigma = True,
# print out
# out = optimize.leastsq(errfuncECPL, pinit,
# args=(energy, flux, fluxEr / flux),
# full_output=1)
else:
pinit = [-26, -2.25] # nb ln
out = optimize.leastsq(errfunc, pinit,
args=(logx, logy, logyerr),
full_output=1)
pfinal = out[0]
covar = out[1]
print (pfinal)
print (np.diag(covar))
N0 = np.exp(pfinal[0])
gamma = pfinal[1]
E = np.linspace(energy[0],energy[-1],num=100)
if ECPL:
beta = pfinal[2]
F = f_ecpl(E,N0,E0, -1. * gamma, beta)
chi2 = np.sum((flux - f_ecpl(energy,N0,E0, -1. * gamma, beta))**2/fluxEr**2) / (len(energy) - 3)
print(chi2)
beta_err = np.sqrt( covar[2][2] ) * chi2 #* N0
gamma_err = np.sqrt( covar[0][0] ) * chi2
N0_err = np.sqrt( covar[1][1] ) * N0 * chi2
cov_gg = gamma_err**2
cov_bb = beta_err**2
Ferr = ferr_ecpl(E,F,N0,N0_err,E0,cov_gg,beta,cov_bb)
fitTitle = (cuts + ' - N0: {0:.2e} +- {2:.2e}, '\
'G: {1:.2f} +- {3:.2f}, '\
'Ec: {4:.2f} +- {5:.2f}, '
'E0: {6:.0f}').format(float(N0), float(gamma),
float(N0_err), float(gamma_err),
float(beta), float(beta_err), float(E0))
else:
F = f(E,N0,E0, -1. * gamma)
chi2 = np.sum((flux - f(energy, N0, E0, -1. * gamma))**2/fluxEr**2) / (len(energy) - 2)
print (chi2)
gamma_err = np.sqrt( covar[0][0] ) * chi2
N0_err = np.sqrt( covar[1][1] ) * N0 * chi2
cov_gg = gamma_err**2
Ferr = ferr(E,F,N0,N0_err,E0,cov_gg)
fitTitle = (cuts + ' - N0: {0:.2e} +- {2:.2e}, '\
'gamma: {1:.2f} +- {3:.2f}, '\
'E0: {4:.2f}').format(float(N0), float(gamma),
float(N0_err), float(gamma_err), float(E0))
plt.plot(E, F * (E)**power, color=color, ls="--", marker="", label = fitTitle)
plt.fill_between(E, (E)**power*(F+Ferr), (E)**power*(F-Ferr), color=color, alpha='0.25')
In [22]:
power = 2.
err_bar_red = sqrt(2) # projection for doubled dataset
# my VEGAS points
#plotPythonFit(V_Points['col1'], V_Points['col2'], V_Points['col3'], name="M2016", color="gray", ECPL=True, E0=1.0)
plt.errorbar(V_Points['col1'], V_Points['col2'] * V_Points['col1']**power,
yerr = V_Points['col3'] * V_Points['col1']**power,
label = "Matt Buchovecky 2016 Update", ls="", marker="_", color="gray")
# HESS points
plt.errorbar(H_energy[:-3], H_energy[:-3]**2*H_flux[:-3],
yerr=H_energy[:-3]**2*(H_err_dwn[:-3], H_err_up[:-3]),
marker="+", ls="", color="blue", label="HESS")
plt.errorbar(H_energy[-3:], H_energy[-3:]**2*H_flux[-3:],
yerr=(H_energy[-3:]**2*H_err_up[-3:], H_energy[-3:]**2*H_err_dwn[-3:]),
marker="_", ls="", uplims=True, color="blue")
# Andy's points
And_p = np.genfromtxt(gc_dir+"/spectralPoints/SgrA_Andy2016_E2flux_erg-m2.csv")
plt.plot(And_p[:,0], And_p[:,1] * And_p[:,0]**(power-2.),
# yerr = And_p[:,] * And_p[:,0]**(power-2.),
label = "VERITAS 2016 Paper", ls="", marker="+", color="green")
# plot format and save
plt.title("Sgr A* Spectrum ")
plt.loglog(nonposy="clip")
#plt.ylim(ymin=1e-15)
plt.legend(loc="best")
plt.xlabel("Energy [TeV]")
plt.ylabel(r"E$^2$ dN/dE [TeV cm$^{-2}$ s$^{-1}$]")
print(log10(And_p[:,0]))
plt.savefig(gc_dir+"/plots/spectra/SgrA_spectra_HESSoverlay_wAndy_ECPL_projected.png")
In [55]:
def extract_spectral_points(logfile, power=2.):
"""supply a stage 6 log file and get the spectral points from it"""
points = [[], [], []]
bin = 0
begin = False
with open(logfile) as file:
for line in file:
split = line.split()
if len(split) > 1 and split[0] == "Bin" and split[1] == "Energy":
begin = True
elif "FCN=" in line:
begin = False
elif begin and split[0] == '+':
points[0].append(float(split[2]))
points[1].append(float(split[4])*float(split[2])**power)
points[2].append(float(split[5])*float(split[2])**power)
bin += 1
return points
V5 = extract_spectral_points(gc_dir+"/log/stage6/SgrA_V5_disp5t_4tel_stage6.txt")
V6 = extract_spectral_points(gc_dir+"/log/stage6/SgrA_V6_disp5t_4tel_stage6.txt")
allOff = extract_spectral_points(gc_dir+"/log/stage6/SgrA_test_allOff_stage6.txt")
mine_rl = extract_spectral_points(gc_dir+"/log/stage6/stage6_Andy_SgrA_spectrum.txt") # me running Andy's runlist
msb = extract_spectral_points(gc_dir+"/log/stage6/SgrA_bin_Andy_no69365_stage6.txt")
# comparison with Andy
A_c_a = np.genfromtxt(gc_dir+"spectralPoints/SgrA_spectrum_Andy_runlist_comparison_E2flux_ergs-cm2.csv")
A_c_m = np.genfromtxt(gc_dir+"spectralPoints/SgrA_spectrum_Matt_runlist_comparison_E2flux_ergs-cm2.csv")
print(A_c_a[1]*1e4*erg2TeV)
#" Bin Energy error Flux error Non Noff Nexcess RawOff Alpha Sig Low Edge High Edge":