.. _event_weighting:

Weighting of genie-reader simulation
====================================

Because of technical reasons, simulation produced with ``genie-reader`` has to be weighted differently then ``GENIE-IceTray`` or NuGen simulations.
This document provides instructions on how to properly weight genie-reader events for analysis.

General idea of weighting IceCube simulation
--------------------------------------------

Usually it is convenient to generate simulation with an artificial (non-physical) spectrum and in a specific geometry (one of possible reasons could be to populate phase-space regions with low selection efficiency) and then reweight it.

Goal of weighting is then to "undo" event distributions resulting from generation process, such as generated flux and also geometrical effects, so

   a) we can combine sets of events generated with different geometries and/or energy ranges;

   b) can get realistic rate when weighting events with physical flux.

In case of typical IceCube neutrino simulation,

.. math::

   w = P_{int}\cdot S\cdot \Omega \cdot \frac{\Phi_{phys}(E,[\theta, ...])}{\Phi_{gen}(E)},

where :math:`P_{int}` is interaction probability, :math:`S` is the cross section of generated cylinder, :math:`\Omega` is a solid angle in which neutrinos were generated, :math:`\Phi_{phys}(E,[\theta, ...])` is physical flux and :math:`\Phi_{gen}(E)` is generated flux.

In IceCube simulation generated flux is generally a power law spectrum

.. math::

   \Phi_{gen}(E) \sim E^{-\gamma},

so

.. math::

   \Phi_{gen}(E) = A\cdot E^{-\gamma}.

One can see that amplitude :math:`A` is then needed to calculate event weight. The amplitude can be found if we know number of flux events (initially generated events which haven't necessarily interacted in the detector),

.. math::

   N_{flux~ev.} = \int_{E_{min}}^{E_{max}} \Phi_{gen}(E) dE = A\cdot\int_{E_{min}}^{E_max} E^{-\gamma} dE,

then amplitude is 

.. math::

   A = \frac{N_{flux~ev.}}{\int_{E_{min}}^{E_max} E^{-\gamma} dE}.

Event weight then can be calculated as 

.. math::

   w = P_{int}\cdot S\cdot \Omega \cdot \frac{\int_{E_{min}}^{E_max} E^{-\gamma} dE}{N_{flux~ev.}\cdot E^{-\gamma}}\cdot \Phi_{phys}(E,[\theta, ...]).

.. caution::

   This should be done for each of neutrino types and/or generation volume/energy range setups separately!

Number of events vs. Number of flux events
------------------------------------------

It is important to use **number of flux events** and **not number of interactions for weighting**, as they are not the same.

.. math::

   N_{flux~ev.} = \int_{E_{min}}^{E_{max}} \Phi_{gen}(E) dE,

.. math::

   N_{int.~ev.} = \int_{E_{min}}^{E_{max}} P_{int}(E)\cdot \Phi_{gen}(E) dE,

.. math::

   \Rightarrow N_{flux~ev.} \ne N_{int.~ev.}

Weighting genie-reader simulation
---------------------------------

In ``genie-reader`` simulation number of flux neutrinos :math:`N_{flux~ev.}` is stored inside an :cpp:class:`icetray::simclasses::I3GenieInfo` object in S-frame (typically named ``I3GenieInfo``) and can be obtained as 

.. code-block:: python

       n_flux_events = frame['I3GenieInfo'].n_flux_events

In order to correctly weight simulation user should sum up ``n_flux_events`` values for all used files **for a given neutrino type** (NuE, NuEBar, NuMu, NuMuBar, NuTau, NuTauBar) and **with the same energy (zenith, azimuth) range and generation volume** and use this number to divide the weight:

.. math::

   N_{flux~ev}^{total} = \sum_{files} N_{flux~ev.},

then

.. math::

   w = \frac{OneWeight}{N_{flux~ev}^{total}}\cdot \Phi_{phys}(E,[\theta, ...]),

where :math:`OneWeight = P_{int}\cdot S\cdot \Omega \cdot \frac{\int_{E_{min}}^{E_max} E^{-\gamma} dE}{E^{-\gamma}}`, as a standard for IceCube neutrino simulation, is precalculated and stored in ``I3MCWeightDict`` in Q-frames.

.. caution::

   Each set of simulation with specific settings (energy range, neutrino type, etc.) should be weighted using it's corresponding :math:`N_{flux~ev}^{total}` value!