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.