Pure python library for calculating the weights of Monte Carlo simulation for IceCube.
SimWeights was designed with goal of calculating weights for IceCube simulation in a way that it is easy to combine combine datasets with different generation parameters into a single sample. It was also designed to be a stand alone project which does not depend on IceTray in any way so that it can be installed easily on laptops. SimWeights gathers all the information it needs form information in the hdf5 file so there is no need for access to the simulation production database. SimWeights works with files produced with corsika-reader, neutrino-generator, and genie-reader.
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
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 HzThe 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 HzTo 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 HzYou 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 HzYou 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.0Simulation 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 HzAlso note that these examples use pandas. SimWeights will work equally well with
pandas, h5py, or pytables.
Full documentation is available on the IceCube Documentation Server.
Please direct any questions to @kjm on the slack channel
#software.
See the contributing guide