Getting Started

The source code for simweights is available on github.

Prerequisites

Required: numpy, scipy

Optional: h5py, pytables, pandas, uproot, nuflux

Installation

To install from pypi run:

pip install simweights

Alternatively, if you need to install unreleased code from main you can run:

pip install git+https://github.com/icecube/simweights.git

On certain installs of python on cvmfs you might get the following error: ModuleNotFoundError: No module named 'glob'. If this happens you can add the following option --no-build-isolation to the above command.

If you want to develop simweights you can install directly with flit. The -s option will symlink the module into site-packages rather than copying it, so that you can test changes without reinstalling the module:

pip install flit
git clone git@github.com:icecube/simweights.git
cd simweights
flit install [--user] -s

Basic Usage

For triggered CORSIKA or CORSIKA produced by corsika-reader with S-Frames files use CorsikaWeighter() without any additional arguments:

>>> import simweights, pandas
>>> simfile = pandas.HDFStore("Level2_IC86.2016_corsika.021889.000000.hdf5", "r")
>>> flux_model = simweights.GaisserH4a()
>>> weight_obj = simweights.CorsikaWeighter(simfile)
>>> weights = weight_obj.get_weights(flux_model)
>>> print(f"Rate = {weights.sum():5.2f} Hz")
Rate = 122.84 Hz

The value returned by get_weights() is the rate of events in Hz

For traditional CORSIKA files made with corsika-reader you will also use simweights.CorsikaWeighter(), but you need to know the number of .i3 files that contributed to create this hdf5 file and pass it as the nfiles parameter.

For neutrino-generator you can use NuGenWeighter() which also requires you to know the number of files. Flux models from nuflux can be used:

>>> import nuflux
>>> simfile = pandas.HDFStore("Level2_IC86.2016_NuMu.020878.000000.hdf5")
>>> flux_model = nuflux.makeFlux("CORSIKA_GaisserH3a_QGSJET-II")
>>> weight_obj = simweights.NuGenWeighter(simfile, nfiles=1)
>>> weights = weight_obj.get_weights(flux_model)
>>> print(f"Rate = {weights.sum():5.2e} Hz")
Rate = 1.41e-02 Hz

To weight a spectrum with a function you can also pass a callable to get_weights()

>>> weights = weight_obj.get_weights(lambda energy: 7.2e-8 * energy**-2.2)
>>> print(f"Rate = {weights.sum():5.2e} Hz")
Rate = 2.34e-05 Hz

You can also pass flux values as a numpy array with the same length as the sample

>>> fluxes = 7.2e-8 * simfile["I3MCWeightDict"]["PrimaryNeutrinoEnergy"] ** -2.2
>>> weights = weight_obj.get_weights(fluxes)
>>> print(f"Rate = {weights.sum():5.2e} Hz")
Rate = 2.34e-05 Hz

You can also pass a scalar to weight all events with the same flux. Passing a value of 1.0 will result in the well known quantity OneWeight divided by the number of events.

>>> OneWeight = weight_obj.get_weights(1.0)
>>> OldOneWeight = simfile["I3MCWeightDict"]["OneWeight"] / (simfile["I3MCWeightDict"]["NEvents"] / 2)
>>> (OneWeight - OldOneWeight).median()
0.0

Simulation created with genie-reader can be weighted with GenieWeighter():

>>> simfile = pandas.HDFStore("genie_reader_NuE.hdf5")
>>> flux_model = nuflux.makeFlux("IPhonda2014_spl_solmax")
>>> weight_obj = simweights.GenieWeighter(simfile)
>>> weights = weight_obj.get_weights(flux_model)
>>> print(f"Rate = {weights.sum():5.2e} Hz")
Rate = 3.78e+00 Hz

Also note that these examples use pandas. SimWeights will work equally well with pandas, h5py, or pytables.

Documentation

Full documentation is available on the IceCube Documentation Server.

Getting Help

Please direct any questions to @kjm on the slack channel #software.