Neutrino Generator Tutorial

The easiest way to use simweights is to book your data to hdf5files using tableio.

#!/usr/bin/env python3

# SPDX-FileCopyrightText: © 2022 the SimWeights contributors
#
# SPDX-License-Identifier: BSD-2-Clause

from pathlib import Path

from I3Tray import I3Tray
from icecube import hdfwriter, simclasses

FILE_DIR = Path("/data/sim/IceCube/2016/filtered/level2/neutrino-generator/21217/0000000-0000999/")
files = sorted(str(f) for f in FILE_DIR.glob("Level2_IC86.2016_NuMu.021217.00000*.i3.zst"))

tray = I3Tray()
tray.Add("I3Reader", FileNameList=files)
tray.Add(
    hdfwriter.I3HDFWriter,
    SubEventStreams=["InIceSplit"],
    keys=["PolyplopiaPrimary", "I3MCWeightDict"],
    output="Level2_IC86.2016_NuMu.021217.hdf5",
)

tray.Execute()

Note that one of the booked keys is I3MCWeightDict, the key which contains the information necessary to calculate the weights.

You can check that the hdf5 file was created correctly by running h5ls. The output should look something like this:

$ h5ls Level2_IC86.2016_NuMu.021217.hdf5
I3MCWeightDict           Dataset {7485/Inf}
PolyplopiaPrimary        Dataset {7485/Inf}
__I3Index__              Group

Now we can run a our script which calculates the weights and make a histogram.

#!/usr/bin/env python3

# SPDX-FileCopyrightText: © 2022 the SimWeights contributors
#
# SPDX-License-Identifier: BSD-2-Clause

import pandas as pd
import pylab as plt
from numpy.typing import ArrayLike

import simweights

# load the hdf5 file that we just created using pandas
hdffile = pd.HDFStore("Level2_IC86.2016_NuMu.021217.hdf5", "r")

# instantiate the weighter object by passing the pandas file to it
weighter = simweights.NuGenWeighter(hdffile, nfiles=10)


def northern_track(energy: ArrayLike) -> ArrayLike:
    """This function to represent the IceCube northern track limit.
    Note that the units are GeV^-1 * cm^-2 * sr^-1 * s^-1 per particle type.
    """
    return 1.44e-18 / 2 * (energy / 1e5) ** -2.2


# get the weights by passing the flux to the weighter
weights = weighter.get_weights(northern_track)

# print some info about the weighting object
print(weighter.tostring(northern_track))

# create equal spaced bins in log space
bins = plt.geomspace(1e2, 1e8, 50)

# get energy of the primary cosmic-ray from `PolyplopiaPrimary`
primary_energy = weighter.get_column("PolyplopiaPrimary", "energy")

# histogram the primary energy with the weights
plt.hist(primary_energy, weights=weights, bins=bins)

# make the plot look good
plt.loglog()
plt.xlabel("Primary Energy [GeV]")
plt.ylabel("Event Rate [Hz]")
plt.xlim(bins[0], bins[-1])
plt.ylim(1e-8, 2e-6)
plt.tight_layout()
plt.show()

Note that we need to pass the number of files to NuGenWeighter and that the model is a function that returns a value in units of \(\mathrm{GeV}^{-1}\cdot\mathrm{cm}^{-2}\cdot\mathrm{sr}^{-1}\cdot{s}^{-1}\) per neutrino flavor.

The output should look something like this: