Interface for flux and rate calculation

Note

Users should consider using high-level interface, provided by snewpy.snowglobes This low-level interface is not guaranteed to be stable and may change at any time without warning.

Flux containers

The module snewpy.flux defines the class snewpy.flux.Container - a container for the neutrino flux, fluence, event rates etc.

This object wraps a 3D array, and its dimensions: flavor, time and energy.

Usage

The flux container will be produced by the SupernovaModel.get_flux(), and consumed by the RateCalculator.run()

A example of usage:

import astropy.units as u
import numpy as np
from snewpy.models import ccsn
from snewpy.flavor_transformation import AdiabaticMSW

#get the suernova model
model = ccsn.Bollig_2016(progenitor_mass=27<<u.Msun)
#define the sampling points
energy = np.linspace(0,100,51)<<u.MeV
time = np.linspace(0,10,1000)<<u.s

#calculate the flux
flux = model.get_flux(time, energy, distance=10<<u.kpc, flavor_xform=AdiabaticMSW())

#optionally: integrate the flux over time bins to obtain the fluence
fluence = flux.integrate('time', limits=np.linspace(0,1,21)<<u.s)

#calculate the event rates using the RateCalculator
from snewpy.rate_calculator import RateCalculator
rc = RateCalculator()
event_rates = rc.run(fluence, detector='icecube')
ibd_rate = event_rates['ibd'] #will be also a Container instance
#get the total values vs. time bin
N_ibd_vs_T = ibd_rate.sum('energy')
#or total values vs. energy bin
N_ibd_vs_E = ibd_rate.sum('time')
#or total number of interactions
N_ibd_tot = ibd_rate.sum('time').sum('energy')
#to retrieve the number we can access the `array` member
print(N_ibd_tot.array.squeeze()) #312249.3525844854

Reference

class snewpy.flux.Container(data: Quantity, flavor: List[Flavor], time: Unit('s'), energy: Unit('MeV'), *, integrable_axes: Set[Axes] | None = None)[source]

A container class storing the physical quantity (flux, fluence, rate…), which depends on flavor, time and energy.

Parameters:

  • data (astropy.Quantity) – 3D array of the stored quantity, must have dimensions compatible with (flavor, time, energy)

  • flavor (list of snewpy.neutrino.Flavor) – array of flavors (should be len(flavor)==data.shape[0]

  • time (array of astropy.Quantity) – sampling points in time (then len(time)==data.shape[1]) or time bin edges (then len(time)==data.shape[1]+1)

  • energy (array of astropy.Quantity) – sampling points in energy (then len(energy)=data.shape[2]) or energy bin edges (then len(energy)=data.shape[2]+1)

  • integrable_axes (set of Axes or None) – List of axes which can be integrated. If None (default) this set will be derived from the axes shapes

can_integrate(axis)

return true if can be integrated along given axis

can_sum(axis)

return true if can be summed along given axis

integrate(axis: Axes | str, limits: ndarray | None = None) → Container

Integrate along given axis, producing a Container with the integral quantity.

Parameters:

  • axis (Axes or str) – An axis to integrate over. String should be one of "time" or "energy" (check Container.can_integrate())

  • limits (np.ndarray or None) – A sorted array (or astropy.Quantity consistent with the units of given axis) of integration limits If limits are None (default), then integrate over the the whole range of this axis Otherwise limits are treated as bin edges - the integration happens within each bin limits

Returns:

Container with integrated value

Raises:

ValueError – if the given axis cannot be integrated over (i.e. it must be summed instead, see Container.sum()).

Example

The resulting data array will be 3D array, but the dimension, corresponding to axis parameter will be reduced to len(limits)-1.

For an examplar Container a of a given shape:

>>> a.shape
(4, 10, 20)
>>> a.integrate('time').shape
(4, 1, 20)
>>> a.integrate('time', limits=[0, 0.5, 1]<<u.s).shape
(4, 2, 20)

The axis in the class will also be modified, keeping only the integration limits:

>>> a.time
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.] s
>>> a.integrate('time').time
[0., 9.] s
>>> a.integrate('time', limits=[0, 0.5, 1]<<u.s).time
[0., 0.5, 1.] s

All the other axes and dimensions will be kept the same

classmethod load(fname: str) → Container

Load container from a given file

save(fname: str) → None

Save container data to a given file (using numpy.savez)

sum(axis: Axes | str) → Container

Sum along given axis, producing a Container with the summary quantity.

Parameters:

axis (Axes or str) – An axis to sum over. String should be one of "flavor", "time" or "energy" (check Container.can_sum())

Returns:

Container with summed value

Raises:

ValueError – if the given axis cannot be summed over (i.e. it must be integrated instead, see Container.integrate()). One can check summable axes with Container._sumable_axes

Example

The resulting data array will be 3D array, but the dimension, corresponding to axis parameter will be reduced to 1.

For an examplar Container a of a given shape:

>>> a.shape
(4, 10, 20)
>>> a.sum('flavor').shape
(1, 10, 20)

The axis in the class will also be modified, keeping only the first and last points of the summation:

>>> a.flavor
array([0, 1, 2, 3])
>>> a.sum('flavor').flavor
array([0, 3])

All the other axes and dimensions will be kept the same

class snewpy.flux.Axes(value)[source]

Enum to keep the number number of the array dimension for each axis

Rate calcuation

The module snewpy.rate_calculator defines a Python interface for calculating event rates using data from SNOwGLoBES

Reference

class snewpy.rate_calculator.RateCalculator(base_dir='')[source]

Simple rate calculation interface. Computes expected rate for a detector using SNOwGLoBES data.

Use RateCalculator.run() method to run the simulation for specific detector and flux object.

Input flux can either be energy differential flux \(F_\nu(E_i) = \frac{d N}{d_E}(E_i)\) or integral flux in the energy bins \(N_{\nu}(E_i)\).

If the detector_effects=False the result is the number of neutrino interactions \(N\) - a product with with cross-section \(\sigma(E)\) and the number of targets in the detector \(N_{tgt}\):

\[N (E_i) = N_{tgt} \cdot F_\nu(E_i)\cdot \sigma(E_i)\]

If the detector effects are included, it is integrated within the energy bins and multiplied by the smearing matrix and detection efficiency to get number of interactions in energy bin:

\[N_i = N_{tgt} \cdot \int\limits_{E_i}^{E_{i+1}}~F_\nu(E_i)~\sigma(E_i)~d E\]

and this is convoluted with the detector smearing matrix \(M_{ij}\) and efficiency \(\varepsilon(E)\):

\[N^{rec}_i = \sum\limits_{j} M_{ij} \cdot N_j \cdot \varepsilon(E_i)\]
Parameters:

base_dir (Path or None) – Path to the directory where the cross-section, detector, and channel files are located If empty, try to get it from $SNOWGLOBES environment var

run(flux: Container, detector: str, material: str | None = None, detector_effects: bool = True) → Dict[str, Container][source]

Run the rate calculation for the given detector.

Parameters:

  • flux (Container) – The incoming neutrino flux (or fluence).

  • detector (str) – Name of the detector to calculate the rate. Check RateCalculator.detectors for the list of options

  • material (str or None) – Name of the detector material. Check RateCalculator.materials for the list of options If None (default) it will be guessed based on the name of the detector

  • detector_effects (bool) – If true (default), account for efficiency and smearing. If false, consider a perfect detector.

Returns:

dict[str, Container] – A dictionary with interaction rates (as instances of snewpy.flux.Container) for each channel.