notebook.community

Edit and run 





In [ ]:



    


import os
import sys

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc as pltrc
%matplotlib inline

import ROOT
from ROOT import gROOT, gSystem, TFile, TGraphAsymmErrors, TH1, TF1, TFitResultPtr
import rootnotes
import rootprint

import plottingFunctions as pF
from astropy import units as u
from astropy.io import fits, ascii



    











In [ ]:



    


import pandas
import root_numpy
from scipy.interpolate import spline

# set dirs 
home = os.path.expanduser("~")
crabdir = home + "/VERITAS/Crab/"
gcdir = home + "/VERITAS/GC/results/"
plotdir = home + "/Dropbox/GalacticCenter/plots/comparison/"
specdir = home + "/Dropbox/GalacticCenter/spectralPoints/pyformat/"
logdir = home + "/Dropbox/GalacticCenter/log/stage6/"



    











In [ ]:



    


def autoplot(filename, **kwargs):
    """chooses read function and label based on filename"""
    path = os.path
    if not path.isfile(filename):
        raise IOError("file does not exist!")
    base, ext = path.splitext(path.basename)
    
    sP = pF.spectrumPlotter(**kwargs)
    if ext == '.root':
        sP.readVEGASs6Root(filename)
    elif ext == '.csv':
        sP.readCSV()
    elif ext == '.txt':
        sP.readVEGASLog()
    else:
        print "Extenstion ", ext, "not known!"
        
    sP.plotSpectrum(**kwargs)
    
    return sP



    











In [ ]:



    


font = {'family' : 'DejaVu Sans',
        'weight' : 'bold',
        'size'   : 17}

pltrc('font', **font)

#file1 = "SgrA_disp5t_medium_okay_stg5-v255_EA-v255_s6.root"
file1 = "SgrA_disp5t_medium_okay_stg5-v255_EA-v254_s6.root"
#file1 = "SgrA_v255_disp5t_Andy_spectrum_newEA_s6.root"
#file1 = "stage6_Andy_SgrA_spectrum_SgrA_s6.root"
#file2 = "HESS_SgrA_spectral-points_formatted.csv"
file2 = "HESS_SgrA_spectral-points_2009_pulled.csv"

#file3 = 

base, ext = os.path.splitext(file1)
sP1 = pF.spectrumPlotter(c='red')
sP1.readVEGASs6Root(gcdir + file1)
sP1.plotSpectrum(label="VERITAS v254", c='red')

base, ext = os.path.splitext(file2)
sP2 = pF.spectrumPlotter(c='blue')
sP2.readCSV(specdir + file2)
sP2.plotSpectrum(pltFit=False, fontsize=20)
sP2.fitPlot("HESS 2009 fit", energyRange=sP1.energyRange, pltPts=False, fontsize=20)
sP2._plotCommon(fontsize=20)

fig = plt.gcf()
fig.set_size_inches(16, 9)
plt.legend(loc=3)
plt.title("Sgr A*")

#plt.savefig(plotdir + "SgrA_HESS-fit-2009_vs_VERITAS_old-table.png")

V5 vs V6





In [ ]:



    


font = {'family' : 'DejaVu Sans',
        'weight' : 'bold',
        'size'   : 17}
pltrc('font', **font)
EAv = 'v254'

file1 = "SgrA_disp5t_medium_okay_na_stg5-v255_EA-"+EAv+"_s6.root"
file2 = "SgrA_disp5t_medium_okay_ua_stg5-v255_EA-"+EAv+"_s6.root"


base, ext = os.path.splitext(file1)
sP1 = pF.spectrumPlotter(c='red')
sP1.readVEGASs6Root(gcdir + file1)
sP1.plotSpectrum(label="V5")

base, ext = os.path.splitext(file2)
sP2 = pF.spectrumPlotter(c='blue')
sP2.readVEGASs6Root(gcdir + file2)
sP2.plotSpectrum(label="V6")

fig = plt.gcf()
fig.set_size_inches(16, 9)
plt.legend(loc=3)
plt.title("Sgr A*")

#plt.savefig(plotdir + "SgrA_V5-V6_EA"+EAv+"v255_E0-4.png")



    











In [ ]:



    


file1 = "SgrB2_disp5t_medium_okay_stg5-v255_EA-v254_s6.root"
file2 = "HESS_diffuse_spectral-points_formatted.csv"

#file3 = 

base, ext = os.path.splitext(file1)
sP1 = pF.spectrumPlotter(c='red')
sP1.readVEGASs6Root(gcdir + file1)
sP1.plotSpectrum(label="VERITAS")

base, ext = os.path.splitext(file2)
sP2 = pF.spectrumPlotter(c='blue')
sP2.readCSV(specdir + file2)
sP2.plotSpectrum(label="HESS", pltFit=False)

fig = plt.gcf()
fig.set_size_inches(16, 9)
plt.legend(loc=3)
plt.title("Sgr A*")

#plt.savefig(plotdir + "SgrB2_vs-HESS_EAv254.png")



    











In [ ]:



    


file1 = "SgrA_disp5t_medium_okay_stg5-v255_EA-v254_Enorm4_s6.root"
file2 = "SgrA_disp5t_medium_okay_stg5-v255_EA-v254_s6.root"

base, ext = os.path.splitext(file1)
sP1 = pF.spectrumPlotter(c='red')
sP1.readVEGASs6Root(gcdir + file1)
sP1.plotSpectrum(label="norm4")

base, ext = os.path.splitext(file2)
sP2 = pF.spectrumPlotter(c='blue')
sP2.readVEGASs6Root(gcdir + file2)
sP2.plotSpectrum(label="norm1")

fig = plt.gcf()
fig.set_size_inches(16, 9)
plt.legend(loc=3)
plt.title("Sgr A*")

#plt.savefig(plotdir + "SgrB2_vs-HESS_EAv254.png")



    











In [ ]:



    


# Andy vs full runlist 
filename = gcdir + "/SgrA_disp5t_v255_okay-obs-4tel_medium_both_s6.root"
spectrum = pF.spectrumPlotter(c='orange')
spectrum.readVEGASs6Root(filename)
spectrum.plotSpectrum(label="full runlist") #xmin=0.1

spectrum = pF.spectrumPlotter(c='purple')
spectrum.readVEGASs6Root(home+"/VERITAS/GC/results/stage6_Andy_SgrA_spectrum_SgrA_s6.root")
spectrum.plotSpectrum(label="Andy runlist")

plt.title("Sgr A*")
plt.savefig(home+"/Downloads/specPlot.png")

del spectrum



    











In [ ]:



    


file1 = "SgrA_disp5t_medium_okay_stg5-v255_EA-v255_s6.root"
files = '''SgrA_disp5t_medium_okay_stg5-v254_EA-v255_s6.root
SgrA_disp5t_medium_okay_stg5-v255_EA-v254_s6.root
SgrA_disp5t_medium_okay_stg5-v255_EA-v255_Andy-bin-fine4_s6.root'''
#SgrA_v255_disp5t_Andy_spectrum_newEA_s6.root

for f in files.splitlines():

    base, ext = os.path.splitext(file1)
    spec1 = pF.spectrumPlotter(c='red')
    spec1.readVEGASs6Root(gcdir + file1)
    spec1.plotSpectrum(label=base)
    
    base, ext = os.path.splitext(f)
    pathname = gcdir + f
    spec2 = pF.spectrumPlotter(c='blue')
    spec2.readVEGASs6Root(pathname)
    spec2.plotSpectrum(label=base)
    
    plt.title(base)
    plt.savefig(plotdir + base + ".png")
    plt.cla()

    #print raw_input("Press any key to continue..")
    #sys.stdin.read(1)

Log vs ROOT input





In [ ]:



    


spRoot = pF.spectrumPlotter(c='red')
spRoot.readVEGASs6Root(gcdir+"/SgrA_disp5t_v255_okay-obs-4tel_medium_both_s6.root")
spRoot.plotSpectrum(label="ROOT")


spLog = pF.spectrumPlotter(c='blue')
spLog.readVEGASLog(logdir+"/SgrA_disp5t_v255_okay-obs-4tel_medium_both_stage6.txt")
#spLog.readVEGASLog(home+"/Dropbox/GalacticCenter/log/stage6/stage6_Andy_SgrA_spectrum.txt")
spLog.plotSpectrum(label="log")

# larger butterfly due to covariance being present in ROOT file but not log 
#help(plt.errorbar)

Crab comparison for disp 5t





In [ ]:



    


crabdir = home + "/VERITAS/Crab"
#pltrc('font', **font)

for array in ('V5',):
    sp1 = pF.spectrumPlotter(c='red')
    plt.title("Crab "+array)
    fn = crabdir + "/Crab_LZA_disp5t_"+array+"_spectrum_fit1-10TeV_s6.root"
    sp1.readVEGASs6Root(fn)
    sp1.plotSpectrum(label='LZA')
    fn = crabdir + "/Crab_SZA_std_"+array+"_spectrum_fit1-10TeV_s6.root"
    
    sp2 = pF.spectrumPlotter(c='blue')
    sp2.readVEGASs6Root(fn)
    sp2.plotSpectrum(label='SZA')
    
    #fn = crabdir + "/Crab_LZA_disp5t_"+array+"_spectrum_newEA_s6.root"
    #sp.readVEGASs6Root(fn)
    #sp.plotSpectrum(label='newEA', c='orange')
    
fig = plt.gcf()
fig.set_size_inches(16, 9)
plt.legend(loc=3)

plt.savefig(plotdir + "Crab_SZA-LZA_V5.png")

ROOT s6SpectralAnl tests 





In [ ]:



    


def get_hist_points(h):
    """return E[TeV], flux[1/TeV*m^2*s], and flux_err as lists"""
    x, y, yerr = [], [], []
    for pt in range(1, h.GetNbinsX()+1):
        tmpX, tmpY = ROOT.Double(0), ROOT.Double(0)
        #h.GetPoint(pt, tmpX, tmpY)
        x.append(h.GetBinCenter(pt))
        y.append(h.GetBinContent(pt))
        yerr.append(h.GetBinError(pt))
        
    #print np.array(x), np.array(y), np.array(yerr)
    return np.power(10, np.array(x)), np.array(y), np.array(yerr)
# get_hist_points



    











In [ ]:



    


%%rootprint
vegasPath = os.path.expandvars("$VEGAS")
#vegasPath = "/home/mbuchove/Downloads"

# test with ROOT6
#gROOT.Reset()
gSystem.Load("libTreePlayer.so")
gSystem.Load("libPhysics.so")
gSystem.Load(vegasPath + "/common/lib/libSP24sharedLite.so")
gSystem.Load(vegasPath + "/resultsExtractor/lib/libStage6shared.so")
gSystem.AddIncludePath("-Wno-unused -Wno-shadow -Wno-unused-parameter")
gROOT.ProcessLine(".L " + vegasPath + "/common/include/VACommon.h")
gROOT.ProcessLine(".include " + vegasPath + "/common/include/")
gROOT.ProcessLine(".include " + vegasPath + "/resultsExtractor/include/")
gROOT.ProcessLine(".include " + vegasPath + "/cfitsio/include/")
    

#vacomm = VASpectrumAnl()
    
try:
    vegas_class = VACommon()
except:
    print "Unexpected error:", sys.exc_info()[0]
    _use_vegas = False
else:
    _use_vegas = True
    
    
print _use_vegas



    











In [ ]:



    


%%rootprint
import root_numpy

#root_numpy.hist2array
filename = home + "/VERITAS/GC/results/SgrA_disp5t_v255_okay-obs-4tel_medium_both_s6.root"

s6F = TFile(filename, "read")
if not s6F.IsOpen():
    print "Could not open file! ", rootfilename
specAn = s6F.Get("Spectrum/VASpectrumAnl")
specGraph = specAn.GetSpectrumGraph()
xaxis = specGraph.GetXaxis()
#specHist = specAn.GetSpectrumHist()
specHist = specAn.GetRebinnedSpectrumHist()

alpha = specAn.GetAlphaHist()
sig = specAn.GetSigmaHist()
alpha = specAn.GetRebinnedAlphaHist()
#specAn.Rebin(2)

xaxis = specGraph.GetXaxis()

alphaArray = root_numpy.hist2array(alpha)
print alphaArray
sig = root_numpy.hist2array(sig)
print sig

print specGraph.GetN()
print len(alphaArray)
print specHist.GetNbinsX()

print type(specAn)
#specAn.MakeSpectrumGraph()

#s6F.Close()



    











In [ ]:



    


xlab = specHist.GetXaxis().GetTitle()
ylab = specHist.GetYaxis().GetTitle()
print xlab

E, flux, flux_err = get_hist_points(specHist)
#rebinned 

print E
print flux
print flux_err

print type(specGraph)



    











In [ ]:



    


# specGraph 
npoints     = specGraph.GetN()

x2, y2 = [], []
y_err = []
for i in range(npoints):
    tmpX, tmpY = ROOT.Double(0), ROOT.Double(0)
    specGraph.GetPoint(i, tmpX, tmpY)
    x2.append(tmpX)
    y2.append(tmpY)
    y_err.append((specGraph.GetErrorYlow(i), specGraph.GetErrorYhigh(i)))
    
x2 = np.array(x2)
y2 = np.array(y2)
y_err = np.array(y_err)

print x2
print y2
print y_err



    











In [ ]:



    


#%%rootprint
tf1 = specGraph.GetFunction("fFitFunction")
fitnorm = tf1.GetParameter(0)
fitindex = tf1.GetParameter(1)
normenergy = tf1.GetParameter(2)

r = specGraph.Fit(tf1, "S") #TFitResultPtr
#for i in range(3):
#    print tf1.GetParameter(i)
    
cov = r.GetCovarianceMatrix() #TMatrixTSym<double>
cov.Print()

var_norm = cov(0, 0)
var_index = cov(1, 1)
cov_normindex = cov(0, 1) # == (1, 0)

print np.sqrt(var_norm)
print np.sqrt(var_index)
print cov_normindex

#covarr = cov.GetMatrixArray()



    











In [ ]:



    


%%rootprint

cov.Print()

s6F.ls()
s6F.Close()

adjust binning





In [ ]:



    


rebin = 4
binning="""
0.45
0.55
0.65
0.75
0.85
0.95
1.1
1.3
1.5
"""

lastnum = None 
for num in binning.splitlines():
    try:
        num = float(num)

        if lastnum != None:
            delta = ( num - lastnum ) / 4 
            print lastnum
            for i in range(1, rebin):
                print lastnum + i*delta
        
    except ValueError:
        continue 
    
    lastnum = num
    
print num

Misc tests.. 





In [ ]:



    


# pandas 
df = pandas.read_csv(specdir+"HESS_SgrA_spectral-points_formatted.csv", delimiter='\t')
names = [n.split()[0] for n in df.columns]
#units = [u.Unit(n.split()[1]) for n in df.columns]
df.columns = names
print type(df['Energy'])
print np.asarray(df['Energy'])



    











In [ ]:



    


from memory_profiler import memory_profile
#import memory_profiler

#@memory_profile
def firstn(n):
    num = 0
    while num < n:
        yield num
        num += 1
        
#memory_profiler.memory_profile(sum(firstn(100)))
        

firstn(1000000)

print float("inf") >= float("inf")
z = ([0]*5, [1]*5, [2]*5, [3]*5)
print z
z = np.asarray(z)
print z
print (z[2],z[3])
print np.asarray((z[2],z[3])) ### this one!



    











In [ ]:



    


a = [5,2]
try:
    x = 5 
    print len(5)
except Exception as ex:
    template = "An exception of type {0} occurred. Arguments:\n{1!r}"
    message = template.format(type(ex).__name__, ex.args)
    print message

xy = 5
print eval('xy + 5')
print xy

Content source: mbuchove/analysis-tools-m

Similar notebooks:

notebook.community | gallery | about