IceTop CORSIKA Tutorial

IceTop simulations are typically performed by throwing CORSIKA showers uniformly within a circular area defined by a radius (the “sampling radius”), using the “topsimulator” project. Although the high-energy muons in the showers can also be propagated to the deep in-ice detector, it is the sampling radius at the surface that is important for weighting. Each CORSIKA shower is also typically thrown (“resampled”) onto the IceTop array some number of times (usually 100) to ensure that some will trigger the detector. So, IceTop simulations have their own class of S-frame object: the I3TopInjectorInfo, which records:

  • The number of events generated

  • The primary nuclear type

  • The sampling radius

  • The min and max zenith angle

  • The min and max primary energy

  • The power law index

As with the other kinds of weighters, the easiest way to use simweights is to book your data to hdf5files using tableio, if this has not already been done for you:

#!/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

# The following is an example of some "Level3" files for which S-frames were retrofitted
# after production.  But the input can be any IceTop simulation files containing S-frames.
FILE_DIR = Path("/data/ana/CosmicRay/IceTop_level3/sim/IC86.2012/SIBYLL2.1/p/12360_v1s/")
files = sorted(str(f) for f in FILE_DIR.glob("Level3_IC86.2012_SIBYLL2.1_p_12360_E6.0_*.i3.bz2"))

tray = I3Tray()
tray.Add("I3Reader", FileNameList=files)
tray.Add(
    hdfwriter.I3HDFWriter,
    SubEventStreams=["IceTopSplit"],
    keys=["I3TopInjectorInfo"],  # <--- and of course whatever other keys you want... MCPrimary, etc...
    output="Level3_IC86.2012_SIBYLL2.1_p_12360_E6.0.i3.bz2",
)

tray.Execute()

Note that one of the booked keys is I3TopInjectorInfo, which lives in the S-frame and contains the information necessary to calculate the weights. You should see I3TopInjectorInfo listed as a Dataset if you inspect your output file using the h5ls command.

Once you’ve got hdf5 files booked, the simweights project can be used to calculate the weights for an analysis. Cosmic ray analyses typically deal with large numbers of events, stored in not just one file but collections of files which must be combined together. Here is an example:

#!/usr/bin/env python3

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

from pathlib import Path

import numpy as np
import pylab as plt
import tables

import simweights

# —- Define which files you want to be part of the dataset ——
# (This particular example is all low-energy protons)
FILE_DIR = Path("/data/ana/CosmicRay/IceTop_level3/sim/IC86.2012/SIBYLL2.1/p/12360_v1s/h5files/")
bigfilelist = sorted(str(f) for f in FILE_DIR.glob("Level3_IC86.2012_SIBYLL2.1_p_12360_E*.h5"))

# —- Create the “IceTopWeighter” ——
# In this case, it's from a whole bunch of files strung together using the "+" operator.
# (This example uses tables, but one could also use pandas if you prefer)
weighter = None
for filename in bigfilelist:
    file_obj = tables.open_file(filename, "r")
    if weighter is None:
        weighter = simweights.IceTopWeighter(file_obj)
    else:
        weighter += simweights.IceTopWeighter(file_obj)

# -- Choose a flux model --
# This particular one is the "four-component" Gaisser H4a model (p, He, O, Fe).
# Note the "_IT" in the name of the function.  This distinguishes it from the five-component
# (p, He, N, Al, Fe) version of this model.
flux = simweights.GaisserH4a_IT()

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

# Dump some info to the screen about the weighting object, if you like
print(weighter.tostring(flux))

# —- Make a plot of your favorite variable ——
# ...such as primary energy.
primary_energy = weighter.get_column("MCPrimary", "energy")
plt.hist(
    np.log10(primary_energy),
    weights=weights,
    bins=46,
    range=[5.0, 9.6],
    log=True,
    histtype="step",
)
plt.xlabel("log10(Primary Energy [GeV])")
plt.ylabel("Event Rate [Hz]")
plt.show()

Note that many of the cosmic ray flux models in simweights, such as Hoerendel, GaisserH3a, GaisserH4a, GlobalFitGST, etc., are five-component (p, He, N, Al, Fe) models. Many IceTop simulations however use four components (p, He, O, Fe). So not all of these models are usable out of the box with four-component IceTop simulations… two of them (Hoerandel_IT, and GaisserH4a_IT) have been coded up as four-component versions.

To weight to a fixed spectral index rather than a model, use something like this:

powerlawindex = -2.7
weights = weighter.get_weights(lambda energy: energy**powerlawindex)

For the GaisserH4a_IT model, the output should look something like this. This plot was made by combining low- and high-energy samples from each nuclear type (for instance, dataset 12360 + 20143 for the protons), so that it spans the full energy range from 10^5 -> 10^9.6 GeV.

_images/icetop_eprimary_combinedLEHE_2012_H4a.svg

Note that the calculated effective area by simweights is the projected effective area, which can be corrected by dividing by cos(zenith).