SNEWS2 / snewpy

A Python package for working with supernova neutrinos
https://snewpy.readthedocs.io
BSD 3-Clause "New" or "Revised" License
24 stars 17 forks source link

Detector and DetectionChannel classes for rate calculator #286

Closed Sheshuk closed 5 months ago

Sheshuk commented 8 months ago

This is a refactoring of the RateCalculator internals in a more OOP style. I define new classes Detector and DetectionChannel which contain the configuration info, read from SNOwGLoBES, and capable of calculating the event rates.

This is an internal change, which does not change the existing interface.

However it defines a lower level interface, which enables the following use case:

Use case

As a user I want to be able to run the event rate estimation on a custom detector, not the ones which are present in SNOwGLoBES. I want to be able to make a copy of the existing detector, but set different mass, efficiencies, etc, and run the rate calculation for this case. Or define a new detector with my own parameters (providing the detection channels etc)

Usage example

See Detector_demo.ipynb for more expanded example.

import astropy.units as u
from snewpy.models import ccsn

#calculate the neutrino flux
m = ccsn.Bollig_2016(progenitor_mass=11.2<<u.M_sun)
t = np.linspace(0,10,101)<<u.s
E = np.linspace(0,100,51)<<u.MeV
flux = m.get_flux(t, E, distance=10<<u.kpc)
#get the detector configuration
from snewpy.rate_calculator import RateCalculator
rc = RateCalculator()
#read the standard detector configuration
det  = rc.read_detector('icecube')

#make a modified copy of the detector configuration
from copy import copy
det1 = copy(det)
det1.mass = 10000<<u.kt #different mass
det1.channels['ibd'].efficiency=1. #100% efficiency for IBD

#function to plot the rates
def plot_rate(rate, *args, **kwargs):
    x = rate.energy.to_value('MeV')
    y = rate.integrate_or_sum('time').array.squeeze()
    if rate.can_integrate('energy'):
        plt.plot(x,y, *args, **kwargs)
    else:
        plt.stairs(y,x, *args, **kwargs)

#calculate and plot right away
plot_rate(det.run(flux)['ibd'], label='scint20kt default', ls='-',color='k')
plot_rate(det1.run(flux)['ibd'], label='scint10kt modified', ls='--',color='k')
plt.legend()
plt.show()

Outlook

In the future this can be a more convenient interface for the rate calculation: access and separate classes for Smearing and Efficiencies will allow to define them from the code more in a more flexible way.

TODO:

Allow defining smearing matrices for:

Sheshuk commented 7 months ago

I think I addressed all the issues, could you review again @mcolomermolla?