Add "calculate neutrino flux at given distance" to SupernovaModel

[Sheshuk]

What I want

As a user I would like to get the supernova neutrino flux, predicted by a certain model (subclass of SupernovaModel) at a given distance at given times and energies (in neutrinos/cm^2/MeV/s).

What we have now

Currently SupernovaModel class allows only to get neutrino flux in 4pi (get_initial_spectra, get_transformed_spectra), and has no means for calculating flux at a given distance.

We have these calculations only in generate_fluence and generate_time_series: https://github.com/SNEWS2/snewpy/blob/b35379952fe18aedc28cdb106f1ee9e7ad08b082/python/snewpy/snowglobes.py#L122 https://github.com/SNEWS2/snewpy/blob/b35379952fe18aedc28cdb106f1ee9e7ad08b082/python/snewpy/snowglobes.py#L288

which will finally save the flux/fluence to the disk, which is not always needed.

While these calculations are not complicated, we shouldn't make user I think they belong to the SupernovaModel class.

What I suggest

1. Create a method in the base class SupernovaModel.get_initial_flux (t, E, distance, **kwargs) which by default would call self.get_initial_spectra. If there are any additional input parameters for a specific model (direction? flavors?), this method should be reimplemented in the models's class, making use of **kwargs argument.

2. This method should return an object, containing the flux for each of the (Flavor,T,E) points, and optionally pack

   • Energy points E
   • Time points T
   • Available flavors

   In a simplest case it could be Dict[Flavor, np.array], but I propose to discuss this separately in #214

[Sheshuk]

After that I suggest to use this new method in generate_fluence and generate_time_series

[JostMigenda]

I think for the initial flux/fluence it only really makes sense to get that at one distance (namely, the edge of the SN simulation boundary); since at larger distances from the star, the neutrinos have already undergone some FlavourTransformation(s). Therefore, I think a get_transformed_flux (with FlavourTransformation already applied) might be more appropriate?

Also, while the 1/(4 π d^2) factor is generic, I think the factors for time/energy bins are specific to SNOwGLoBES input files? So wouldn’t those need to remain in the generate_* functions?

[Sheshuk]

I think for the initial flux/fluence it only really makes sense to get that at one distance (namely, the edge of the SN simulation boundary); since at larger distances from the star, the neutrinos have already undergone some FlavourTransformation(s). Therefore, I think a get_transformed_flux (with FlavourTransformation already applied) might be more appropriate?

I agree from the physical point of view. However I suggest to keep the ability to apply the FlavorTransformation on the already existing flux - for example when checking various oscillation hypotheses without the need of calling the model calculation again and again, or for the case of having additional transformation layer, for example, Earth matter effect.

So how about a compromise?

class SupernovaModel:
    ...
    def get_flux (self, t, E, distance, flavor_transform=NoTransformation, **kwargs):
        ...

Also, while the 1/(4 π d^2) factor is generic, I think the factors for time/energy bins are specific to SNOwGLoBES input files? So wouldn’t those need to remain in the generate_* functions?

I totally agree, I only meant the distance factor here.

[JostMigenda]

That get_flux suggestion looks fine to me, yes. The implementation would then simply be

def get_flux (self, t, E, distance, flavor_transform=NoTransformation, **kwargs):
    return get_transformed_spectra(self, t, E, flavor_transform) / (4 * np.pi * (d*1000*3.086e+18)**2) 

, I assume? I’m not sure it’s worth adding another function just to multiply by this constant factor, but I’m not actively opposed to it, either.

[Sheshuk]

The implementation would then simply be

yes, something like that. However I suggest using astropy.units for conversion rather than hardcoded values. Also since now get_transformed_spectra returns dict, it's not as easily scalable, so it would be like this:

from astrpy import units as u

def get_flux (self, t, E, distance, flavor_transform=NoTransformation, **kwargs):
    distance = distance << u.kpc #assume that provided distance is in kpc, or convert
    factor = 1/(4*np.pi*(distance.to('cm'))**2)
    flux = get_transformed_spectra(self, t, E, flavor_transform)
    return {flavor: f*factor for flavor,f in flux.items()}

I’m not sure it’s worth adding another function just to multiply by this constant factor, but I’m not actively opposed to it, either.

Do you mean a function inside SupernovaModel? I don't think that would make sense. But scaling function would be useful in the NeutrinoFlux object proposed in #214. This would allow to evade awkward code like return {flavor: f*factor for flavor,f in flux.items()}