Albrecht keeps asking about this geometry cut, convinced that background is getting in at certain angles, and that a simple radius cut on the position won't take care of the problem. This notebook is to make L8 plots to see if the data/MC disagreement is positional dependent.
In [1]: 
from icecube import icetray, dataclasses, dataio
from icecube.icetray import I3Units
import tables, dashi
import numpy as n
import matplotlib.pyplot as plt
dashi.visual()  

pylab.rcParams['savefig.dpi']=90
In [2]: 
from icecube.OscillationFitter import MixGenieNugen, pathDicts
from icecube.OscillationFitter.MixGenieNugen import VectWeightMod
In [8]: 
data = '/Users/gladstone/IceCube/data/' #copied to laptop
# data = '/Users/gladstone/Desktop/netuserpoint/data/' # connect by sshfs
data = '/Volumes/home/netuser/gladstone/data/' #connect by VPN & mac Finder
print data

/Volumes/home/netuser/gladstone/data/
In [9]: 
def cat(n1, n2, g1, g2, g3):
    return n.concatenate((n1, n2, g1, g2, g3),axis=0)

Get files, weights, and energy (reco & true)

In [10]: 
dfile = tables.openFile(data+"exp/L10_exp0And1.hdf5")
datalivetime = (781.19 + 785.2) * I3Units.hour
eRecoE = dfile.root.Monopod.cols.energy[:]
eRecoZen = cos(dfile.root.Monopod.cols.zenith[:])
eWeights = n.array([1. for x in range(len(eRecoE))])
In [12]: 
# nugenMuf = tables.openFile(data+"baseline/Base_nugenNumu/L8calc_Base_nugenNumu_29Files.hdf5")
# nfiles = 290. 
# nfiles = 100. 

nugenMuf = tables.openFile(data+"baseline/Base_nugenNumu/L8calc_Base_nugenNumu.hdf5")
nfiles = 4000. # only if you used all the available files; they're huge.

nugenMuWeights = n.array(nugenMuf.root.AtmWeightsOsc_mod.cols.honda2006_fogli[:]/nfiles)*datalivetime
nugenMuWeightsNoOsc = n.array(nugenMuf.root.AtmWeights_mod.cols.numu[:]/nfiles)*datalivetime
nugenMuRecoE = nugenMuf.root.Monopod.cols.energy[:]
nugenMuRecoZen = cos(nugenMuf.root.Monopod.cols.energy[:])
nugenMuTrueE = nugenMuf.root.PrimaryNu.cols.energy[:]
In [14]: 
# nugenEf = tables.openFile(data+"baseline/Base_nugenNue/L8calc_Base_nugenNue_38Files.hdf5")
# nfiles = 380.

nugenEf = tables.openFile(data+"baseline/Base_nugenNue/L8calc_Base_nugenNue.hdf5")
nfiles = 2000.

nugenEWeights = n.array(nugenEf.root.AtmWeightsOsc_mod.cols.honda2006_fogli[:]/nfiles)*datalivetime
nugenEWeightsNoOsc = n.array(nugenEf.root.AtmWeights_mod.cols.nue[:]/nfiles)*datalivetime
nugenERecoE = nugenEf.root.Monopod.cols.energy[:]
nugenERecoZen = cos(nugenEf.root.Monopod.cols.zenith[:])
nugenETrueE = nugenEf.root.PrimaryNu.cols.energy[:]
In [15]: 
genieMuf = tables.openFile(data+"baseline/Base_genieNumu/L8calc_Base_genieNuMuAndBar.hdf5")
nfiles = 10000.
genieMuWeights = n.array(genieMuf.root.AtmWeightsOsc_mod.cols.honda2006_fogli[:]/nfiles)*datalivetime
genieMuWeightsNoOsc = n.array(genieMuf.root.AtmWeights_mod.cols.numu[:]/nfiles)*datalivetime
genieMuRecoE = genieMuf.root.Monopod.cols.energy[:]
genieMuRecoZen = cos(genieMuf.root.Monopod.cols.zenith[:])
genieMuTrueE = genieMuf.root.PrimaryNu.cols.energy[:]
In [16]: 
genieEf = tables.openFile(data+"baseline/Base_genieNue/L8calc_Base_genieNuEAndBar.hdf5")
nfiles = 10000.
genieEWeights = n.array(genieEf.root.AtmWeightsOsc_mod.cols.honda2006_fogli[:]/nfiles)*datalivetime
genieEWeightsNoOsc = n.array(genieEf.root.AtmWeights_mod.cols.nue[:]/nfiles)*datalivetime
genieERecoE = genieEf.root.Monopod.cols.energy[:]
genieERecoZen = cos(genieEf.root.Monopod.cols.zenith[:])
genieETrueE = genieEf.root.PrimaryNu.cols.energy[:]
In [17]: 
genieTf = tables.openFile(data+"baseline/Base_genieNutau/L8calc_Base_genieNuTauAndBar.hdf5")
nfiles = 2000.
genieTWeights = n.array(genieTf.root.AtmWeightsOsc_mod.cols.honda2006_fogli[:]/nfiles)*datalivetime
genieTWeightsNoOsc = n.array(genieTf.root.AtmWeights_mod.cols.nutau[:]/nfiles)*datalivetime
genieTRecoE = genieTf.root.Monopod.cols.energy[:]
genieTRecoZen = cos(genieTf.root.Monopod.cols.zenith[:])
genieTTrueE = genieTf.root.PrimaryNu.cols.energy[:]
In [18]: 
NuWeight = cat(nugenMuWeights,nugenEWeights,genieMuWeights,genieEWeights,genieTWeights)
NuWeightNoOsc = cat(nugenMuWeightsNoOsc,nugenEWeightsNoOsc,genieMuWeightsNoOsc,genieEWeightsNoOsc,genieTWeightsNoOsc)
NuRecoE = cat(nugenMuRecoE, nugenERecoE, genieMuRecoE, genieERecoE, genieTRecoE)
NuRecoZen = cat(nugenMuRecoZen, nugenERecoZen, genieMuRecoZen, genieERecoZen, genieTRecoZen)
NuTrueE = cat(nugenMuTrueE, nugenETrueE, genieMuTrueE, genieETrueE, genieTTrueE)
NuTrueCZen = cos(cat(nugenMuf.root.Monopod.cols.zenith[:],nugenEf.root.Monopod.cols.zenith[:],genieMuf.root.Monopod.cols.zenith[:],genieEf.root.Monopod.cols.zenith[:],genieTf.root.Monopod.cols.zenith[:]))

print len(NuWeight)

181460
In [19]: 
cfile = tables.openFile(data+"baseline/Base_corsika100/L8calc_DOM100_corsika.hdf5")
nfiles = 100000.
cRecoE = cfile.root.Monopod.cols.energy[:]
cRecoZen = cos(cfile.root.Monopod.cols.zenith[:])
cWeights = n.array([datalivetime/((10.036*I3Units.s)*(nfiles)) for x in range(len(cRecoE))])
# No "true" energy here because it's a bitch to calculate

Grab cut variables

Reco vertex position: x and y (radius)

In [20]: 
#String 36 coordinates: the center string
#String 36 is at x = 46.29, y = -34.88
#Cut at radius 250 from st36

sX = 46.29
sY = -34.88
In [21]: 
x = cat(nugenMuf.root.Monopod.cols.x[:],nugenEf.root.Monopod.cols.x[:],genieMuf.root.Monopod.cols.x[:],genieEf.root.Monopod.cols.x[:],genieTf.root.Monopod.cols.x[:])
y = cat(nugenMuf.root.Monopod.cols.y[:],nugenEf.root.Monopod.cols.y[:],genieMuf.root.Monopod.cols.y[:],genieEf.root.Monopod.cols.y[:],genieTf.root.Monopod.cols.y[:])
r2 = (x-sX)*(x-sX) + (y-sY)*(y-sY)
RecoR = numpy.sqrt(r2)
Nux = x
Nuy = y

print len(r2), len(x), len(Nux)

181460 181460 181460
In [22]: 
ex = n.array(dfile.root.Monopod.cols.x[:])
ey = n.array(dfile.root.Monopod.cols.y[:])
er2 = (ex-sX)*(ex-sX) + (ey-sY)*(ey-sY)
eRecoR = numpy.sqrt(er2)
print len(er2), len(ex)

1142 1142
In [23]: 
cx = n.array(cfile.root.Monopod.cols.x[:])
cy = n.array(cfile.root.Monopod.cols.y[:])
cr2 = (cx-sX)*(cx-sX) + (cy-sY)*(cy-sY)
cRecoR = numpy.sqrt(cr2)
print len(cr2), len(cx)

248 248

Reco vertex depth: z

In [24]: 
#deployed at depths between 1450 m and 2450 m. We thus want the origin close to a depth of 1950 m. 
#(deployed heights are therefore +500 to -500

RecoZ = cat(nugenMuf.root.Monopod.cols.z[:],nugenEf.root.Monopod.cols.z[:],genieMuf.root.Monopod.cols.z[:],genieEf.root.Monopod.cols.z[:],genieTf.root.Monopod.cols.z[:])
eRecoZ = n.array(dfile.root.Monopod.cols.z[:])
cRecoZ = n.array(cfile.root.Monopod.cols.z[:])

Direct Hits:

In [25]: 
NuNDirC = cat(nugenMuf.root.Monopod_RecoParams_SRT.cols.ndirC[:],nugenEf.root.Monopod_RecoParams_SRT.cols.ndirC[:],genieMuf.root.Monopod_RecoParams_SRT.cols.ndirC[:],genieEf.root.Monopod_RecoParams_SRT.cols.ndirC[:],genieTf.root.Monopod_RecoParams_SRT.cols.ndirC[:])
eNDirC = n.array(dfile.root.Monopod_RecoParams_SRT.cols.ndirC[:])
cNDirC = n.array(cfile.root.Monopod_RecoParams_SRT.cols.ndirC[:])

NString: testing.

In [26]: 
NuNString = cat(nugenMuf.root.NString_OfflinePulses_STW_cRT_DCFid.cols.value[:],nugenEf.root.NString_OfflinePulses_STW_cRT_DCFid.cols.value[:],genieMuf.root.NString_OfflinePulses_STW_cRT_DCFid.cols.value[:],genieEf.root.NString_OfflinePulses_STW_cRT_DCFid.cols.value[:],genieTf.root.NString_OfflinePulses_STW_cRT_DCFid.cols.value[:])
eNString = n.array(dfile.root.NString_OfflinePulses_STW_cRT_DCFid.cols.value[:])
cNString = n.array(cfile.root.NString_OfflinePulses_STW_cRT_DCFid.cols.value[:])

N-early: testing.

In [27]: 
NuNearly = cat(nugenMuf.root.Monopod_RecoParams.cols.nearly[:],nugenEf.root.Monopod_RecoParams.cols.nearly[:],genieMuf.root.Monopod_RecoParams.cols.nearly[:],genieEf.root.Monopod_RecoParams.cols.nearly[:],genieTf.root.Monopod_RecoParams.cols.nearly[:])
eNearly  = n.array(dfile.root.Monopod_RecoParams.cols.nearly[:])
cNearly  = n.array(cfile.root.Monopod_RecoParams.cols.nearly[:])

Define the L8 cut set:

In [28]: 
Emin = 10.
RMax = 1500.
ZMin = -510.
ZMax = -200.
NDirMin = 5
# New cut! Experimental. 
# NStringMin = 2
# NearlyMin = -1.
 
cut  = (NuRecoE>Emin)&( RecoZ>ZMin)&( RecoZ<ZMax)&( RecoR<RMax)&(NuNDirC>NDirMin)&(NuNString>NStringMin)&(NuNearly>NearlyMin)
ecut = ( eRecoE>Emin)&(eRecoZ>ZMin)&(eRecoZ<ZMax)&(eRecoR<RMax)&( eNDirC>NDirMin)&( eNString>NStringMin)&( eNearly>NearlyMin)
ccut = ( cRecoE>Emin)&(cRecoZ>ZMin)&(cRecoZ<ZMax)&(cRecoR<RMax)&( cNDirC>NDirMin)&( cNString>NStringMin)&( cNearly>NearlyMin)

nugenMucut = (nugenMuf.root.Monopod.cols.energy[:]>Emin)&(nugenMuf.root.Monopod.cols.z[:]>ZMin)&(nugenMuf.root.Monopod.cols.z[:]<ZMax)&(sqrt((nugenMuf.root.Monopod.cols.x[:]-sX)*(nugenMuf.root.Monopod.cols.x[:]-sX) + (nugenMuf.root.Monopod.cols.y[:]-sY)*(nugenMuf.root.Monopod.cols.y[:]-sY))<RMax)&(nugenMuf.root.Monopod_RecoParams.cols.ndirC[:]>NDirMin)&(nugenMuf.root.NString_OfflinePulses_STW_cRT_DCFid.cols.value[:]>NStringMin)&(nugenMuf.root.Monopod_RecoParams_SRT.cols.nearly[:]>NearlyMin)
nugenEcut  = ( nugenEf.root.Monopod.cols.energy[:]>Emin)&( nugenEf.root.Monopod.cols.z[:]>ZMin)&( nugenEf.root.Monopod.cols.z[:]<ZMax)&(sqrt(( nugenEf.root.Monopod.cols.x[:]-sX)*( nugenEf.root.Monopod.cols.x[:]-sX) + ( nugenEf.root.Monopod.cols.y[:]-sY)*( nugenEf.root.Monopod.cols.y[:]-sY))<RMax)&( nugenEf.root.Monopod_RecoParams.cols.ndirC[:]>NDirMin)&( nugenEf.root.NString_OfflinePulses_STW_cRT_DCFid.cols.value[:]>NStringMin)&( nugenEf.root.Monopod_RecoParams_SRT.cols.nearly[:]>NearlyMin)
genieMucut = (genieMuf.root.Monopod.cols.energy[:]>Emin)&(genieMuf.root.Monopod.cols.z[:]>ZMin)&(genieMuf.root.Monopod.cols.z[:]<ZMax)&(sqrt((genieMuf.root.Monopod.cols.x[:]-sX)*(genieMuf.root.Monopod.cols.x[:]-sX) + (genieMuf.root.Monopod.cols.y[:]-sY)*(genieMuf.root.Monopod.cols.y[:]-sY))<RMax)&(genieMuf.root.Monopod_RecoParams.cols.ndirC[:]>NDirMin)&(genieMuf.root.NString_OfflinePulses_STW_cRT_DCFid.cols.value[:]>NStringMin)&(genieMuf.root.Monopod_RecoParams_SRT.cols.nearly[:]>NearlyMin)
genieEcut  = ( genieEf.root.Monopod.cols.energy[:]>Emin)&( genieEf.root.Monopod.cols.z[:]>ZMin)&( genieEf.root.Monopod.cols.z[:]<ZMax)&(sqrt(( genieEf.root.Monopod.cols.x[:]-sX)*( genieEf.root.Monopod.cols.x[:]-sX) + ( genieEf.root.Monopod.cols.y[:]-sY)*( genieEf.root.Monopod.cols.y[:]-sY))<RMax)&( genieEf.root.Monopod_RecoParams.cols.ndirC[:]>NDirMin)&( genieEf.root.NString_OfflinePulses_STW_cRT_DCFid.cols.value[:]>NStringMin)&( genieEf.root.Monopod_RecoParams_SRT.cols.nearly[:]>NearlyMin)
genieTcut  = ( genieTf.root.Monopod.cols.energy[:]>Emin)&( genieTf.root.Monopod.cols.z[:]>ZMin)&( genieTf.root.Monopod.cols.z[:]<ZMax)&(sqrt(( genieTf.root.Monopod.cols.x[:]-sX)*( genieTf.root.Monopod.cols.x[:]-sX) + ( genieTf.root.Monopod.cols.y[:]-sY)*( genieTf.root.Monopod.cols.y[:]-sY))<RMax)&( genieTf.root.Monopod_RecoParams.cols.ndirC[:]>NDirMin)&( genieTf.root.NString_OfflinePulses_STW_cRT_DCFid.cols.value[:]>NStringMin)&( genieTf.root.Monopod_RecoParams_SRT.cols.nearly[:]>NearlyMin)
In [29]: 
cutNuWeight = NuWeight[cut]
cutNuWeightNoOsc = NuWeightNoOsc[cut]
cutCWeight = cWeights[ccut]
cutEWeight = eWeights[ecut]

cutNuRecoE = NuRecoE[cut]
cutCRecoE = cRecoE[ccut]
cutERecoE = eRecoE[ecut]

cutNuTrueE = NuTrueE[cut]

print len(NuWeight), len(cutNuWeight)

181460 102342
In [30]: 
print "statistics check: # neutrinos passing ", sum(cut)
print "prediction check: rate per year passing:", 10.*sum(cutNuWeight)
print
print "statistics check: # corsika passing ", sum(ccut)
print "prediction check: rate per year passing:", 10.*sum(cutCWeight)
print
print "statistics check: # exp passing ", sum(ecut)
print "prediction check: rate per year passing:", 10.*sum(cutEWeight)

statistics check: # neutrinos passing  102342
prediction check: rate per year passing: 9923.529923

statistics check: # corsika passing  51
prediction check: rate per year passing: 2865.57596652

statistics check: # exp passing  1057
prediction check: rate per year passing: 10570.0

Plot ALL THE THINGS After Cuts

Alternate: CoG (defaults to Monopod if commented out)

In [24]: 
x = cat(nugenMuf.root.clast.cols.x[:],nugenEf.root.clast.cols.x[:],genieMuf.root.clast.cols.x[:],genieEf.root.clast.cols.x[:],genieTf.root.clast.cols.x[:])
y = cat(nugenMuf.root.clast.cols.y[:],nugenEf.root.clast.cols.y[:],genieMuf.root.clast.cols.y[:],genieEf.root.clast.cols.y[:],genieTf.root.clast.cols.y[:])

ex = n.array(dfile.root.clast.cols.x[:])
ey = n.array(dfile.root.clast.cols.y[:])

cx = n.array(cfile.root.clast.cols.x[:])
cy = n.array(cfile.root.clast.cols.y[:])
In [25]: 
bins = (linspace(-300,300,101),linspace(-300.0,300.0,101)) # small: 6 m / bin
bins = (linspace(-300,300,34),linspace(-300.0,300.0,34))   #big: ~18 m / bin

bins = (linspace(-150,210,61),linspace(-210.0,150.0,61)) # small:  6 m / bin
#bins = (linspace(-150,210,31),linspace(-210.0,150.0,31)) #big:    12 m / bin
#bins = (linspace(-150,210,16),linspace(-210.0,150.0,16)) #bigger: 24 m / bin
In [26]: 
c = figure(figsize=(5,4))

Nu = dashi.factory.hist2d((Nux[cut],Nuy[cut]),bins,weights=cutNuWeight)
ax = Nu.pcolor(norm=Normalize(vmin=1e-2))
#ax = Nu.pcolor(log=True,norm=Normalize(vmin=1e-1, vmax=10))
ylabel("y [m]")
xlabel("x [m]")
xlim(bins[0][0], bins[0][-1])
ylim(bins[1][0], bins[1][-1])
title("Neutrinos")
grid(linestyle="-", color="lightgrey")
cbar = c.colorbar(ax)
In [27]: 
c = figure(figsize=(5,4))
corsika = dashi.factory.hist2d((cx[ccut],cy[ccut]),bins,weights=cutCWeight)
ax = corsika.pcolor()
#ax = corsika.pcolor(log=True)
ylabel("y [m]")
xlabel("x [m]")
xlim(bins[0][0], bins[0][-1])
ylim(bins[1][0], bins[1][-1])
title("Corsika")
grid(linestyle="-", color="lightgrey")
cbar = c.colorbar(ax)
In [28]: 
c = figure(figsize=(5,4))
MC = corsika + Nu
ax = MC.pcolor()
ylabel("y [m]")
xlabel("x [m]")
xlim(bins[0][0], bins[0][-1])
ylim(bins[1][0], bins[1][-1])

title("Total MC")
grid(linestyle="-", color="lightgrey")
cbar = c.colorbar(ax)
In [29]: 
c = figure(figsize=(5,4))
exp = dashi.factory.hist2d((ex[ecut],ey[ecut]),bins)
ax = exp.pcolor()
ylabel("y [m]")
xlabel("x [m]")
xlim(bins[0][0], bins[0][-1])
ylim(bins[1][0], bins[1][-1])

title("Experimental Burn Sample")
grid(linestyle="-", color="lightgrey")
cbar = c.colorbar(ax)
In [30]: 
#from matplotlib.colors import Normalize, LogNorm
c = figure(figsize=(5,4))

myratio = dashi.histfuncs.histratio(exp,MC)
myratio.pcolor(norm=Normalize(vmin=0.5, vmax=2))
#myratio.pcolor(log=True)

#ax = ratio.pcolor()
ylabel("y [m]")
xlabel("x [m]")
xlim(bins[0][0], bins[0][-1])
ylim(bins[1][0], bins[1][-1])

title("Data / MC")
grid(linestyle="-", color="lightgrey")
cbar = c.colorbar(ax)
#cbar.ax.set_yticklabels(['< -1', '0', '> 1'])# vertically oriented colorbar

#cbar = c.colorbar(ax,ticks=[-1, 0, 1])
#cbar.ax.set_yticklabels(['< -1', '0', '> 1'])# vertically oriented colorbar

dashi/histfuncs.py:259: RuntimeWarning: divide by zero encountered in divide
  ratio = h1.bincontent / h2.bincontent
dashi/histfuncs.py:259: RuntimeWarning: invalid value encountered in divide
  ratio = h1.bincontent / h2.bincontent
dashi/histfuncs.py:248: RuntimeWarning: divide by zero encountered in divide
  df_lin = lambda x,xerr,y,yerr:  n.sqrt( (xerr/y)**2 + (x*yerr/(y**2))**2)
dashi/histfuncs.py:248: RuntimeWarning: invalid value encountered in divide
  df_lin = lambda x,xerr,y,yerr:  n.sqrt( (xerr/y)**2 + (x*yerr/(y**2))**2)
In [33]: 















In [49]:




















In []: