Weighting Simulations
refrences:
- https://wiki.icecube.wisc.edu/index.php/Weights_in_CORSIKA_Simulation_Data
- http://icecube.wisc.edu/~jvansanten/docs/misc/index.html
- http://software.icecube.wisc.edu/simulation_trunk/projects/weighting/index.html
Dependencies:
Headers
load headers
In [1]:
from I3Tray import *
import icecube
from icecube import icetray, dataio, dataclasses, weighting, simclasses
import numpy as np
import matplotlib as mpl
import pylab
import matplotlib.pyplot as plt
import sys, tables, glob, os, copy
import nuflux
from icecube.icetray import I3Units
import matplotlib.gridspec as gridspec
from icecube.weighting.fluxes import Hoerandel5,GaisserH3a
from icecube.weighting import weighting
%matplotlib inline
weighting nugen:
Load data:
In [14]:
import glob
dataset_dir = '/data/sim/IceCube/2016/filtered/level2/neutrino-generator/21217/'
filelist = list(glob.glob(dataset_dir+'0000000-0000999/Level2_IC86.2016_NuMu.021217.00000*.i3.zst'))
print(len(filelist))
print(filelist)
To weight nugen simulations, type, primary energy, zenith angle
and OneWeight variables of the particle are needed.
They can be read from .i3 file(s).
In [58]:
from icecube.dataclasses import I3Particle
def get_most_energetic_muon(mmclist):
emax = 0
for m in list(mmclist):
p = m.particle
if (p.type in (dataclasses.I3Particle.MuMinus,dataclasses.I3Particle.MuPlus) and p.total_energy > emax):
emax = p.total_energy
return emax
MCtype_nugen=np.array([])
MCenergy_nugen=np.array([])
MCzenith_nugen=np.array([])
OneWeight_nugen=np.array([])
MCmuonEnergy_nugen=np.array([])
for f in filelist:
inFile_nugen = dataio.I3File(f)
#go to the begining of the file
inFile_nugen.rewind()
while(inFile_nugen.more()):
frame=inFile_nugen.pop_physics()
#if(dataclasses.get_most_energetic_muon(frame['I3MCTree_preMuonProp'])!=None):
if('FilterMask' in frame):
if(frame['FilterMask']['MuonFilter_13'].condition_passed):
MCmuonEnergy_nugen=np.append(MCmuonEnergy_nugen,get_most_energetic_muon(frame['MMCTrackList']))
MCtype_nugen=np.append(MCtype_nugen,frame['PolyplopiaPrimary'].type)
MCenergy_nugen=np.append(MCenergy_nugen,frame['PolyplopiaPrimary'].energy)
MCzenith_nugen=np.append(MCzenith_nugen,frame['PolyplopiaPrimary'].dir.zenith)
OneWeight_nugen=np.append(OneWeight_nugen,frame['I3MCWeightDict']['OneWeight'])
You can read off the "NEvents" (number of generated particle in a single file) from "I3MCWeightDict"
In [16]:
nevts=frame['I3MCWeightDict']['NEvents']
print(nevts)
Load the neutrino flux model from nuflux
In [5]:
conventional = nuflux.makeFlux("honda2006")
conventional.knee_reweighting_model = "gaisserH3a_elbert"
Calculate the weights
In [17]:
nfiles = len(filelist) #number of files (in this case only 1)
weights_nugen = np.divide( conventional.getFlux( np.int32(MCtype_nugen) , MCenergy_nugen,np.cos(MCzenith_nugen) ),
(0.5*nevts*nfiles) )
weights_nugen = np.multiply(weights_nugen,OneWeight_nugen)
Plot the true primary and muon energy histogram
In [56]:
fig= plt.figure(figsize=(12,6))
ax = plt.subplot(111)
bins=np.logspace(1,6,25)
null=plt.hist(MCenergy_nugen,bins=bins,histtype='step', weights=weights_nugen, label='Primary Energy')
null=plt.hist(MCmuonEnergy_nugen,bins=bins,histtype='step',weights=weights_nugen, label='Muon Energy')
#null=plt.hist(MCmuonEnergy_nugen,bins=bins)
plt.loglog()
plt.legend()
plt.xlabel("True Energy", fontsize=14)
plt.ylabel("events", fontsize=14)
plt.ylim(1e-6,1e-2)
ax.grid()
for item in (ax.get_xticklabels() + ax.get_yticklabels()):
item.set_fontsize(16)
weighting CORSIKA
load CORSIKA files
In [28]:
corsika_dataset_dir = '/data/sim/IceCube/2016/filtered/level2/CORSIKA-in-ice/20904/'
corsika_filelist = list(glob.glob(corsika_dataset_dir+'0000000-0000999/Level2_IC86.2016_corsika.020904.00000*.i3.zst'))
print(len(corsika_filelist))
print(corsika_filelist)
In [29]:
MCtype_corsika=np.array([])
MCenergy_corsika=np.array([])
#need this for plots
MCmuonEnergy_corsika=np.array([])
for f in corsika_filelist:
inFile_corsika = dataio.I3File(f)
#go to the begining of the file
inFile_corsika.rewind()
while(inFile_corsika.more()):
frame=inFile_corsika.pop_physics()
if('FilterMask' in frame):
if(frame['FilterMask']['MuonFilter_13'].condition_passed):
#Frame may contain coincident events so select injected primary shower 'PolyplopiaPrimary'
MCtype_corsika=np.append(MCtype_corsika,frame['PolyplopiaPrimary'].type)
MCenergy_corsika=np.append(MCenergy_corsika,frame['PolyplopiaPrimary'].energy)
MCmuonEnergy_corsika=np.append(MCmuonEnergy_corsika,get_most_energetic_muon(frame['MMCTrackList']))
Get flux model from weighting library
In [33]:
from icecube.weighting.fluxes import Hoerandel5,GaisserH3a
fluxG = GaisserH3a()
fluxH = Hoerandel5()
Define weighting generator
In [41]:
showers = 1000000 # Generated showers
oversampling = 1
prescalefactor = 500
# Sampling cylinder volume
length = 1200*I3Units.m
radius = 600*I3Units.m
#Energy range
eprimarymin = 600
eprimarymax = 100000000
#Flux normalization and spectral indices
pnorm = np.array([10, 5, 3, 2, 1])
gamma = [-2,-2,-2,-2,-2]
#Energy cutoff
LowerCutoffType = 'EnergyPerNucleon'
UpperCutoffType = LowerCutoffType
generator = weighting.FiveComponent(oversampling*showers,
emin=eprimarymin*I3Units.GeV,
emax=eprimarymax*I3Units.GeV,
normalization=pnorm, gamma=gamma,
LowerCutoffType=LowerCutoffType, UpperCutoffType=UpperCutoffType,
height=length, radius=radius)
weights_GaisserH3a=fluxG(MCenergy_corsika, MCtype_corsika)/generatorHE(MCenergy_corsika, MCtype_corsika)
weights_Hoerandel=fluxH(MCenergy_corsika, MCtype_corsika)/generatorHE(MCenergy_corsika, MCtype_corsika)
weights_GaisserH3a *= prescalefactor
weights_Hoerandel *= prescalefactor
Plot the true primary and muon energy histogram
In [57]:
fig= plt.figure(figsize=(12,6))
ax = plt.subplot(111)
bins=np.logspace(3,9,50)
null=plt.hist(MCenergy_corsika,bins=bins,histtype='step',weights=weights_GaisserH3a, label='GaisserH3a')
null=plt.hist(MCenergy_corsika,bins=bins,histtype='step',weights=weights_Hoerandel, label='Hoerandel')
plt.loglog()
plt.legend()
plt.xlabel("Primary Energy [GeV]",fontsize=14)
plt.ylabel("Event Rate [Hz]",fontsize=14)
plt.ylim(1e-3,3e3)
ax.grid()
for item in (ax.get_xticklabels() + ax.get_yticklabels()):
item.set_fontsize(16)
Calculate the livetime (sum of weights/sum of squared weights)
\begin{equation} T_\textrm{live} = \frac{\sum{w_i}}{\sum{w^2_i}} \end{equation}In [54]:
weightssqr_GaisserH3a=np.power(weights_GaisserH3a,2)
livetime_GaisserH3a=weights_GaisserH3a/weightssqr_GaisserH3a
weightssqr_Hoerandel=np.power(weights_Hoerandel,2)
livetime_Hoerandel=weights_GaisserH3a/weightssqr_Hoerandel
fig= plt.figure(figsize=(12,6))
ax = plt.subplot(111)
bins=np.logspace(3,9,50)
histweightsG,cpuedges=np.histogram(MCenergy_corsika,bins=bins,weights=weights_GaisserH3a)
histweights2G,cpuedges2=np.histogram(MCenergy_corsika,bins=bins,weights=weightssqr_GaisserH3a)
plt.plot(bins[1:],histweightsG/histweights2G,label='Livetime (H3a)')
histweightsH,cpuedges=np.histogram(MCenergy_corsika,bins=bins,weights=weights_Hoerandel)
histweights2H,cpuedges2=np.histogram(MCenergy_corsika,bins=bins,weights=weightssqr_Hoerandel)
plt.plot(bins[1:],histweightsH/histweights2H,label='Livetime (Hoerandel)')
plt.loglog()
plt.legend()
plt.xlabel("Primary Energy [GeV]", fontsize=14)
plt.ylabel("Livetime [s]", fontsize=14)
plt.ylim(1e-3,1e3)
ax.grid()
for item in (ax.get_xticklabels() + ax.get_yticklabels()):
item.set_fontsize(16)
In [ ]: