Triggered CORSIKA 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/CORSIKA-in-ice/21889/0000000-0000999")
files = sorted(str(f) for f in FILE_DIR.glob("Level2_IC86.2016_corsika.021889.00000*.i3.zst"))

tray = I3Tray()
tray.Add("I3Reader", FileNameList=files)
tray.Add(
    hdfwriter.I3HDFWriter,
    SubEventStreams=["InIceSplit"],
    keys=["PolyplopiaPrimary", "I3PrimaryInjectorInfo", "I3CorsikaWeight"],
    output="Level2_IC86.2016_corsika.021889.hdf5",
)

tray.Execute()

Note that two of the booked keys are I3PrimaryInjectorInfo and I3CorsikaWeight, these are the keys which contain 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_corsika.021682.hdf5
I3CorsikaWeight          Dataset {6948/Inf}
I3PrimaryInjectorInfo    Dataset {40/Inf}
PolyplopiaPrimary        Dataset {6948/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

import simweights

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

# instantiate the weighter object by passing the pandas file to it
weighter = simweights.CorsikaWeighter(hdffile)

# create an object to represent our cosmic-ray primary flux model
flux = simweights.GaisserH4a()

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

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

# create equal spaced bins in log space
bins = plt.geomspace(3e4, 1e6, 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(0.1, 10)
plt.savefig("triggered_corsika_tutorial.svg")
plt.show()

The output should look something like this: