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

bug with xent detector #318

Open mcolomermolla opened 2 months ago

mcolomermolla commented 2 months ago

running xent detector causes an error in snewpy when doing:

rT = r.sum('energy') ValueError: Cannot sum over energy! Valid axes are {<Axes.flavor: 0>}

I suspect this is because xent is missing smearing matrices (only detector I have run that does not have them), but I am not sure.

JostMigenda commented 2 months ago

Hi Marta, thanks for reporting this!

Just to make sure I understand correctly: You're running code from FluxContainer_demo.ipynb? That line looks like it's from the plot_rates function there; but that isn't actually run by the notebook currently? I think I'd need a bit more context of what code you're running, to make sure I can reproduce this issue.

mcolomermolla commented 2 months ago

Hi Jost, It is a code of my own (based indeed on the notebook you mention). I just do the usual steps to get the expected rates for the detector. And then integrate in energy to build the lightcurve. It does work for other detectors without any error. model = Bollig_2016(progenitor_mass=27<<u.Msun) times = np.linspace(0,10,10000)<<u.second; energies = np.linspace(0.5,100,200)<<u.MeV flux = model.get_flux(t = times, E = energies, distance=10<<u.kpc) from snewpy.flux import Container rc = RateCalculator(base_dir="/home/marta/snowglobes_git/") rates = rc.run(flux, detector='xent', material='xenon_coh') times=np.zeros(len(times)) for ch,r in rates.items(): rT = r.sum('energy') for i in range(len(rT.time)): times[i] = round(rT.time.value[i],3)

JostMigenda commented 2 months ago

Thanks for this code sample; that helps a lot! I’ve cleaned it up a bit (see below, including nicer formatting), but the only substantive change I’ve made is removing the base_dir argument from RateCalculator(), so it uses the files from the snowglobes_data package instead of your SNOwGLoBES repo.

With that change, the code runs without errors for me. snowglobes_data contains the exact same files as the SNOwGLoBES v1.3.1 release. Can you try to run the code below on your machine, so we can identify whether the issue is related to modifications in your local SNOwGLoBES repo or to something else?

import numpy as np
from astropy import units as u

from snewpy.flavor_transformation import AdiabaticMSW
from snewpy.models.ccsn import Bollig_2016
from snewpy.neutrino import MassHierarchy
from snewpy.rate_calculator import RateCalculator

model = Bollig_2016(progenitor_mass=27<<u.Msun)
times = np.linspace(0,10,10000)<<u.second
energies = np.linspace(0.5,100,200)<<u.MeV
flux = model.get_flux(t = times, E = energies, distance=10<<u.kpc, flavor_xform = AdiabaticMSW(mh=MassHierarchy.NORMAL))
rc = RateCalculator()  # No base_dir set, so I’m using the `snowglobes_data` package.
rates = rc.run(flux, detector='xent', material='xenon_coh')

times = np.zeros(len(times))
for ch, r in rates.items():
    rT = r.sum('energy')
    for i in range(len(rT.time)):
        times[i] = round(rT.time.value[i], 3)