Weighting

In order to weight generated muon bundles to a flux, you need to know both the distribution of events in the flux model (parameterized in the tables provided with MuonGun) and the distribution of events that you generated (calculated by the Generator, stored in a special “S” frame at the beginning of every generated file). Given those, you can calculate weights either within IceTray using the WeightCalculatorModule I3Module or from a standalone Python script using the WeightCalculator class.

First, you should collect the generators for all of the files you plan to use:

def harvest_generators(infiles):
    """
    Harvest serialized generator configurations from a set of I3 files.
    """
    from icecube.icetray.i3logging import log_info as log
    generator = None
    for fname in infiles:
        f = dataio.I3File(fname)
        fr = f.pop_frame(icetray.I3Frame.Stream('S'))
        f.close()
        if fr is not None:
            for k in fr.keys():
                v = fr[k]
                if isinstance(v, MuonGun.GenerationProbability):
                    log('%s: found "%s" (%s)' % (fname, k, type(v).__name__), unit="MuonGun")
                    if generator is None:
                        generator = v
                    else:
                        generator += v
    return generator

The generators can simply be added together, or multiplied by an integer to represent a larger number of identically-configured generators.

Weighting to a cosmic ray flux model

You can pass this combined generator to WeightCalculatorModule to calculate a weight appropriate for the combined set of files:

model = MuonGun.load_model('GaisserH4a_atmod12_SIBYLL')
generator = harvest_generators(infiles)
tray.AddModule('I3MuonGun::WeightCalculatorModule', 'MuonWeight', Model=model,
    Generator=generator)

This will put an I3Double called “MuonWeight” into the frame that represents a weight in events per second. Alternatively, you can use the provided tableio converter to write the parameters needed for the weight calculation to a table:

from icecube.hdfwriter import I3HDFWriter
tray.AddSegment(I3HDFWriter, 'scribe',
    Output=outfile,
    Keys=[dict(key='I3MCTree', name='BundleParameters',
             converter=MuonGun.converters.MuonBundleConverter(1, generator.surface))],
    Types=[],
    SubEventStreams=['nullsplit'],
)

and then use the standalone WeightCalculator class to calculate a weight:

model = MuonGun.load_model('GaisserH4a_atmod12_SIBYLL')
generator = harvest_generators(infiles)
weighter = MuonGun.WeightCalculator(generator.surface, model, generator)

with tables.openFile(outfile) as hdf:
    axis = hdf.root.MCPrimary.read()
    bundle = hdf.root.BundleParameters.read()
    weights = weighter(axis['x'], axis['y'], axis['z'], axis['zenith'], axis['azimuth'],
        bundle['multiplicity'], bundle['energy'], bundle['radius'])

Note

The weighter will only be able to accept Numpy arrays if you have boost::numpy installed. If you do not have boost::numpy it will simply be exposed as a scalar function.

Calculating a muon effective area

To calculate a muon effective area for single muons, simply sum up the inverse of the generated fluences for each event, e.g.:

mctree = frame['I3MCTree']
primary = mctree.primaries[0]
muon = mctree.get_daughters(primary)[0]
bundle = BundleConfiguration([BundleEntry(0, muon.energy)])
area_weight = 1./generator.generated_events(primary, bundle)

area_weight has units of \(GeV \, m^{2} \, sr\), and is analogous to NeutrinoGenerator’s OneWeight. To obtain the effective area in units of \(m^{2}\) as a function of muon energy and direction, fill area_weight into a histogram and divide each bin by its width in energy and solid angle.